Logo Search packages:      
Sourcecode: fseries version File versions

A2-2-GarchModelling.c

/* Copyright (C) 1997-1999  Adrian Trapletti
  
   This library is free software; you can redistribute it and/or
   modify it under the terms of the GNU Library General Public
   License as published by the Free Software Foundation; either
   version 2 of the License, or (at your option) any later version.
  
   This library is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   Library General Public License for more details.
  
   You should have received a copy of the GNU Library General Public
   License along with this library; if not, write to the Free
   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

   GARCH estimation 
   
   Reference: T. Bollerslev (1986): Generalized Autoregressive Conditional 
   Heteroscedasticity, Journal of Econometrics 31, 307-327. */


#include <R.h>


extern void F77_NAME(dsumsl) ();
extern void F77_NAME(dsmsno) ();
extern void F77_NAME(ddeflt) ();


#define BIG 1.0e+10  /* function value if the parameters are invalid */


static double dsqrarg;
#define DSQR(a) ((dsqrarg=(a)) == 0.0 ? 0.0 : dsqrarg*dsqrarg)

static double dmaxarg1,dmaxarg2;
#define DMAX(a,b) (dmaxarg1=(a),dmaxarg2=(b),(dmaxarg1) > (dmaxarg2) ?\
        (dmaxarg1) : (dmaxarg2))


struct garch_handler  /* used to set up the additional parameters used in calcf and calcg */
{
  double* y;  /* the time series to fit */
  double* h;  /* the conditional variance (cv) */
  double* dh;  /* dh_i/dp_j */
  int n;  /* the number of observations */
  int p, q;  /* GARCH(p,q) */
};

static struct garch_handler garch_h;


static void F77_SUB(calcf) (int *pq, double *p, int *nf, double *f, int *uiparm, 
                      double *urparm, void (*F77_SUB(ufparm))(void))
     /* compute negative log likelihood apart from the constant and the pre-sample values */
{
  int i, j, ok;
  int maxpq = DMAX(garch_h.p,garch_h.q); 
  double temp = 0.0;
  double sum = 0.0;
  
  ok = 1;
  if (p[0] <= 0.0) ok = 0;
  for (i=1; i<(*pq); i++)
    if (p[i] < 0.0) ok = 0;
  if (ok)  /* parameters are valid */
  {
    for (i=maxpq; i<garch_h.n; i++)  /* loop over time */
    {                                    /* compute cv at time i */
      temp = p[0];  /* compute ARCH part of cv */
      for (j=1; j<=garch_h.q; j++)
      temp += p[j]*DSQR(garch_h.y[i-j]);
      for (j=1; j<=garch_h.p; j++)  /* compute GARCH part of cv */
      temp += p[garch_h.q+j]*garch_h.h[i-j];      
      sum += log(temp)+DSQR(garch_h.y[i])/temp;  /* compute eq. 18 */
      garch_h.h[i] = temp;  /* assign cv at time i */
    }
    (*f) = 0.5*sum;
  }
  else  /* parameters are invalid */
    (*f) = BIG;
}

static void F77_SUB(calcg) (int *pq, double *p, int *nf, double *dp, int *uiparm, 
                      double *urparm, void (*F77_SUB(ufparm))(void))
     /* compute derivative of negative log likelihood */
{
  int i, j, k;
  int maxpq = DMAX(garch_h.p,garch_h.q); 
  double temp1, temp2, temp3;
  
  for (k=0; k<(*pq); k++)  /* initialize */
    dp[k] = 0.0;
  for (i=maxpq; i<garch_h.n; i++)  /* loop over time */
  {                                   /* compute cv at time i and derivatives dh_i/dp_j */
    temp1 = p[0];  /* compute ARCH part of cv */
    for (j=1; j<=garch_h.q; j++)
      temp1 += p[j]*DSQR(garch_h.y[i-j]);
    for (j=1; j<=garch_h.p; j++)  /* compute GARCH part of cv */
      temp1 += p[garch_h.q+j]*garch_h.h[i-j];   
    garch_h.h[i] = temp1;  /* assign cv at time i */
    temp2 = 0.5*(1.0-DSQR(garch_h.y[i])/temp1)/temp1;  /* compute dl_i/dh_i, eq. 19 */
    temp3 = 1.0;  /* compute dh_i/dp_0, eq. 21 */
    for (j=1; j<=garch_h.p; j++)
      temp3 += p[garch_h.q+j]*garch_h.dh[(*pq)*(i-j)];
    garch_h.dh[(*pq)*i] = temp3;  /* assign dh_i/dp_0 */
    dp[0] += temp2*temp3;  /* assign dl_i/dp_0 = dl_i/dh_i * dh_i/dp_0 */
    for (k=1; k<=garch_h.q; k++)  /* compute dl_i/dp_k for the ARCH part, eq. 19 */
    {
      temp3 = DSQR(garch_h.y[i-k]);  /* compute dh_i/dp_k, eq. 21 */
      for (j=1; j<=garch_h.p; j++)
      temp3 += p[garch_h.q+j]*garch_h.dh[(*pq)*(i-j)+k];
      garch_h.dh[(*pq)*i+k] = temp3;  /* assign dh_i/dp_k */
      dp[k] += temp2*temp3;  /* assign dl_i/dp_k = dl_i/dh_i * dh_i/dp_k */
    }
    for (k=1; k<=garch_h.p; k++)  /* compute dl_i/dp_k for the GARCH part, eq. 19 */
    {
      temp3 = garch_h.h[i-k];  /* compute dh_i/dp_k, eq. 21 */
      for (j=1; j<=garch_h.p; j++)
      temp3 += p[garch_h.q+j]*garch_h.dh[(*pq)*(i-j)+garch_h.q+k];
      garch_h.dh[(*pq)*i+garch_h.q+k] = temp3;  /* assign dh_i/dp_k */ 
      dp[garch_h.q+k] += temp2*temp3;  /* assign dl_i/dp_k = dl_i/dh_i * dh_i/dp_k */
    }
  } 
}

