@       PROC GARCH                          @
@   This proc will estimate a GARCH(p,q) model  @
@   for the variance and a user-specified model @
@   for the mean of a series.                   @
@                                               @
@   INPUT: Tx1 dependent variable               @
@     TxK independent variables                 @
@     number of lags of resids^2 in GARCH model @
@     number of lags of conditional variance    @
@     limits: 1=calculate Wald CL's 0 = do not  @
@     startvalues = vector of inital values     @
@   OUTPUT:                                     @
@     b, a (K+p+q+1)x1 vector of parameters     @
@     f, scalar, function at minimum            @
@     g, (K+p+q+1)x1 vector, gradient at min    @
@     c11, Wald confidence limits               @
@     ret, scalar, CML return code              @
@     hhat, Tx1 vector of conditional variance  @
@           estimates at each point in time     @

/*  GARCH Code found at:
    http://faculty.washington.edu/rons/garch.html
    A site connected to Ronald Schoenberg's site

    Modified by Andrew Patton
    email: apatton@econ.ucsd.edu
    Fri, Sep 3, 1999.

    This code is provided gratis, without performance warrantee.
    */



library cml;
cmlset;

__output=0;

_p=0;_q=0;_numx=0;  @ setting these up for future reference @
proc (6) =  garch(depvar,indepvar,garchp,garchq,limits,startvalues);
    local p,q,z,numx,r,start,b,f,g,vcv,ret,cl1,_ww_, y, _numx, _p, _q, u2, hhat, u;

    p = garchp;
    q = garchq;
    y = varput(p,"_p"); @ variable made global so it can be used in other proc's @
    y = varput(q,"_q"); @ variable made global so it can be used in other proc's @


/*  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.                                            */

    if indepvar/=0;
        z = depvar~indepvar;
        numx = cols(indepvar);
        _cml_ParNames = ("indep"$+ftocv(seqa(1,1,cols(indepvar)),1,0))  | "kappa" |
            (0$+"delta"$+ftocv(seqa(1,1,p),1,0)) |
           (0$+"alpha"$+ftocv(seqa(1,1,q),1,0));
    else;
        z = depvar;
        numx = 0;
        _cml_ParNames = "omega" |
            (0$+"beta"$+ftocv(seqa(1,1,p),1,0)) |
            (0$+"alpha"$+ftocv(seqa(1,1,q),1,0));
    endif;
    y = varput(numx,"_numx"); @ variable made global so it can be used in other proc's @

/*  constraints
**  kappa > 0.0001, 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] = .0001;
    _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 = -.9999;

/*  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 indepvar/=0;
       start = z[.,1] / z[.,2:cols(z)];
       start = start | 1 | .01 * ones(p+q,1);
    else;
        start = startvalues;
    endif;

    _cml_GradProc = &grd;
    { b,f,g,vcv,ret } = cml( z, 0, &garchll, start );

    if indepvar/=0;
        u = (depvar - indepvar * b[1:cols(indepvar)]);
        _numx = cols(indepvar);
    else;
        u = (depvar);
        _numx=0;
    endif;
    u2 = u.^2;
    hhat = recserar(b[_numx+1]+_shape(u2,garchq,1)*b[_numx+garchp+2:_numx+garchp+garchq+1],
                (meanc(u2)-meanc(u).^2)*ones(garchp,1),b[_numx+2:_numx+garchp+1]);
            @ initial value for hhat is unconditional error variance @

    if limits==1;
        __title = "Wald confidence limits";
        cl1 = cmlclimits(b,vcv);
        call cmlclprt(b,f,g,cl1,ret);
    else;
        cl1 = miss(0,0);
@       print "parameter estimates" b;          @
@       print "variance covariance matrix for parameter estimates" vcv;  @
    endif;
    retp(b,f,g,cl1,ret,hhat);
endp;


@ %%%%%%%%%% some supplementary procedures %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @

proc garchll( b, x );
  local u,u2,kappa,delta,alpha,h;
    if _numx>0;
        u = (x[.,1] - x[.,2:cols(x)] * b[1:_numx]);
    else;
        u = (x[.,1]);
    endif;
    u2 = u.^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)-meanc(u).^2)*ones(_p,1),delta);

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

proc grd(b,x);      @ original code for the gradient @
  local j,u,u2,kappa,delta,alpha,r,h,s,g;
  local p,q,numx;
  p = _p; q = _q; numx = _numx;

  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)-meanc(u).^2)*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;

@ Modified gradient code, trying to deal with the case where    @
@ there is NOTHING in the model for the mean. Useful for        @
@ estimating autoregressive conditional duration (ACD) models   @
proc grd2(b,x);
    local j,u,u2,kappa,delta,alpha,r,h,s,g;

    if _numx>0;
        u = x[.,1] - x[.,2:cols(x)] * b[1:_numx];
    else;
        u = x[.,1];
    endif;
    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)-meanc(u).^2)*ones(_p,1),delta);

    if _numx>0;
        g[.,1:_numx] = -2*_dshape(u.*x[.,2:cols(x)],_q,alpha);
    endif;
    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)));

    if _numx>0;
        g[.,1:_numx] = g[.,1:_numx] + (u./h).*x[.,2:cols(x)];
    endif;

    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;
