Estimate GARCH parameters using maximum likelihood pseudocode I have to estimate the GARCH parameters using maximum likelihood in Scilab. I have tried many ways and so far nothing works properly. I have
$$ x_t = \sigma_t y_t, \ \ \ \ \ y_t \sim N(\mu, \sigma) $$
$$ \sigma^2_t = a_0 + a_1 x^2_{t-1} + b \sigma^2_{t-1} $$


*

*I know that the conditional distribution of of $x_t$ is normal. Do I have to find the unconditional distribution and how?

*Should I write in my code the joint likelihood function?

*How does the estimation process work? I have not found any article about that. First, we maximize the log likelihood function and find sigma at $t=0$. Second, we use our data and sigma, to iterate and find the other sigmas at $t= 1,2,3$ etc. Third, we use this to find the approximate estimates of $a_0$, $a_1$, $b$?

*It would be great if someone could provide me a link on how this process works and I could write it in my code.

 A: I can offer you my thoughts on how the general MLE procedure for GARCH parameters works, not the concrete implementation in SciLab, however. 
Assume $\sigma^2 = h_t$ and $\theta = (a_0,a_1,...,b_p)$. 
First, you need to determine the log-likelihood function $L(h_t)$. In the case of a normal distribution, it is pretty easy to derive. See here for example.
Second, you need to set its derivative to zero in order to find a maximum. This is more easily said than done. Then take the sample variance as initial value (as Richard Hardy suggested). 
Third, in order to do step 2 it is usually necessary to employ some kind of numerical optimization method. One widely used method is the quasi-Newton method. The gist is (without proof) to approximate the 'true values' with
$$\theta = \theta_t - \nabla L(h_t)(L''(h_t))^{-1}.$$
Fourth, as you can see, it is necessary to compute the Hessian matrix $L''(h_t)$ for each iteration. Since this is time-consuming it is standard procedure to turn to Hessian approximation (hence QUASI-Newton). Choose either the BFGS or BHHH algorithm. 
All in all, I am not familiar with SciLab, but I guess you need to implement some optimization method with approximations for the second derivative according to steps 1-4 I sketched above. In Stata, for example, these routines are already built in. Maybe you are able to draw from open source code somewhere. Otherwise you need to implement the procedures by yourself. For further reading, see e.g. here and here. 
