Consider the data generating process (DGP):
$y_{t}= x_{t}^{\prime}\beta+\varepsilon_{t}$ (1)
$\varepsilon_{t}= \sigma_{t}z_{t},\quad z_{t}\sim i.i.d.\, N\left(0,\,1\right)$ (2)
$\sigma_{t}^{2}= \sigma^{2}+\alpha\varepsilon_{t-1}^{2}+\beta\sigma_{t-1}^{2}$ (3)
Equation (1) is a model for the conditional mean of the process Eq. (2) and (3) define a model for the conditional variance of the process (in this case the residuals are Gaussian). You could estimate these equations step-by-step using OLS (same as testing for ARCH(k) effect) since OLS is consistent but OLS will be inefficient and there will be non-linear estimators such as the Maximum Likelihood estimator (ML) which will produce a lower variance. The likelihood function for the model above would look like:
$l_{T}\left(\theta\right)=-\frac{T}{2}\log\left(2\pi\right)-\frac{1}{2}\sum_{t=1}^{T}\left(\log\sigma_{t}^{2}\left(\theta\right)+\frac{\varepsilon_{t}^{2}}{\sigma_{t}^{2}\left(\theta\right)}\right)\Leftrightarrow$
$l_{T}\left(\theta\right)=-\frac{T}{2}\log\left(2\pi\right)-\frac{1}{2}\sum_{t=1}^{T}\left(\log\left(\sigma^{2}+\alpha\varepsilon_{t-1}^{2}+\beta\sigma_{t-1}^{2}\right)+\frac{\left(y_{t}-x_{t}^{\prime}\beta\right)^{2}}{\sigma^{2}+\alpha\varepsilon_{t-1}+\beta\sigma_{t-1}^{2}}\right)$
You could drop the first term when maximizing since its constant. In practice you would start by estimating Eq. (1) and saving the residuals. If all misspecification tests are okay and if you have do have ARCH effect then you would estimate Eq. (3) (or a similar ARCH/GARCH family model) on your saved residuals. Furthermore you could calculate s.e.'s from the negaative of the expected value of the Hessian matrix although robust s.e.'s are recommended in case your distributional assumption does not hold (QMLE).
To illustrate I have simulated an AR-GARCH process:
$y_{t}=\mu+\theta y_{t-1}+\varepsilon_{t}$ (4)
$\varepsilon_{t}=\sigma_{t}z_{t},\quad z_{t}\sim i.i.d.\, N\left(0,\,1\right)$ (5)
$\sigma_{t}^{2}=\sigma^{2}+\alpha\varepsilon_{t-1}^{2}+\beta\sigma_{t-1}^{2}$ (6)
with parameter values $\mu=0.01
$, $\theta=0.6
$, $\sigma^{2}=0.04
$, $\alpha=0.2
$ and $\beta=0.5
$. First I fit an AR(1) model to the simulated series:
I estimate the model using ML and get estimated values of: $y_{t}=0.0013538+0.59745y_{t-1}+\varepsilon_{t}$
Then I save the residuals an dplot the residuals and squared residuals. We see there are ARCH effects as indicated by Residual2 although I perhaps should have simulated a more persistent series.
Then I estimate a GARCH(1,1) model, Eq. (6) on the series called Residual. This gives us estimated values of: $\sigma_{t}^{2}=0.03720+0.23753\varepsilon_{t-1}^{2}+0.50754\sigma_{t-1}^{2}$
Which we can use to make forecasts of the conditional variance. Note that when dealing with financial series you will often end up with a model for the conditional mean corresponding with: $y_{t}=\mu+\varepsilon_{t}
$, i.e. only a constant which correspond with the efficient market hypothesis since if it was possible to forecast the conditional mean of some givens share price it would be easy to make a profit.
Another way to estimate the model is to estimate in one go. Doing that I get a model: $y_{t}=0.003444+0.605931+\varepsilon_{t}$
$\varepsilon_{t}=\sigma_{t}z_{t}$
$\sigma_{t}^{2}=0.037121+0.237526\varepsilon_{t-1}^{2}+0.508122\sigma_{t-1}^{2}$
Which is very similar to the estimated model above which was estimated in two stages. Forecasting 10 periods into the future gives us the figure below.
We see that the first graph shows the simulated series, AR(1), and the fitted model while the second graph depicts the fitted values from the GARCH(1,1) model on the residuals from the conditional mean model.
Usually you would first estimate the mean of the model and then when you would have well specified model for the conditional mean you would proceed to test for ARCH effects. The reason for this is that if your DGP follows an AR(2) process but you estimate an AR(1) model then your residuals will exhibit autocorrelation. If your residuals exhibit autocorrealtion this implies that your squared residuals will exhibit autocorrelation but in general the converse is not true. Therefore the ARCH test will also have power against residual autocorrelation and this is the reason to make sure that your model is well specified before testing for ARCH effects.
If this did not answer your question then let me know and I will amend my answer.