Today, numerical simulation has become a common tool for various engineering activities, such as design and optimization. Nevertheless, standard numerical methods can prove limited when it comes to simulating complex multi-dimensional models in real time. In this context, model reduction methods have been widely developed in recent decades, providing an effective solution when fast predictions are required. If the efficiency of simulation techniques is an important objective, so is the reliability of the results obtained. The Proper Generalized Decomposition (PGD) is a model reduction method based on a modal representation of the solution with separation of variables. It has proven its efficiency in many applications. The approximation error in PGD comes mainly from two sources: the truncation of the modal representation and the discretization error linked to the underlying numerical method used to compute modes. In order to assess the accuracy of PGD solutions, some a posteriori error estimation tools have been developed, using extensions of classical verification procedures in finite element analysis. These tools lead to error evaluation strategies for parametric PGD solutions or PGD solutions based on the separation of space and time. In this work, applications where space variables are separated one from another in the PGD decomposition are considered. This strategy is particularly relevant for applications related to plate or shell geometries. In this context, a first error indicator based on recovery techniques has been proposed. Only the discretization error is estimated for the purpose of mesh adaptation. The present work introduces an a posteriori error estimator, taking all error sources into account, with guaranteed bounds for the verification of PGD reduced models with separation of space variables. We restrict ourselves to 2D or 3D linear elliptic problems (such as the diffusion or the elasticity problems) on rectangular or parallelepiped domains. The error bounds are derived from the Constitutive Relation Error (CRE) method. The main difficulty of this approach lies in the construction of equilibrated fluxes (or stresses), and we show here how they can be obtained. The estimator and indicators that we introduce make it possible to certify the reduced model as well as to drive an adaptive strategy. We provide several numerical examples illustrating the efficiency of our procedure.