I have been working on manually implementing an ARMA GARCH (1,1) model but have been running into a few issues, namely a very large forecasted variance. I am hoping by outlining my process someone can catch a mistake.
First for ARMA, $$X_t = c + \phi X_{t-1} + \theta \epsilon_{t-1} + \epsilon_t$$ I assume some $c, \phi$, and $\theta$ and that $\epsilon_0 = 0$. Then by writing the above as $$\epsilon_t = X_t - c - \phi X_{t-1} -\theta\epsilon_{t-1}$$ I marched forward to construct a vector of $\epsilon_i$. The above process is optimized using OLS to return the set of $c, \phi$ and $\theta$, such that the vector of $\epsilon_i$ is minimized. I then use the optimized set of parameters to construct a forecast so that I can obtain a vector of ARMA residuals. Using my code, this is what an example output looks like:
Do these residuals seem reasonable? Is my assumption that $\epsilon_0 = 0$ fair, or should it be excluded afterwards? Also, I currently have no constraints on my ARMA parameters, some sources claim that $|\phi| \leq 1$ should be imposed. Will this make a difference?
I then feed these residuals into a GARCH model, $$\sigma_{t}^2 = \omega + \alpha\epsilon_{t-1}^2 + \beta\sigma_{t-1}^2\\\epsilon_i = \sigma_i p_i $$ where $p_i$ is a process associated with some (not necessarily Gaussian) distribution with mean 0 and scale 1.
To calculate the GARCH parameters I use MLE. I assume $\sigma^2_0 = 0$. For each time $n \in [0, t]$ I march forward using the GARCH equation and the ARMA residuals until I reach $\sigma_n^2$. I then evaluate $\frac{\epsilon_n}{\sigma_n}$ at the distribution's log PDF to obtain a probability. The sum of these values from $T=1$ to $T = t$ is my overall MLE. I also assume some parameters for the distribution (if there are any beyond the mean and scale). I then repeat this process so that I optimize both my distribution parameters and also my GARCH parameters. I require that all GARCH parameters be non-negative and that $\alpha + \beta = 1$
As noted, this process leads to incredibly large values for the forecasted variances (around 4000-5000). I notice it also chooses a rather large value for $\omega$ (somewhere around 8-30), and almost always chooses $\alpha = 0, \beta = 1$. For what it is worth, I am using the NLOPT package for the optimization and the COBYLA algorithm.
Does anyone see where I might have made a mistake or if I have made a false assumption? Thanks for reading!