static void F77_SUB(ufparm) ()
{
  error ("fatal error in fit_garch ()\n");
}


void fit_garch (double *y, int *n, double *par, int *p, int *q, int *itmax, 
            double *eps, double *fret, int *agrad, int *trace)
     /* fit a GARCH (p, q) model 
      
      Input: 
      
      y[0..n-1]      series to fit
      par[0..p+q]    initial parameter estimates
      p, q           model orders
      itmax          maximum number of iterations
      eps            precision
      agrad          estimation with analytical/numerical gradient
      trace          output yes/no  
      
      Output:
      
      par[0..p+q]    parameter estimates at minimum
      fret           function value at minimum 
     */   
{
  int i, j, pq, liv, lv, alg;
  double *d, *v;
  int *iv;
  double var;
  
  /* set up general optimizer parameters to default values */ 
  pq = (*p)+(*q)+1; 
  d = Calloc (pq, double);
  for (i=0; i<pq; i++)
    d[i] = 1.0;
  liv = 60;
  iv = Calloc (liv, int);
  lv = 77+pq*(pq+17)/2;
  v = Calloc (lv, double);
  alg = 2;  
  F77_CALL(ddeflt) (&alg, iv, &liv, &lv, v);
  iv[0] = 12;  
  
  /* set up user defined optimizer parameters */ 
  iv[16] = 2*(*itmax);  
  iv[17] = (*itmax);  
  if (*trace)
    iv[20] = 6;  
  else
    iv[20] = 0;  
  v[30] = DMAX(1.0e-20, DSQR((*eps)));
  v[31] = DMAX(1.0e-10, pow((*eps),2.0/3.0));
  v[32] = sqrt((*eps));
  v[33] = 100.0*(*eps);
  
  /* set handler values */
  garch_h.p = (*p); garch_h.q = (*q); garch_h.n = (*n); garch_h.y = y;  
  garch_h.h = Calloc ((*n), double); 
  garch_h.dh = Calloc ((*n)*pq, double);
  var = 0.0;
  for (i=0; i<(*n); i++)  /* estimate unconditional variance (uv) */
    var += DSQR(y[i]);
  var /= (double) (*n); 
  for (i=0; i<DMAX((*p),(*q)); i++)  /* initialize */
  {
    garch_h.h[i] = var;  /* with uv */
    garch_h.dh[pq*i] = 1.0;  /* dh_i/dp_0 with 1 */
    for (j=1; j<pq; j++)
      garch_h.dh[pq*i+j] = 0.0;  /* dh_i/dp_j with 0 */
  }
  
  if (*agrad)  /* estimation with analytical gradient */
  {
    if (*trace) Rprintf ("\n ***** ESTIMATION WITH ANALYTICAL GRADIENT ***** \n\n");
    F77_CALL(dsumsl) (&pq, d, par, F77_SUB(calcf), F77_SUB(calcg), 
                  iv, &liv, &lv, v, NULL, NULL, F77_SUB(ufparm));
    if (*trace) Rprintf ("\n");
  }
  else  /* estimation with numerical gradient */
  {
    if (*trace) Rprintf ("\n ***** ESTIMATION WITH NUMERICAL GRADIENT ***** \n\n");
    F77_CALL(dsmsno) (&pq, d, par, F77_SUB(calcf), iv, 
                  &liv, &lv, v, NULL, NULL, F77_SUB(ufparm));
    if (*trace) Rprintf ("\n");
  }
  
  /* return function value */
  (*fret) = v[9];

  /* free memory */
  Free (d);
  Free (iv); 
  Free (v);
  Free (garch_h.h);
  Free (garch_h.dh);
}


