library cml;
  #include cml.ext;
  cmlset;

  dsn = "garch";

  depvar = { INDEX };
  indvar = { };  

/* set order here for a GARCH-M(p) model */
  p = 1;


output file = garchm.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$+"theta"$+ftocv(seqa(1,1,p),1,0)) | "gamma";
  else;
    z = z[.,dv]~ones(rows(z),1);
    numx = 1;
    _cml_ParNames = "CONST" | "kappa" |
        (0$+"theta"$+ftocv(seqa(1,1,p),1,0)) | "gamma";
  endif;



/*
**  constraints
**
**  kappa > 0.01, theta .>= 0, sumc(theta) < 1
**
*/


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

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


/*
**  Set weights for first p observations to zero
**  to keep them out of the log-likelihood
*/

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


/*
** start values
*/

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

  _cml_Options = { XPROD };

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


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

  __title = "LR confidence limits";
  cl1 = cmlpflclimits(b,f,z,0,&garchm);
  call cmlclprt(b,f,g,cl1,ret);


  _cml_BootFname = "gmboot";
  _cml_NumSample = 500;

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

  __title = "bootstrap confidence limits";
  clb = cmlblimits("gmboot");
  call cmlclprt(b,f,g,clb,ret);


  output off;



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

proc garchm( b, x );
  local u,beta,gamm,kappa,delta,theta,h,i,rx,cx;

  beta = b[1:numx];
  kappa = b[numx+1];
  theta = b[numx+2:numx+p+1];
  gamm = b[numx+p+2];
  
  rx = rows(x);
  cx = cols(x);

  u = zeros(rx,1);
  h = zeros(rx,1);
  
  h[1:p] = meanc((x[.,1] - x[.,2:cols(x)] * beta)^2) * ones(p,1);
  u[1:p] = h[1:p];

  i = p + 1;
  do until i > rx;
      h[i] = kappa + h[i-p:i-1]'*theta;
      u[i] = x[i,1] - x[i,2:cx]' * beta - gamm * h[i];
      i = i + 1;
  endo;

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