Numerical methods for parabolic homogenization problems combining finite element methods (FEMs) in space with Runge-Kutta methods in time are proposed. The space discretization is based on the coupling of macro and micro finite element methods following the framework of the Heterogeneous Multiscale Method (HMM). We present a fully discrete analysis in both space and time. Our analysis relies on new (optimal) error bounds in the norms L-2(H-1), C-0(L-2), and C-0(H-1) for the fully discrete analysis in space. These bounds can then be used to derive fully discrete spacetime error estimates for a variety of Runge-Kutta methods, including implicit methods (e.g. Radau methods) and explicit stabilized method (e.g. Chebyshev methods). Numerical experiments confirm our theoretical convergence rates and illustrate the performance of the methods.