void pred_garch (double *y, double *h, int *n, double *par, int *p, int *q, int *genuine)
     /* predict cv with a GARCH (p, q) model 
      
      Input: 
      
      y[0..n-1]    series to predict
      par[0..p+q]  parameters of the GARCH (p, q)
      p, q         model orders
      genuine      logical indicating if a genuine prediction is computed
      
      Output:
      
      h[0..N]      predicted cv, where N = n for genuine prediction, and N = n-1 otherwise
     */   
{
  double var, temp;
  int i, j, maxpq, N;
  
  if (*genuine) N = (*n)+1;
  else N = (*n);
  maxpq = DMAX((*p),(*q));
  var = 0.0;
  for (i=1; i<=(*p)+(*q); i++)  /* compute uv */
    var += par[i];
  var = par[0]/(1.0-var); 
  for (i=0; i<maxpq; i++)  /* initialize with uv */
    h[i] = var;
  for (i=maxpq; i<N; i++)  /* loop over time */
  {                                    /* compute cv at time i */
    temp = par[0];  /* compute ARCH part of cv */
    for (j=1; j<=(*q); j++)
      temp += par[j]*DSQR(y[i-j]);
    for (j=1; j<=(*p); j++)  /* compute GARCH part of cv */
      temp += par[(*q)+j]*h[i-j]; 
    h[i] = temp;  /* assign cv at time i */
  }
}


void ophess_garch (double *y, int *n, double *par, double *he, int *p, int *q)
     /* Compute outer product approximation of the hessian of the 
      negative log likelihood of a GARCH (p, q) model at given parameter
      estimates
      
      Input: 
      
      y[0..n-1]              time series
      par[0..p+q]            parameter estimates of the GARCH (p, q)
      p, q                   model orders
      
      Output:
      
      he[0..(p+q+1)*(p+q+1)-1]      predicted cv
     */   
{
  double var, temp1, temp2, temp3;
  int i, j, k, pq;
  double *h, *dh, *dpar;
  
  pq = (*p)+(*q)+1; 
  h = Calloc ((*n), double); 
  dh = Calloc ((*n)*pq, double);
  dpar = Calloc (pq, double);
  var = 0.0;
  for (i=0; i<(*n); i++)  /* estimate uv */
    var += DSQR(y[i]);
  var /= (double) (*n); 
  for (i=0; i<DMAX((*p),(*q)); i++)  /* initialize */
  {
    h[i] = var;  /* with uv */
    dh[pq*i] = 1.0;  /* dh_i/dp_0 with 1 */
    for (j=1; j<pq; j++)
      dh[pq*i+j] = 0.0;  /* dh_i/dp_j with 0 */
  }
  for (k=0; k<pq; k++)  /* initialize */
    for (j=0; j<pq; j++)
      he[pq*k+j] = 0.0;
  for (i=DMAX((*p),(*q)); i<(*n); i++)  /* loop over time */
  {                                   /* compute cv at time i and derivatives dh_i/dp_j */
    temp1 = par[0];  /* compute ARCH part of cv */
    for (j=1; j<=(*q); j++)
      temp1 += par[j]*DSQR(y[i-j]);
    for (j=1; j<=(*p); j++)  /* compute GARCH part of cv */
      temp1 += par[(*q)+j]*h[i-j];   
    h[i] = temp1;  /* assign cv at time i */
    temp2 = 0.5*(1.0-DSQR(y[i])/temp1)/temp1;  /* compute dl_i/dh_i, eq. 19 */
    temp3 = 1.0;  /* compute dh_i/dp_0, eq. 21 */
    for (j=1; j<=(*p); j++)
      temp3 += par[(*q)+j]*dh[pq*(i-j)];
    dh[pq*i] = temp3;  /* assign dh_i/dp_0 */
    dpar[0] = temp2*temp3;  /* assign dl_i/dp_0 = dl_i/dh_i * dh_i/dp_0 */
    for (k=1; k<=(*q); k++)  /* compute dl_i/dp_k for the ARCH part, eq. 19 */
    {
      temp3 = DSQR(y[i-k]);  /* compute dh_i/dp_k, eq. 21 */
      for (j=1; j<=(*p); j++)
      temp3 += par[(*q)+j]*dh[pq*(i-j)+k];
      dh[pq*i+k] = temp3;  /* assign dh_i/dp_k */
      dpar[k] = temp2*temp3;  /* assign dl_i/dp_k = dl_i/dh_i * dh_i/dp_k */
    }
    for (k=1; k<=(*p); k++)  /* compute dl_i/dp_k for the GARCH part, eq. 19 */
    {
      temp3 = h[i-k];  /* compute dh_i/dp_k, eq. 21 */
      for (j=1; j<=(*p); j++)
      temp3 += par[(*q)+j]*dh[pq*(i-j)+(*q)+k];
      dh[pq*i+(*q)+k] = temp3;  /* assign dh_i/dp_k */ 
      dpar[(*q)+k] = temp2*temp3;  /* assign dl_i/dp_k = dl_i/dh_i * dh_i/dp_k */
    }
    for (k=0; k<pq; k++)  /* compute outer product approximation, p. 317 */
      for (j=0; j<pq; j++)
      he[pq*k+j] += dpar[k]*dpar[j];
  }
  Free (h);
  Free (dh);
  Free (dpar);
}


Generated by  Doxygen 1.6.0   Back to index