Let $(X_n)_{n\in \mathbb Z}$ be a GARCH process with $E(X_0^4)<\infty$, and let $μ_n$ denote the distribution of $\frac 1{\sqrt n}\sum_{i=1}^n [X_i^2-\mathbb E(X_0^2)]$. We derive a numerical approximation of $μ_n$ when $x_1,...,x_n$ are observed. This yields the derivation of confidence intervals for $μ= E(X_0^2)$ and we investigate the accuracy of these confidence intervals in comparison with standard ones based on normal approximation. Moreover, when the innovation process has heavy tail distribution, we improve the method using a new resampling method.