Is there a way to use a GARCH model to prewhiten data (the way you can use ARIMA to prewhiten data)? Is there a way to do it in MATLAB (MATLAB gives two predictions epsilon and sigma, but I'm not sure which one to use to calculate the residuals)?Thanks

I am trying to understand what you're trying to do here. If you're trying to extract the white noise shocks (that is, the residuals) that are driving your model, then you can certainly do this for an ARMA-GARCH model. For an AR(1) model: y_t = B*y_t-1 + e_t, it's just e_t = y_t - B*y_t-1If you have an AR(1)-GARCH(1,1) model, then you do this with e_t = (y_t - B*y_t-1)/sigma_t where sigma_t is the fitted conditional standard deviation coming from the GARCH component of the model.Matlab can estimate an ARMA-GARCH model parameters jointly. It has a whole bunch of functions in the econometrics toolbox to work with these models (garchfit, garchsimulate) etc all of which have dozens of options.. I am 90% sure there is one that can extract the residuals. But even if there is no, it should be trivial to extract them yourself once you have the fitted parameters.

Thanks a lot APS, that's really helpful. I also had one other question. I haven't been able to found anything about how to choose the appropriate values for P, Q. I've seen that there are tests to determine if a ARCH process exists, but I couldn't find anything about what type of process (values for lags) exist.And finally, if you know of any way to do this in R, that would be super helpful. I actually don't have access to the Econometrics toolbox in MATLAB. Thanks again for your help.

A lot of the time people just estimate a GARCH(1,1) model and ignore other possibilities, perhaps following the reasoning that a very simple but possibly misspecified model is better than a complex misspecified model to avoid overfitting. There are probably some formal papers on p,q selection in GARCH models, but even the original GARCH paper (Bollerslev, 1986) has some good model hints in it. For one, GARCH(p,q) model is equivalent to an ARMA(max{p,q},p) model in e^2 such that the autoregressive coefficients of the GARCH model become corresponding MA coefficients in the ARMA model (that's why the number of those lags is the same equal to p). Therefore, I believe, you can apply all the model selection methodology for ARMA processes to the squares of variable you're studying. In particular, you can use the autocorrelation function of e^2 to select the value of p, and the partial autocorrelation function to select the value of max(p,q). You can apply the AIC or BIC criteria to select the model. You can also do something like a formal likelihood ratio test to test something like GARCH(p,q) against GARCH(p+1,q). Finally, you can also apply the tests for the ARCH effects to the residuals of the model. If the residuals look like white noise with no arch effects, then you're done. Obviously, this is a lot of things to test for. It would be nice if R did it, but I haven't looked into it myself. If you can't do it, just estimate GARCH(1,1) based on my own GARCH model selection philosophy that states roughly:1) Any ARCH model is better than no ARCH model.2) A complex ARCH model is not necessarily better than a simple GARCH model like GARCH(1,1) or GJR-GARCH(1,1) or EGARCH(1,1)Generally, testing and accounting for asymmetric effects in volatility is more important than messing around with lag selection, so I'd suggest to look into the decision between GJR-GARCH(1,1) (which is GARCH(1,1) with asymmetric effect of shocks depending on the sign) vs plain GARCH(1,1) first.