/**********    t-GARCH (Schoenberg 95)       *******************
Below is a CML program for a t-distributed
GARCH(p,q).  It introduces a new parameter, 'nu', which is the "degrees of
freedom".  In this situation, I presume, large values of 'nu' may be
interpreted as an absence of kurtosis, and small values as some degree of
leptokurtosis (Again CML comes through here because 'nu' must be
constrained to be greater than 2 since moments don't exist for 'nu' <= 2)

Below are comparative results for a GARCH(1,1) using the t-distribution
versus the normal distribution as applied to a time series of a stock
index (EAFE).  The estimate of 'nu' is 8.88 which suggests, perhaps, some
'fatness' of the tail.  However, it has a rather wide confidence interval.

I'm not sure how one would go about showing there's kurtosis.
The likelihood ratio chi-squared statistic, 267 * (3.05438 - 3.04723), is
decidely not significant and therefore we are not able to conclude that
the t-distributed model is better than the Normal model.  On the other
hand, the upper bound of the 'nu' parameter is than the 30 one could
choose as indistinguishable from the Normal.


TGARCH.PRG below doesn't have gradients, and, possibly as a result of
this, I had difficulty getting estimates for GARCH(3,3) and above.  I
will try to remedy the absence of a gradient in the near future.


************************************************************************
*    Ronald     *      Aptech Systems        *   rons@Aptech.com       *
*  Schoenberg   * University of Washington   *   rons@u.washington.edu *
************************************************************************




===============================================================================
                   Student's t distribution results
===============================================================================
 CML Version 1.0.14                                          9/07/95  10:47 am
===============================================================================

return code =    0
normal convergence

Mean log-likelihood        -3.04723
Number of cases     267


                             0.95 confidence limits
Parameters    Estimates     Lower Limit   Upper Limit   Gradient
------------------------------------------------------------------
CONST            1.1304        0.5311        1.7296      0.0000
kappa            4.3345        0.0010       12.9626     -0.0000
delta1           0.7599        0.3709        0.9182     -0.0000
alpha1           0.0808        0.0000        0.2017     -0.0000
nu               8.8867        2.0010       19.0241      0.0000

Number of iterations    13
Minutes to convergence     0.18767



===============================================================================
                     Normal distribution results
===============================================================================
 CML Version 1.0.14                                          9/07/95  10:47 am
===============================================================================

return code =    0
normal convergence

Mean log-likelihood        -3.05438
Number of cases     267


                             0.95 confidence limits
Parameters    Estimates     Lower Limit   Upper Limit   Gradient
------------------------------------------------------------------
CONST            1.0872        0.4729        1.7015      0.0000
kappa            4.8618        0.0100       13.4253     -0.0000
delta1           0.7436        0.3668        0.9155     -0.0000
alpha1           0.0745        0.0000        0.1776     -0.0000

Number of iterations    22
Minutes to convergence     0.17667


*****************************************************************/




/*
**  TGARCH.PRG  -  August 15, 1995
**
**  The following GAUSS source code solves the GARCH(p,q) problem via
**  constrained maximum likelihood estimation based on Student's
**  t-distribution.
**
**  It is in the public domain, and the author accepts no responsibility
**  for its actual use.  It requires GAUSS and the Constrained Maximum
**  Likelihood module.
**
************************************************************************
*    Ronald     *      Aptech Systems        *   rons@Aptech.com       *
*  Schoenberg   * University of Washington   *   rons@u.washington.edu *
************************************************************************
*/


library cml;
#include cml.ext;
cmlset;

dsn = "indx";

depvar = { EAFE };
indvar = { };

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


output file = tgarch.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 at each call to the function.
*/

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)) |
         "nu";
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)) |
         "nu";
endif;


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

_ww_ = { -1e250 1e250 };
_cml_Bounds = ones(numx+p+q+2,2).*_ww_;
_cml_Bounds[numx+1,1] = .001;                      /* kappa  >= .001   */
_cml_Bounds[numx+2:numx+p+q+1,1] = zeros(p+q,1);   /* alpha,delta >= 0 */
_cml_Bounds[numx+p+q+2,1] = 2.001;                 /* nu >= 2.001      */

/*
**  GARCH type constraint - parameters constrained to sum to less than one
**  comment out if IGARCH wanted
*/
_cml_c = zeros(1,rows(_cml_Parnames));
_cml_c[1,numx+2:numx+p+q+1] = -ones(1,p+q);
_cml_d = -.999;                               /* sum(alpha,delta) < .999 */

/*
**  IGARCH type constraint - parameters sum to exactly one
**  uncomment below if IGARCH wanted
**
** _cml_a = zeros(1,rows(_cml_Parnames));
** _cml_a[1,numx+2:numx+p+q+1] = ones(1,p+q);
** _cml_b = 1;                                /* sum(alpha,delta) = 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
*/

start = z[.,1] / z[.,2:cols(z)] | 1 | .25*ones(p+q,1) | 4;

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

__title = "ordinary Wald confidence limits";
call cmlclprt(b,f,g,cmltlimits(b,h),ret);
print;
print;
__title = "inequality constrained Wald confidence limits";
call cmlclprt(b,f,g,cmlclimits(b,h),ret);

output off;


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

proc garch( b, x );
  local u2,kappa,delta,alpha,h,nu,C;
if ismiss(b); retp(error(0)); endif;
  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];
  nu = b[numx+p+q+2];

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

  C = lnfact(0.5*(nu+1)-1) - 0.5*ln(pi*(nu-2)) - lnfact(0.5*nu-1);

  retp( C - 0.5*ln(h) - ((nu+1)/2)*ln(1 + u2 ./ ((nu -2)*h)) );
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;
