(* This program provides the Yule-Walker type estimators for alpha and beta *)
(* concerning the  Active Data generated by a GARCH(1,1) process.           *)
(* An active sample is required.As alternative you can create a process by  *)
(* the function GenerateGarch.Then the main program has to be modified.     *)
(* Date: 21th of April 1998                                                 *)


program YuleWalkerEstimators;

var x:realvector;
    
   


(* 'GenerateGarch' Generates a GARCH(1,1) process with parameters gammaf, *)                                    
(* alphaf,betaf ,samplesize ssize and normal-distributed innovations      *)

function GenerateGarch(gammaf,alphaf,betaf:real;ssize:integer):realvector;
var s,halpha,hbeta:real;
    epsilon,sigma:realvector;
    i:integer;
begin
  while (alphaf+betaf)>=1.0 do begin
    writeln('Please enter values for alpha and beta so that alpha+beta<1');
    write('alpha= ');
    read(halpha);
    write('beta= ');
    read(hbeta);
    alphaf:=halpha;
    betaf:=hbeta;
  end;
  s:=gammaf/(1-(alphaf+betaf));
  epsilon:=normaldata(ssize);
  for i:=1 to ssize do begin
    sigma:=combine(sigma,s);
    s:=gammaf+s*(betaf+alphaf*sqr(epsilon[i]));
  end;
  return sigma*epsilon*epsilon;     (*  the squared process is returned!  *)  

end;




(* 'Gamma' calculates estimator for gamma *)

function Gamma(y:realvector):real;
var n:integer;
    z1,z2,z3,z4,num,denum:real;
begin
  n:=size(y);
  z1:=sum(y[2..(n-1)]*y[1..(n-2)])/n;
  z2:=sum(y[3..n]*y[1..(n-2)])/n;
  z3:=sum(y[1..(n-1)])/n;
  z4:=z3-y[n-1]/n;
  num:=z1-z2;
  denum:=z1/z3-z4;
  return num/denum;
end;



(* 'AlphaPlusBeta' calculates estimator for (alpha + beta) *)
(* fgamma:estimator for gamma                              *)

function AlphaPlusBeta(y:realvector;fgamma:real):real;
var n:integer;
    num,denum:real;

begin
 
  n:=size(y);
  num:=sum(y[2..n])/(n-1)-fgamma;
  denum:=sum(y[1..(n-1)])/(n-1);
  return num/denum

end;


(* 'OneDividedByBetaMinusAlpha' calculates estimator for 1/beta-alpha *)
(* h: estimator for (alpha+beta),fgamma: estimator for gamma          *)

function OneDividedByBetaMinusAlpha(y:realvector;fgamma,h:real):real; 
                                                
var h1,h2,num,denum:real;
    n:integer;
begin

  n:=size(y);
  h1:=sum(((y-fgamma)*y)[2..n]);
  h2:= sum((y-fgamma)[3..n]*y[1..(n-2)]);
  num:=(h1-h2)/(n-1);
  h1:=sum((y-fgamma)[2..n]*y[1..(n-1)]);
  h2:=sum((y*y)[1..(n-1)]);
  denum:=(h1-h2*h)/(n-1);
  return -num/denum

end;



(*  'Beta' calculates estimator for beta  *)

function Beta(v:realvector):real;
var c,z,hgamma:real;
    
begin 
  hgamma:=Gamma(v);
  z:=AlphaPlusBeta(v,hgamma);           
  c:=z+OneDividedByBetaMinusAlpha(v,hgamma,z);
  if c>=2.0 then
  return c/2.0-sqrt(sqr(c)/4.0-1.0)
  else return 0.0;
end;


 
(* 'Alpha' calculates estimator for alpha *)

function Alpha(v:realvector):real;
var z:real;
begin
  
  z:=AlphaPlusBeta(v,Gamma(v));
  return z-Beta(v);

end;



begin (* main program *)


x:=columndata(2)*columndata(2);

messagebox('The following parameters were estimated:' + chr(13) + chr (13) + 'gamma='+str(gamma(x))+',   alpha='+str(alpha(x))+',   '+'beta='+str(Beta(x)));


end.