Is a GARCH process serial-uncorrelated? A GARCH process is defined as

When it is wide-sense stationary, I wonder if $Cov(a_t, a_{t+h})=0$ for all $h>0$?
 A: The GARCH process that you have included in your question implicitly assumes no mean dynamics. Let me rewrite it:
$$
Y_t = \sigma_t \varepsilon_t 
$$
With some manipulations, we can derive that the ACF (recall that this is a mean zero process)
$$
\gamma_Y(k) = \mathbb{E}(Y_t Y_{t-k})
$$
is zero for $k>0$
$$
\begin{align}
\mathbb{E}(Y_t Y_{t-k}) &= \mathbb{E}(\mathbb{E}(Y_t Y_{t-k}\mid \mathcal{I_{t-k+1}})) \\
&= \mathbb{E}(Y_{t-k}\mathbb{E}(Y_t \mid \mathcal{I_{t-k+1}})) \\
&= \mathbb{E}(Y_{t-k}\mathbb{E}(\mathbb{E}(Y_t\mid \mathcal{I}_{t}) \mid \mathcal{I_{t-k+1}})) \\
&= \mathbb{E}(Y_{t-k}\mathbb{E}(\sigma_t\underbrace{\mathbb{E}(\varepsilon_t\mid \mathcal{I}_{t})}_{=0} \mid \mathcal{I_{t-k+1}})) \\
&= 0\qquad\forall\; k> 0 
\end{align}
$$
where all the results follow from the law of iterated expectations, and $\mathcal{I}_t\equiv \sigma(Y_{t-1}, \ldots)$ is the information set at time $t$.
If however, you have any mean dynamics included, the shape fo the ACF is the same as it would be in the case of the mean model without the GARCH innovations. 
$$
\phi(L)\Delta^dY_t =  \theta(L)( \sigma_t \varepsilon_t )
$$
where such autocorrelation is driven purely by the mean dynamics.
The key difference however is that inference about the ACF values is impacted when the innovations are heteroskedastic as in the case of GARCH(1, 1) innovations. This is because the variance of the ACF is no longer just a function of the sample size. This is covered in Mikosch and Starica (2000).  
You can see this using some simulations.
library(fGarch)
library(boot)

set.seed(1234)

#================================================
# simulate the rejection from a GARCH process
#================================================
fnGarchACF = function() {
  gspecSO = garchSpec(model = list(omega = 1e4,
                       alpha = c(0.1), beta = c(0.8)),
          cond.dist = 'norm')

  gsimSO = garchSim(gspecSO, n = 1000)
  invisible(abs(acf(gsimSO, plot = FALSE)$acf) > 1.96/sqrt(1000))
}

#================================================
# simulate the rejection from a white noise process
#================================================
fnWNACF = function() {
  abs(acf(rnorm(1000), plot = FALSE)$acf) > 1.96/sqrt(1000)
}

#================================================
# GARCH ACF rejections
#================================================
mResultsGarch = do.call(cbind, lapply(1:1000, function(x) fnGarchACF()))

# plot the rejection rates
plot(apply(mResultsGarch, 1, mean), main = 'Rejection rates for GARCH(1, 1) white noise',
     xlab = 'Lag', ylab = 'Rejection Rate', pch = 16, col = 'red')
abline(h = 0.05, col = 'blue')

#================================================
# WN ACF rejections
#================================================
mResultsWN = do.call(cbind, lapply(1:1000, function(x) fnWNACF()))

# plot the rejection rates
plot(apply(mResultsWN, 1, mean), main = 'Rejection rates for independent white noise',
     xlab = 'Lag', ylab = 'Rejection Rate', pch = 16, col = 'red')
abline(h = 0.05, col = 'blue')

Note how the null of white noise is systematically over-rejected in the case of GARCH innovations, compared to the rejection rate in the case of an independent white noise process.
The conclusion is that you should be careful in testing serial correlation in time series with variance dynamics such as the GARCH process using standard tests and visual diagnostics.


In order to test serial correlation in the presence of variance dynamics, such as GARCH, you can follow the conclusions in Silvapulle & Evans (1998).
