library cml;
  #include cml.ext;
  cmlset;

  dsn = "garch";

  depvar = { INDEX };
  indvar = { };  

/* set order here for a GARCH(p,q) model */
  p = 2;
  q = 2;


output file = garch.out reset;

/*
**  The dataset is read in using loadd because the GARCH log-likelihood
**  requires the entire data set and passing the data in as a matrix
**  will ensure that the complete data set will be passed to the
**  log-likelihood proc.
*/

  z = loadd(dsn);

  if not scalmiss(indvar);
    { depvar,dv,indvar,iv } = indices2(dsn,depvar,indvar);
  else;
    { depvar,dv } = indices(dsn,depvar);
  endif;

  if scalmiss(dv);
     errorlog "error: variable not found in dataset "$+dsn;
     end;
  endif;

  if not scalmiss(indvar);
    z = z[.,dv]~ones(rows(z),1)~z[.,iv];
    numx = rows(indvar) + 1;
    _cml_ParNames = "CONST" | indvar | "kappa" |
        (0$+"delta"$+ftocv(seqa(1,1,p),1,0)) |
        (0$+"alpha"$+ftocv(seqa(1,1,q),1,0));
  else;
    z = z[.,dv]~ones(rows(z),1);
    numx = 1;
    _cml_ParNames = "CONST" | "kappa" |
        (0$+"delta"$+ftocv(seqa(1,1,p),1,0)) |
        (0$+"alpha"$+ftocv(seqa(1,1,q),1,0));
  endif;



/*
**  constraints
**
**  kappa > 0.01, delta .>= 0, alpha .>= 0, sumc(delta|alpha) < 1
**
*/

  _ww_ = { -1e250 1e250 };
  _cml_Bounds = ones(numx+1+p+q,2).*_ww_;
  _cml_Bounds[numx+1,1] = .01;
  _cml_Bounds[numx+2:numx+p+q+1,1] = zeros(p+q,1);


  _cml_c = zeros(1,rows(_cml_Parnames));
  _cml_c[1,numx+2:numx+p+q+1] = -ones(1,p+q);
  _cml_d = -.99;

/*
**  Constraints for IGARCH model, uncomment these lines and
**  comment out the lines above for _cml_c and _cml_d
**
**  _cml_a = zeros(1,rows(_cml_Parnames));
**  _cml_a[1,numx+2:numx+p+q+1] = ones(1,p+q);
**  _cml_b = 1;
*/

/*
**  Set weights for first maxc(p|q) observations to zero
**  to keep them out of the log-likelihood
*/

  r = maxc(p|q);
  __weight = (rows(z)/(rows(z)-r))*ones(rows(z),1);
  __weight[1:r] = zeros(r,1);


/*
** start values
*/

  if not scalmiss(indvar);
    start = z[.,1] / z[.,2:cols(z)];
  else;
    start = meanc(z[.,1]);
  endif;
  start = start | 1 | .01 * ones(p+q,1);

  _cml_GradProc = &grd;

  { b,f,g,h,ret } = cml( z, 0, &garch, start );

  __title = "Wald confidence limits";
  cl1 = cmlclimits(b,h);
  call cmlclprt(b,f,g,cl1,ret);


  output off;



/*************** procedures *******************/

proc garch( b, x );
  local u2,kappa,delta,alpha,h;

  u2 = (x[.,1] - x[.,2:cols(x)] * b[1:numx])^2;

  kappa = b[numx+1];
  delta = b[numx+2:numx+p+1];
  alpha = b[numx+p+2:numx+p+q+1];

  h = recserar(kappa+_shape(u2,q,1)*alpha,meanc(u2)*ones(p,1),delta);

  retp( -0.5*( (u2 ./ h) + ln(2 * pi) + ln(h) ) );
endp;


proc grd(b,x);
  local j,u,u2,kappa,delta,alpha,r,h,s,g;

  u = x[.,1] - x[.,2:cols(x)] * b[1:numx];
  u2 = u .* u;

  kappa = b[numx+1];
  delta = b[numx+2:numx+p+1];
  alpha = b[numx+p+2:numx+p+q+1];

  g = ones(rows(x),rows(b));
  s = _shape(u2,q,1);
  h = recserar(kappa+s*alpha,meanc(u2)*ones(p,1),delta);

  g[.,1:numx] = -2*_dshape(u.*x[.,2:cols(x)],q,alpha);
  g[.,numx+2:numx+p+1] = _shape(h,p,-1);
  g[.,numx+p+2:numx+p+q+1] =  s;

  g = (0.5*(u2./h - 1)./h) .*
    recserar(g,zeros(p,cols(g)),delta.*ones(1,cols(g)));

  g[.,1:numx] = g[.,1:numx] + (u./h).*x[.,2:cols(x)];

  retp(g);
endp;


proc _shape(z,m,r);
   local y,n,v;
   n = rows(z) - m - 1;

   y = zeros(rows(z),m);
   if r == 1;
       v = seqa(1,1,m)' + seqa(0,1,n+1);
   else;
       v = seqa(m,-1,m)' + seqa(0,1,n+1);
   endif;

   y[m+1:rows(z),.] = reshape(z[v],n+1,m);

   retp(y);
endp;


proc _dshape(z,m,a);
   local y,j,n,v;
   n = rows(z) - m - 1;
   y = zeros(rows(z),cols(z));
   v = seqa(1,1,m)' + seqa(0,1,n+1);

   j = 1;
   do until j > cols(z);
       y[m+1:rows(z),j] = reshape(z[v,j],n+1,m) * a;
       j = j + 1;
   endo;

   retp(y);
endp;