Classical implicit residual type error estimators require using an underlying spatial finer mesh to compute bounds for some quantity of interest. Consequently, the bounds obtained are only guaranteed asymptotically that is with respect to the reference solution computed with the fine mesh. Exact bounds, that is bounds guaranteed with respect to the exact solution, are needed to properly certify the accuracy of the results, especially if the meshes are coarse. The paper introduces a procedure to compute strict upper and lower bounds of the error in linear functional outputs of parabolic problems. In this first part, the bounds account for the error associated with the spatial discretization. The error coming from the time marching scheme is therefore assumed to be negligible in front of the spatial error. The time discretization is performed using the discontinuous Galerkin method, both for the primal and adjoint problems. In the error estimation procedure, equilibrated fluxes at interelement edges are calculated using hybridization techniques.