The use of GARCH I have a conceptual question that I haven't managed to grasp yet and is most likely a econometrics 101 question by here it goes:
If we estimate a GARCH model for a time series, how do we then use this in my model for the returns? For example; I have the return data of an index. I know that I have volatility clustering in this data. I find a suitable GARCH model for the volatility (variance). Now, if I model the returns an a suitable model, i.e. a regression model, and look at the coefficients and the p-values that it spits out, these values are still based on the regular OLS assumptions right? How do I make use of the GARCH in this model so that I can get coefficients and p-values that have accounted for the conditional heteroscedastic variances in the time series?
 A: Consider the data generating process (DGP):
$y_{t}= x_{t}^{\prime}\beta+\varepsilon_{t}$ (1)
$\varepsilon_{t}= \sigma_{t}z_{t},\quad z_{t}\sim i.i.d.\, N\left(0,\,1\right)$ (2)
$\sigma_{t}^{2}= \sigma^{2}+\alpha\varepsilon_{t-1}^{2}+\beta\sigma_{t-1}^{2}$ (3)
Equation (1) is a model for the conditional mean of the process Eq. (2) and (3) define a model for the conditional variance of the process (in this case the residuals are Gaussian). You could estimate these equations step-by-step using OLS (same as testing for ARCH(k) effect) since OLS is consistent but OLS will be inefficient and there will be non-linear estimators such as the Maximum Likelihood estimator (ML) which will produce a lower variance. The likelihood function for the model above would look like:
$l_{T}\left(\theta\right)=-\frac{T}{2}\log\left(2\pi\right)-\frac{1}{2}\sum_{t=1}^{T}\left(\log\sigma_{t}^{2}\left(\theta\right)+\frac{\varepsilon_{t}^{2}}{\sigma_{t}^{2}\left(\theta\right)}\right)\Leftrightarrow$
$l_{T}\left(\theta\right)=-\frac{T}{2}\log\left(2\pi\right)-\frac{1}{2}\sum_{t=1}^{T}\left(\log\left(\sigma^{2}+\alpha\varepsilon_{t-1}^{2}+\beta\sigma_{t-1}^{2}\right)+\frac{\left(y_{t}-x_{t}^{\prime}\beta\right)^{2}}{\sigma^{2}+\alpha\varepsilon_{t-1}+\beta\sigma_{t-1}^{2}}\right)$
You could drop the first term when maximizing since its constant. In practice you would start by estimating Eq. (1) and saving the residuals. If all misspecification tests are okay and if you have do have ARCH effect then you would estimate Eq. (3) (or a similar ARCH/GARCH family model) on your saved residuals. Furthermore you could calculate s.e.'s from the negaative of the expected value of the Hessian matrix although robust s.e.'s are recommended in case your distributional assumption does not hold (QMLE).
To illustrate I have simulated an AR-GARCH process:
$y_{t}=\mu+\theta y_{t-1}+\varepsilon_{t}$ (4)
$\varepsilon_{t}=\sigma_{t}z_{t},\quad z_{t}\sim i.i.d.\, N\left(0,\,1\right)$ (5)
$\sigma_{t}^{2}=\sigma^{2}+\alpha\varepsilon_{t-1}^{2}+\beta\sigma_{t-1}^{2}$ (6)
with parameter values $\mu=0.01
 $, $\theta=0.6
 $, $\sigma^{2}=0.04
 $, $\alpha=0.2
 $ and $\beta=0.5
 $. First I fit an AR(1) model to the simulated series:

I estimate the model using ML and get estimated values of: $y_{t}=0.0013538+0.59745y_{t-1}+\varepsilon_{t}$
Then I save the residuals an dplot the residuals and squared residuals. We see there are ARCH effects as indicated by Residual2 although I perhaps should have simulated a more persistent series.

Then I estimate a GARCH(1,1) model, Eq. (6) on the series called Residual. This gives us estimated values of: $\sigma_{t}^{2}=0.03720+0.23753\varepsilon_{t-1}^{2}+0.50754\sigma_{t-1}^{2}$
Which we can use to make forecasts of the conditional variance. Note that when dealing with financial series you will often end up with a model for the conditional mean corresponding with: $y_{t}=\mu+\varepsilon_{t}
 $, i.e. only a constant which correspond with the efficient market hypothesis since if it was possible to forecast the conditional mean of some givens share price it would be easy to make a profit.
Another way to estimate the model is to estimate in one go. Doing that I get a model: $y_{t}=0.003444+0.605931+\varepsilon_{t}$
$\varepsilon_{t}=\sigma_{t}z_{t}$
$\sigma_{t}^{2}=0.037121+0.237526\varepsilon_{t-1}^{2}+0.508122\sigma_{t-1}^{2}$
Which is very similar to the estimated model above which was estimated in two stages. Forecasting 10 periods into the future gives us the figure below.

We see that the first graph shows the simulated series, AR(1), and the fitted model while the second graph depicts the fitted values from the GARCH(1,1) model on the residuals from the conditional mean model. 
Usually you would first estimate the mean of the model and then when you would have well specified model for the conditional mean you would proceed to test for ARCH effects. The reason for this is that if your DGP follows an AR(2) process but you estimate an AR(1) model then your residuals will exhibit autocorrelation. If your residuals exhibit autocorrealtion this implies that your squared residuals will exhibit autocorrelation but in general the converse is not true. Therefore the ARCH test will also have power against residual autocorrelation and this is the reason to make sure that your model is well specified before testing for ARCH effects.
If this did not answer your question then let me know and I will amend my answer.
A: Suppose the conditional mean of returns is constant. A GARCH model gives you a fitted value of the conditional variance for each data point. These fitted values can be used to weight the data points to construct an efficient estimate of the mean (e.g. using weighted least squares); data points with high fitted conditional variance would be down-weighted relative to data points with low fitted conditional variance.
Now suppose the conditional mean of returns is not constant. Then you would build a model for the conditional mean simultaneously with a GARCH model for the conditional variance. The effect of the GARCH model would again be similar to the case discussed above. The data points with high fitted conditional variance would be down-weighted relative to the points that have low fitted conditional variance when estimating the model for the conditional mean.
One example given by @CadgasOzgenc is an ARMA-GARCH model. A rich choice of specifications of ARIMAX models and different versions of (G)ARCH models can be implemented using "rugarch" package in R (functions ugarchspec, ugarchfit).
Simultaneous estimation is efficient, but two-stage estimation could be done, too, if you can consistently estimate the conditional mean model in presence of conditionally heteroskedastic errors. First you would estimate the conditional mean model ignoring that the errors have a GARCH structure. Second, you would estimate a GARCH model on the residuals from the conditional mean model. Then you would reestimate the conditional mean model using the fitted conditional variances to weight the data points as discussed above. That could be done iteratively until convergence. For example, an AR-GARCH model could be estimated that way as an AR(p) model can be estimated consistently even in presence of GARCH errors. However, estimating an AR-GARCH model in one stage (simultaneously) would be more efficient.
