Conditional Volatility of GARCH squared residuals Motivation
I want to wrap up my own GARCH implementation


*

*to make sure I have understood the underlying model/assumption.

*to leverage forecast::auto.arima to automate model order selection (see below).

*to mitigate the potential issues with parameters, so that I don't have to switch to IGARCH when the sum of parameters are close to 1, etc.


Background
Suppose $a_t$ is a residual term with zero mean and zero conditional mean, the GARCH(p, q) model is given by
$$
\begin{equation}
\left\{
\begin{aligned}
a_t &= \sigma_t \epsilon_t \\
\sigma_t^2 &= \alpha_0 + \sum_{i=1}^p \alpha_i a_{t-i}^2 + \sum_{j=1}^q \beta_j \sigma_{t-j}^2
\end{aligned}
\right.
\end{equation}
$$
Where $\epsilon_t \stackrel{iid}{\sim} \text{WN}(0, 1)$. If we let $\eta_t^2 := a_t^2 - \sigma_t^2$, then an alternative representation is
$$
\begin{equation}
\left\{
\begin{aligned}
a_t &= \sigma_t \epsilon_t \\
a_t^2 &= \alpha_0 + \sum_{i=1}^{\max(p, q)} (\alpha_i + \beta_i)a_{t-i}^2 + \eta_t - \sum_{j=1}^q \beta_j \eta_{t-j}
\end{aligned}
\right.
\end{equation}
$$
According to my professor, $\eta_t$ can be proved to have zero mean, and are uncorrelated (but not independent).
Methodology
From above we can see that the squared residual $a_t^2$ is an ARMA process, so as stated before, I'll just fit an ARMA on $a_t^2$. As far as I know, we don't have xxx::auto.garch but there is forecast::auto.arima, which would help us to find the correct order of ARMA.
This lets us model/forecast the conditional mean of squared residuals. However, what we really care about is the conditional volatility of $a_t$, but I'm lost here. 
$$
\begin{equation}
\begin{aligned}
E(a_t^2 | I_{t-1}) &= \text{var}(a_t | I_{t-1}) \\
&= \text{var}(\sigma_t \epsilon_t | I_{t-1}) \\
&= \text{var}(\sigma_t | I_{t-1}) \text{var}(\epsilon_t | I_{t-1}) + \text{var}(\sigma_t | I_{t-1}) E^2(\epsilon_t | I_{t-1}) + \text{var}(\epsilon_t | I_{t-1}) E^2(\sigma_t | I_{t-1}) \\
&= \text{var}(\sigma_t | I_{t-1}) + E^2(\sigma_t | I_{t-1}) \\
&= \text{???}
\end{aligned}
\end{equation}
$$
Question
How do I get "conditional volatility of residuals" from "conditional mean of squared residuals"?
 A: I just add this piece to the reasoning to complement what you have said and asked, as it seems that your initial objective is to use the ARMA autofit to implement a GARCH autofit. In order to do this, you need to model the process for squared returns, which is not easy due to the fact that the distribution of the squared returns depends on the true distribution of returns. If you assume that the latter is normal, as it usual under normal QMLE, the the distribution of the squared returns is parametric as highlighted later. Now let’s highlight an additional source of complexity: an additional reason why it is not easy to implement a GARCH autofit, and is not correct to use the ARMA autofit, is that the innovations in squared returns cannot be assumed to be iid, as they are not independent. Therefore it is not theoretically possible to transform the specification of a GARCH into the specification of an ARMA on squared returns due to the latter problem. Anyway, if you wish to do so by forcing an iid assumptions on the innovations of squared return to simplify the derivation of the likelihood function, then consider that you need the conditional variance of squared returns. So read the following. 
If you assume returns to be conditionally normal, then the conditional variance of squared returns is the variance of a squared normal distribution (as, if returns are conditionally normally distributed, here you are interested in the second central moment of the conditional distribution of their square). The formula for the calculation of the conditional variance of a squared normal with a given variance (which in your case is the given GARCH conditional variance, or the expected value of your squared ret) is given here along with more complete info on the shape of the conditional distribution of the squared returns. You will need to assume that the fourth central moment of return distribution exists finite, which is ok.
Now you have the conditional mean of squared returns, which is the GARCH conditional variance, the conditional variance of squared returns, and the pdf (from the link provided). Now you can write the likelihood for an ARMA on squared returns. But you can do iff you force the assumptions that the innovations for squared returns are iid. Which we cannot say it is right. So pay attention to the problem that your QMLE estimator may have.
A: Given that $\sigma_t^2$ is part of $I_{t-1}$ because of
$$
\sigma_t^2 = \alpha_0 + \sum_{i=1}^p \alpha_i a_{t-i}^2 + \sum_{j=1}^q \beta_j \sigma^2_{t-j}
$$
and given
$$
\epsilon_t \stackrel{iid}{\sim} \text{WN}(0, 1),
$$
couldn't you do
$$
\begin{equation}
\begin{aligned}
E(a_t^2 | I_{t-1}) &= \text{var}(a_t | I_{t-1}) \\
&= \text{var}(\sigma_t \epsilon_t | I_{t-1}) \\
&= \sigma_t^2 \text{var}(\epsilon_t | I_{t-1}) \\
&= \sigma_t^2 \cdot 1 \\
&=\alpha_0 + \sum_{i=1}^p \alpha_i a_{t-i}^2 + \sum_{j=1}^q \beta_j \sigma^2_{t-j}
\end{aligned}
\end{equation}
$$
?
