!  Maximum Likelihood estimation of a GARCH model.
Given a series of observations, we want to estimate the data generation process, 
allowing variance to vary over time  following a 
first order generalized autoregressive conditional heteroscedasticity model  
(i.e. GARCH(p,q), p=1 & q=1).   
The y(t) are assumed to be generated by the process:
    y(t) = alpha0 + alpha1*x(t) + a(t),
    a(t) = s(t)*e(t),
    s(t)^2 = beta0 + beta1*a(t-1)^2 + beta2*s(t-1)^2.
 where e(t) is white noise.                      
!  Reference: Hamilton, James D. (1994). Time Series Analysis,
 Princeton University Press (see Chapter 21).
Based on a What'sBest! model by Eren Ocakverdi;                     
!  Keywords:  ARIMA, Econometrics, Forecasting, GARCH, Heteroscedasticity,
   Maximum Likelihood, Time Series, Volatility Modeling, 
 ; 
! When solving, it is useful to use the multistart feature.
  Click on:
   Solver | Options... | Global Solver | Multistart solver Attempts | 4
   Solver | Options... | Global Solver | Use Global Solver | unchecked  ;
SETS:
  OBS: x, y, a2, s2, lns2, a2s2;
ENDSETS
DATA:
! Data could alternatively be retrieved from
a spreadsheet if the corresponding ranges in the spreadsheet
are named X and Y, and there is only one spreadsheet open in Excel,
and using the statements:
  ! y = @OLE( );
  ! x = @OLE( );
  Y             X =
-0.3708295  -1.4196
 0.2287753  -2.2396
 0.2081382  -0.4196
 0.0068829  -2.5446
 1.1795391  -1.2646
 3.1387306  -0.4846
 8.2863123   2.7354
 0.1658157  -1.5946
 5.3619457   2.5554
 3.8964039  -0.0046
-1.1897701  -6.1796
 5.2908945  0.4404
 8.2525741  0.1104
 6.8349661  0.0704
 4.2805493  -1.3146
10.4264812  -2.3296
 4.0687621  -2.2696
 1.8502611  -1.3046
-1.3458486  -3.5796
-1.8730369  -4.5696
-2.8381844  -5.5696
 5.3310167  8.2054
 0.0732557  -2.1846
-0.7577384  -1.3696
 4.4930833  7.2404
 2.2265120  3.0254
2.5421704  1.2904
2.9213433  2.3404
5.7213486  2.7404
3.6311287  2.4054
2.7817124  -3.1396
5.9180655  -1.1296
2.2575107  -1.9946
5.7674948  2.7154
9.2050620  1.5054
-0.0209955  -0.5696
5.6606594  6.1704
2.1783104  0.2604
1.5737642  1.2754
0.9454083  -0.5546
1.8044255  -0.5546
2.8915052  2.0454
2.5810686  -0.3846
0.7529759  -0.1596
0.6981101  1.1404
-0.5896699  -1.1146
0.3986052  0.2354
3.7521638  3.0854
1.9010561  -1.8446
1.4289591  -0.8096
5.4507511  -0.5696
7.5936498  0.2554
8.0108612  -0.6396
6.3355794  2.4804
4.2950874  -0.7446
2.1758734  -0.7196
3.3352518  0.0554
0.9927065  -1.9896
4.8748990  2.1504
3.8428777  0.3504
0.9342661  -2.8396
0.7488359  -0.5746
2.7820806  1.6504
2.8280793  4.0604
0.8892052  1.1104
0.2164499  -0.6196
1.8584283  2.8054
1.9861690  2.0654
1.7550617  -0.3346
-2.2903837  -5.4296
1.0867499  1.5404
1.4195743  0.8654
2.4263522  2.3254
0.0010931  -1.2346
3.0172286  3.3704
1.6453181  0.4004
1.2380105  -0.8246
4.1735987  2.4604
2.2087904  0.5204
3.4052056  3.1954
4.2784727  0.2404
3.2013125  -3.6396
8.4369970  3.0254
1.6932748  1.2504
2.2554868  3.5104
1.0315299  0.0404
-3.6637109  -6.1446
2.4601249  2.2954
3.8984784  2.8704
2.8347394  1.7754
2.9695648  3.3704
1.6975146  1.2054
1.1658735  1.3854
1.7135134  1.1904
4.7345067  5.3504
-0.8121917  -2.0596
-0.4192262  -2.2096
1.2585376  0.5704
1.9207934  2.4304
0.1402132  -0.4546
;
 ENDDATA

 NOBS = @SIZE( OBS);
 PI = 3.141592654;
 s2(1) = 5.295;
 ! Compute intermediate terms;
 ! Observation 1 is special;
   a2(1) = ( y(1)- alpha0- alpha1* x(1))^2;
   s2(1) = 5.295;
   
 @FOR( OBS(i) | i #GT# 1 :
   a2( i) = ( y( i)- alpha0- alpha1*x( i))^2;
   s2( i) = beta0 + beta1* a2( i-1) + beta2* s2( i-1);
     );

! Maximize the log likelihood function;
  MAX = -@SUM( OBS( i): @LOG( s2( i)))/2 
        -( NOBS*@LOG(2* PI)/2) 
        - @SUM( OBS(i): a2( i)/ s2( i))/2;

 ! Some apriori constraints on the parameters;
   alpha0 >= 0.2;
   alpha1 >= 0.1;
   beta0 >=  0.1;
   beta1 >= 0.2;
   beta2 >= 0.01;