What's the point of (G)ARCH when you can square the residual and use ARMA? I'm taught that
$$
\begin{equation}
\begin{aligned}
X_t \sim \text{ARCH}(p) & \rightarrow X_t^2 \sim \text{AR}(p) \\
X_t \sim \text{GARCH}(p, q) & \rightarrow X_t^2 \sim \text{ARMA}(\max(p, q), p) \\
\end{aligned}
\end{equation}
$$
When I got a time series whose residual appears to be white noise, instead of fitting GARCH on it, I can simply


*

*Square the residual.

*Fit ARMA on the squared residual.

*Do some prediction, or inference, or whatever with ARMA.

*"Restore" the residual by taking the square root.


Why do I need a (G)ARCH model then?
 A: It is correct that, we can obtain an AR(p) representation for $X_t^2$ if $X_t$ follows an ARCH(p) process and an ARMA(max(p,q),p) representation for $X_t^2$ if $X_t$ follows a GARCH(p,q) process (see this questions - the generalization to the GARCH case is obvious). 
The reason why we are using a GARCH type of model is that we want to model the conditional volatility and not interested in modelling the squared observations. We take the squared observations to be a "realized measure" or signal of volatility. 
Assume that we are interested in forecasting, then for the AR(1) model, we have
$$
\sigma_{t+1}^2  = E_t[x_{t+1}^2] = w + \alpha E_t[x_t^2] + E_t[x_{t+1}^2 - \sigma^2_{t+1}] = w + \alpha x_t^2
$$
Thus, the forecasting equation corresponds to ARCH recursion. 
Typically, one sets up the (quasi)-loglikelihood of the GARCH model and maximizes it to get parameter estimates. Your "strategy" suggests a potential way of estimating ARCH type models described in Chapter 6 of Francq and Zakoian's book "GARCH Models: Structure, Statistical Inference and Financial Applications" (2010). 
To sum up. The simple answer to your question is that we are using GARCH type models because we want to model the conditional volatility!
A: I see two reasons for preferring the GARCH representation to its equivalent ARMA representation.


*

*You typically want to model the conditional distribution or certain distributional characteristics (such as the conditional mean $\mu_t$ and the conditional variance $\sigma^2_t$) of $X_t$ rather than $X^2_t$. A GARCH model does that directly by specifying a conditional mean equation for $\mu_t$, a conditional variance equation for $\sigma^2_t$, and the distribution of the standardized innovations. If your central interest is in the conditional variance of the process, it is more convenient to use the standard GARCH representation.
From GARCH, you get a specific equation for the conditional variance $\sigma^2_t$ directly (both the point estimates and the standard errors), while from the equivalent ARMA representation you have to back it up to get it expressed in terms of the objects of central interest ($\sigma^2_t$). It should be easy to back up the point estimates, but obtaining the relevant standard errors might be more challenging.
 (Thanks to Johan Stax Jakobsen for some elaboration on the topic in the comments to his own answer.) 

*If your process has a nonconstant conditional mean in addition to a nonconstant GARCH-type conditional variance, then you can estimate the model (both the conditional mean part and the conditional variance part) more efficiently if you do it simultaneously rather than doing it in a stepwise fashion (estimating the conditional mean model first, obtaining residuals, and then estimating the conditional variance model). 

