Specifying an ARMA-GARCH model without rugarch There have been probably hundreds of these posts, but I have yet to see any sort of answer to how to actually go about this. I need to fit a ARMA-GARCH model to my data because I cannot assume my data have a constant mean.
I cannot use rugarch as most of the posts addressing this have used so I have to be able to do this myself using Python and PyFlux.
Here are the posts I have looked at:


*

*ARMA-GARCH model selection / fit evaluation

*Fitting an ARCH/GARCH mode basics

*ARMA GARCH estimation in practice
My initial assumption was that I could fit an ARIMA model to the data and then fit a GARCH model to the residuals. This seems to be a common thing people do, but according to this post by @RichardHardy it is incorrect and inconsistent - though I can't find any evidence other than his posts indicating this. Unfortunately this has left me with more questions than answers because I can't just fall back to rugarch to solve this. I also haven't been able to find any literature on the process of fitting an ARMA-GARCH model that doesn't rely on handwaving rugarch as the solution. I'm hoping (perhaps Richard Hardy or someone who understands this stuff) could answer these questions:


*

*What is meant by "fitting simultaneously"? To me this feels like a chicken and egg problem.

*Would it be possible to fit the ARIMA model and then a GARCH model on the residuals if you had an adjusted Lijung-Box test and ACF/PACFs?

*Once the model is fitting, how is a forecast made from the the two models?
 A: You touch upon two main issues: estimation and model selection.
Estimation 
For a given model specification, you may 


*

*Either write down the likelihood function and feed it into a generic optimizer (such as the function optim in R); 

*Or use an existing function that takes the model specification (e.g. ARMA(p,q)-GARCH(s,r)), "writes the likelihood" for you and optimizes it (such as the function ugarchfit in the "rugarch" package in R).


Both ways are fine: 


*

*The first one allows you to dig deep into model fitting and perhaps alter model specification beyond what is available in the standard functions; however, it is time consuming and the way you are fitting might not be optimized for speed (while the standard functions typically are), unless you are a knowledgeable programmer with experience in optimization;

*The second one allows for a convenient and fast estimation without the need to understand how exactly it is done (the latter can also be a drawback depending on what you wish to achieve).


The likelihood functions of ARMA and GARCH are available in time series textbooks such as Hamilton "Time Series Analysis" (2004) or Tsay "Analysis of Financial Time Series" (2005), among other. Hopefully, the likelihood of ARMA-GARCH is also available somewhere, but I do not have a refence handy. You could try textbooks, lecture notes or maybe software documentation (R, Stata, Matlab). (Please post any references here in the comments, I will appreciate it.)
Estimation can also be done 


*

*either simultaneously (maximizing the likelihood function that reflects the dynamics of both the conditional mean and the conditional variance) 

*or in two steps (first estimate the conditional mean equation, implicitly assuming constant conditional variance, then estimate the conditional variance equation for a given estimated conditional mean). 


The problem with two-step estimation is that the first step uses an assumption that is violated in the second step and as such makes the estimators of both steps inefficient and sometimes inconsistent, as discussed in earlier posts (cited e.g. in the OP).
Model selection
Selecting a conditional mean-conditional variance model is not easy because it is a large animal and there are so many choices to make. Selecting sequantially (first the conditional mean model, then the conditional variance model) is suboptimal because the first step depends on assumptions that are violated in the second step, and I am not aware of any theoretical results that guarantee consistent selection or the like, as also discussed in previous posts. Nevertheless, this is sometimes done in practice and even recommended in time series textbooks such as Tsay "Analysis of Financial Time Series" (2005). However, I perceive the recommendation as "a" model selection strategy that is relatively easy, but not necessarily "the best" one.
Among other strategies probably the simplest, although computationally demanding one would be to fix a pool of candidate models (e.g. all submodels of ARMA(4,4)-GARCH(2,2)), estimate them (preferably simultaneously) and select the one with the lowest AIC (if the goal is forecasting) or BIC (if the goal is recovering the "true" model).
Questions 1, 2 and 3


*

*This should be clear from the text above.

*It is possible, but that does not mean it is desirable. As I mentioned, I am not aware of any theoretical results that guarantee consistent selection in the two-step approach. On the other hand, sometimes it makes sense to go for speed and convenience instead of a statistically "clean" approach, and this is what people often do.

*The model estimated on data until time $t$ will immediately yield a forecast for time $t+1$, because the left hand side of the model leads the right hand side by one time period. If you want to forecast further ahead, do that iteratively: consider the $h-1$-step-ahead forecast as real data and this way "extend your sample" by $h-1$ time points, then do exactly the same as when forecasting 1 step ahead (do not refit the model).

