Logo Search packages:      
Sourcecode: gretl version File versions

garch.c

/* gretl - The Gnu Regression, Econometrics and Time-series Library
 * Copyright (C) 1999-2004 by Allin Cottrell
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License 
 * as published by the Free Software Foundation; either version 2
 * of the License, or (at your option) any later version.
 *
 * This software 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
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License 
 * along with this software; if not, write to the 
 * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 * Boston, MA 02111-1307, USA.
 */

/* GARCH plugin for gretl using the Fiorentini, Calzolari and 
   Panattoni mixed-gradient algorithm.
*/

#include "libgretl.h"
#include "libset.h"
#include "var.h"

#include "fcp.h"
#include "mod_garch.h"

#undef VPARM_DEBUG

#define VPARM_MAX 6            /* max number of variance parameters */
#define GARCH_PARAM_MAX 0.999

double vparm_init[VPARM_MAX];

static void add_garch_varnames (MODEL *pmod, const DATAINFO *pdinfo,
                        const int *list)
{
    int p = list[1];
    int q = list[2];
    int r = list[0] - 4;
    int i, j, np = 2 + p + q + r;

    free(pmod->list);
    pmod->list = gretl_list_copy(list);

    gretl_model_allocate_params(pmod, np);
    if (pmod->errcode) {
      return;
    }

    strcpy(pmod->params[0], pdinfo->varname[0]);

    j = 1;
    for (i=0; i<r; i++) {
      if (pmod->list[5+i] > 0) {
          strcpy(pmod->params[j++], pdinfo->varname[pmod->list[5+i]]);
      }
    }

    strcpy(pmod->params[j++], "alpha(0)");

    for (i=0; i<q; i++) {
      sprintf(pmod->params[j++], "alpha(%d)", i + 1);
    }

    for (i=0; i<p; i++) {
      sprintf(pmod->params[j++], "beta(%d)", i + 1);
    }
}

static int make_packed_vcv (MODEL *pmod, double *vcv, int np,
                      int nc, double scale)
{
    const int nterms = np * (np + 1) / 2;
    double sfi, sfj;
    int i, j, k;

    free(pmod->vcv);
    pmod->vcv = malloc(nterms * sizeof *pmod->vcv);
    if (pmod->vcv == NULL) {
      return 1;  
    }

    for (i=0; i<np; i++) {
      if (i < nc) {
          sfi = scale;
      } else if (i == nc) {
          sfi = scale * scale;
      } else {
          sfi = 1.0;
      }
      for (j=0; j<=i; j++) {
          if (j < nc) {
            sfj = scale;
          } else if (j == nc) {
            sfj = scale * scale;
          } else {
            sfj = 1.0;
          }
          k = ijton(i, j, np);
          pmod->vcv[k] = vcv[i + np * j] * sfi * sfj;
      }
    }

    return 0;
}

static int write_garch_stats (MODEL *pmod, const double **Z,
                        double scale, const DATAINFO *pdinfo,
                        const int *list, const double *theta, 
                        int nparam, int pad, const double *res,
                        const double *h)
{
    int err = 0;
    double *coeff, *sderr, *garch_h;
    double den;
    int ynum = list[4];
    int nvp = list[1] + list[2];
    int xvars = list[0] - 4;
    int i;

    coeff = realloc(pmod->coeff, nparam * sizeof *pmod->coeff);
    sderr = realloc(pmod->sderr, nparam * sizeof *pmod->sderr);

    if (coeff == NULL || sderr == NULL) return 1;

    for (i=0; i<nparam; i++) {
      coeff[i] = theta[i+1];
      sderr[i] = theta[i+nparam+1];
    }
    
    pmod->coeff = coeff;
    pmod->sderr = sderr;
   
    pmod->ncoeff = nparam;

    pmod->ess = 0.0;
    for (i=pmod->t1; i<=pmod->t2; i++) {
      pmod->uhat[i] = res[i + pad] * scale;
      pmod->ess += pmod->uhat[i] * pmod->uhat[i];
      pmod->yhat[i] =  Z[ynum][i] * scale - pmod->uhat[i];
    }

    /* set sigma to its unconditional or steady-state value */
    den = 1.0;
    for (i=0; i<nvp; i++) {
      den -= coeff[i+xvars+1];
    }
    pmod->sigma = sqrt(coeff[xvars] / den);

    pmod->adjrsq = NADBL; 
    pmod->fstt = NADBL;

    mle_criteria(pmod, 1);

    pmod->ci = GARCH;
    pmod->ifc = 1;
    
    add_garch_varnames(pmod, pdinfo, list);

    /* add predicted error variance to model */
    garch_h = malloc(pdinfo->n * sizeof *garch_h);
    if (garch_h != NULL) {
      for (i=0; i<pdinfo->n; i++) {
          if (i < pmod->t1 || i > pmod->t2) {
            garch_h[i] = NADBL;
          } else {
            garch_h[i] = h[i + pad] * scale * scale;
          }
      }
      gretl_model_set_data(pmod, "garch_h", garch_h, 
                       MODEL_DATA_DOUBLE_ARRAY,
                       pdinfo->n * sizeof *garch_h);
    }

    return err;
}

static int make_garch_dataset (const int *list, double **Z,
                         int bign, int pad, int nx,
                         double **py, double ***pX)
{
    double *y = NULL, **X = NULL;
    int vx, vy = list[4];
    int i, k, s, t;

    /* If pad > 0 we have to create a newly allocated, padded
       dataset.  Otherwise we can use a virtual dataset, made
       up of pointers into the original dataset, Z. 
    */

    if (pad > 0) {
      y = malloc(bign * sizeof *y);
      if (y == NULL) {
          return E_ALLOC;
      }
      *py = y;
    } 

    if (nx > 0) {
      if (pad) {
          X = doubles_array_new(nx, bign);
      } else {
          X = malloc(nx * sizeof *X);
      }
      if (X == NULL) {
          free(y);
          return E_ALLOC;
      }
    }

    if (pad > 0) {
      /* build padded dataset */
      for (t=0; t<bign; t++) {
          if (t < pad) {
            y[t] = 0.0;
            for (i=0; i<nx; i++) {
                X[i][t] = 0.0;
            }
          } else {
            s = t - pad;
            y[t] = Z[vy][s];
            k = 5;
            for (i=0; i<nx; i++) {
                vx = list[k++]; 
                if (vx == 0) {
                  vx = list[k++];
                }
                X[i][t] = Z[vx][s];
            }
          }
      }
    } else {
      /* build virtual dataset */
      *py = Z[vy];
      k = 5;
      for (i=0; i<nx; i++) {
          vx = list[k++]; 
          if (vx == 0) {
            vx = list[k++];
          }
          X[i] = Z[vx];
      }
    }

    *pX = X;

    return 0;
}

static int get_vopt (int robust)
{
    int vopt = get_garch_vcv_version();
    int ropt = get_garch_robust_vcv_version();

    /* The defaults: QML if "robust" option is in force,
       otherwise negative Hessian */
    if (vopt == VCV_UNSET) {
      if (robust) {
          if (ropt == VCV_UNSET) {
            vopt = VCV_QML;
          } else {
            vopt = ropt;
          }
      } else {
          vopt = VCV_HESSIAN;
      }
    }

    return vopt;
}

static void garch_print_init (const double *coeff, int k,
                        const double *a, int p, int q,
                        int manual, PRN *prn)
{
    int i;

    pputc(prn, '\n');

    if (manual) {
      pputs(prn, "Manual initialization of parameters");
    } else {
      pputs(prn, "Automatic initialization of parameters");
    }

    pputs(prn, "\n\n Regression coefficients:\n");

    for (i=0; i<k; i++) {
      pprintf(prn, "  theta[%d] = %g\n", i, coeff[i]);
    }

    pputs(prn, "\n Variance parameters:\n");

    pprintf(prn, "  alpha[0] = %g\n", a[0]);
    for (i=0; i<p; i++) {
      pprintf(prn, "  alpha[%d] = %g\n", i+1, a[i+3]);
    }
    for (i=0; i<q; i++) {
      pprintf(prn, "   beta[%d] = %g\n", i, a[i+3+p]);
    }

    pputc(prn, '\n');
}

/* pick up any manually set initial values (if these
   have been set via "set initvals") */

static int garch_manual_init (double *a, int p, int q, 
                        double *coeff, double *b,
                        int k, PRN *prn)
{
    const gretl_matrix *m = get_init_vals();
    int mlen = gretl_vector_get_length(m);
    int n = p + q + 1 + k;
    int i, j;

    if (mlen != n) {
      if (mlen > 0) {
          fprintf(stderr, "Number of initvals = %d, but we want %d "
                "values for GARCH\n", mlen, n);
      }
      return 0;
    }

    /* coefficients on regressors */
    for (i=0; i<k; i++) {
      coeff[i] = m->val[i];
      b[i] = 0.0;
    }


    /* variance parameters */
    a[0] = m->val[i++];
    for (j=0; j<p+q; j++) {
      a[j+3] = m->val[i++];
    }

    garch_print_init(coeff, k, a, p, q, 1, prn);

    return 1;
}

int do_fcp (const int *list, double **Z, double scale,
          const DATAINFO *pdinfo, MODEL *pmod,
          PRN *prn, gretlopt opt)
{
    int t1 = pmod->t1, t2 = pmod->t2;
    int ncoeff = pmod->ncoeff;
    int p = list[1];
    int q = list[2];
    double *y = NULL;
    double **X = NULL;
    double *h = NULL;
    double *amax = NULL; 
    double *res = NULL, *res2 = NULL;
    double *coeff = NULL, *b = NULL;
    double *vcv = NULL;
    int fnc = 0, grc = 0, iters = 0;
    int nobs, maxlag, bign, pad = 0;
    int i, nx, nparam, vopt;
    int err = 0;

    vopt = get_vopt(opt & OPT_R);

    nx = ncoeff - 1;
    maxlag = (p > q)? p : q; 
    nparam = ncoeff + p + q + 1;

    nobs = t2 + 1; /* number of obs in full dataset */

    if (maxlag > t1) {
      /* need to pad data series at start */
      pad = maxlag - t1;
    } 

    /* length of series to pass to garch_estimate */
    bign = nobs + pad;
      
    res2 = malloc(bign * sizeof *res2);
    res = malloc(bign * sizeof *res);
    h = malloc(bign * sizeof *h);
    amax = malloc(bign * sizeof *amax);
    if (res2 == NULL || res == NULL || 
      amax == NULL || h == NULL) {
      err = E_ALLOC;
      goto bailout;
    }

    for (i=0; i<bign; i++) {
      res2[i] = res[i] = amax[i] = 0.0;
    }   
 
    coeff = malloc(ncoeff * sizeof *coeff);
    b = malloc(ncoeff * sizeof *b);
    if (coeff == NULL || b == NULL) {
      err = E_ALLOC;
      goto bailout;     
    }

    vcv = malloc((nparam * nparam) * sizeof *vcv);
    if (vcv == NULL) {
      err = E_ALLOC;
      goto bailout;
    }

    for (i=0; i<nparam * nparam; i++) {
      vcv[i] = 0.0;
    } 

    /* create dataset for garch estimation */
    err = make_garch_dataset(list, Z, bign, pad, nx, &y, &X);
    if (err) {
      goto bailout;
    }

    /* a bit odd, but for compatibility with FCP... */
    amax[1] = q;
    amax[2] = p; 

    if (!garch_manual_init(amax, p, q, coeff, b, ncoeff, prn)) {

      /* initial coefficients from OLS */
      for (i=0; i<ncoeff; i++) {
          coeff[i] = pmod->coeff[i];
          b[i] = 0.0;
      }

      /* initialize variance parameters */
      amax[0] = vparm_init[0];
      for (i=0; i<p+q; i++) {
          amax[i+3] = vparm_init[i+1];
      }

      if (opt & OPT_V) {
          garch_print_init(coeff, ncoeff, amax, p, q, 0, prn);
      }
    }

#ifdef USE_FCP
    err = garch_estimate(t1 + pad, t2 + pad, bign, 
                   (const double **) X, nx, coeff, ncoeff, 
                   vcv, res2, res, h, y, amax, b, scale, &iters,
                   prn, vopt);
#else
    if (getenv("FCP_GARCH") != NULL) {
      err = garch_estimate(t1 + pad, t2 + pad, bign, 
                       (const double **) X, nx, coeff, ncoeff, 
                       vcv, res2, res, h, y, amax, b, scale, &iters,
                       prn, vopt);
    } else {
      err = garch_estimate_mod(t1 + pad, t2 + pad, bign, 
                         (const double **) X, nx, coeff, ncoeff, 
                         vcv, res2, res, h, y, amax, b, scale, &fnc,
                         &grc, prn, vopt);
    }
#endif

    if (err != 0) {
      pmod->errcode = err;
    } else {
      int nparam = ncoeff + p + q + 1;

      for (i=1; i<=nparam; i++) {
          if (i <= ncoeff) {
            amax[i] *= scale;
            amax[i + nparam] *= scale;
          } else if (i == ncoeff + 1) {
            amax[i] *= scale * scale;
            amax[i + nparam] *= scale * scale;
          }
          pprintf(prn, "theta[%d]: %#14.6g (%#.6g)\n", i-1, amax[i], 
                amax[i + nparam]);
      }
      pputc(prn, '\n');

      pmod->lnL = amax[0];
      write_garch_stats(pmod, (const double **) Z, scale, pdinfo, 
                    list, amax, nparam, pad, res, h);
      make_packed_vcv(pmod, vcv, nparam, ncoeff, scale);
      if (iters > 0) {
          gretl_model_set_int(pmod, "iters", iters);
      } else if (fnc > 0) {
          gretl_model_set_int(pmod, "fncount", fnc);
          gretl_model_set_int(pmod, "grcount", grc);
      }
      gretl_model_set_int(pmod, "ml_vcv", vopt);
    }

 bailout:

    free(res2);
    free(res);
    free(h);
    free(amax);    
    free(coeff);
    free(b);
    free(vcv); 

    if (pad > 0) {
      free(y);
      doubles_array_free(X, nx);
    } else {
      free(X);
    }

    return err;
}

static int 
add_uhat_squared (const MODEL *pmod, double scale,
              double ***pZ, DATAINFO *pdinfo)
{
    int t, v = pdinfo->v;

    if (dataset_add_series(1, pZ, pdinfo)) {
      return E_ALLOC;
    }

    for (t=0; t<pdinfo->n; t++) {
      double u = pmod->uhat[t];

      if (na(u)) {
          (*pZ)[v][t] = NADBL;
      } else {
          u /= scale;
          (*pZ)[v][t] = u * u;
      }
    }

    strcpy(pdinfo->varname[v], "uhat2");

    return 0;
}

/*
  p and q are the GARCH orders
  ao = max(q,p) is the ar order
  mo = q is the ma order

  it is assumed that armapar contains the arma parameters 
  in the following order:
  armapar[0] : intercept
  armapar[1..ao] : ar terms
  armapar[ao+1..ao+mo] : ma terms
*/

static void
garchpar_from_armapar (const double *armapar, int q, int p)
{
    double x, sum_ab = 0.0;
    int ao = (p > q)? p : q;
    int mo = q;
    int i;

#ifdef VPARM_DEBUG
    for (i=0; i<1+ao+mo; i++) {
      fprintf(stderr, "armapar[%d] = %#12.6g\n", i, armapar[i]);
    }
#endif

    for (i=1; i<=p; i++) {
      x = 0.0;
      if (i <= ao) {
          x += armapar[i];
      } 
      if (i<=mo) {
          x += armapar[p+i];
      } 
      vparm_init[i] = (x < 0.0) ? 0.01 : x;
      sum_ab += vparm_init[i];
    }

    for (i=1; i<=q; i++) {
      x = armapar[p+i];
      vparm_init[p+i] = (x > 0.0) ? 0 : -x;
      sum_ab += vparm_init[p+i];
    }

#ifdef VPARM_DEBUG
    fprintf(stderr, "sum_ab = %#12.6g\n", sum_ab);
#endif

    if (sum_ab > GARCH_PARAM_MAX) {
      for (i=1; i<=p+q; i++) {
          vparm_init[i] *= GARCH_PARAM_MAX / sum_ab;
      }
      sum_ab = GARCH_PARAM_MAX;
    }

    vparm_init[0] = armapar[0];
}

static int 
garch_init_by_arma (const MODEL *pmod, const int *garchlist, 
                double scale, double ***pZ, DATAINFO *pdinfo)
{
    MODEL amod;
    int q = garchlist[1], p = garchlist[2];
    int v = pdinfo->v;
    int *list = NULL;
    int err = 0;

    /* for now we'll try this only for GARCH up to (2,2) */
    if (q > 2 || p > 2) {
      return 0;
    }

    /* add OLS uhat squared to dataset */
    if (add_uhat_squared(pmod, scale, pZ, pdinfo)) {
      return E_ALLOC;
    }

    list = gretl_list_copy(garchlist);
    if (list == NULL) {
      err = E_ALLOC;
      goto bailout;
    }

    list[1] = (q > p)? q : p;
    list[2] = q;
    /* dep var is squared OLS residual */
    list[4] = v;

    amod = arma(list, (const double **) *pZ, pdinfo, OPT_C, NULL);
    if (amod.errcode) {
      err = amod.errcode;
      goto bailout;
    } else {
      int i;

      model_count_minus();
      garchpar_from_armapar(amod.coeff, q, p);

      for (i=0; i<q+p+1; i++) {
          fprintf(stderr, "from ARMA: vparm_init[%d] = %#12.6g\n", i, 
                vparm_init[i]);
      }
    }

 bailout:

    dataset_drop_last_variables(pdinfo->v - v, pZ, pdinfo);
    free(list);

    return err;
}

/* sanity/dimension check */

static int *get_garch_list (const int *list, const double **Z,
                      const DATAINFO *pdinfo, int *err)
{
    int *glist = NULL;
    int i, p = list[1], q = list[2];
    int add0 = 1;

    *err = 0;

    /* rule out pure AR in variance (FIXME, why?) */
    if (p > 0 && q == 0) {
      gretl_errmsg_set(_("Error in garch command"));
      *err = E_DATA;
      return NULL;
    }

    /* rule out excessive total GARCH terms */
    if (p + q > 5) {
      gretl_errmsg_set(_("Error in garch command"));
      *err = E_DATA;
      return NULL;
    }

    /* insert constant if not present */
    for (i=4; i<=list[0]; i++) {
      if (list[i] == 0 || true_const(list[i], Z, pdinfo)) {
          add0 = 0;
          break;
      }
    }

    glist = gretl_list_new(list[0] + add0);
    if (glist == NULL) {
      *err = E_ALLOC;
    } else {
      for (i=1; i<=list[0]; i++) {
          glist[i] = list[i];
      }
      if (add0) {
          glist[i] = 0;
      }
    }

    return glist;
}

/* make regression list for initial OLS */

static int *make_ols_list (const int *list)
{
    int *olist;
    int i;

    olist = malloc((list[0] - 2) * sizeof *olist);
    if (olist == NULL) {
      return NULL;
    }

    olist[0] = list[0] - 3;
    for (i=4; i<=list[0]; i++) {
      olist[i-3] = list[i];
    }

    return olist;
}

#define GARCH_AUTOCORR_TEST 1

#if GARCH_AUTOCORR_TEST

int garch_pretest (MODEL *pmod, double ***pZ, DATAINFO *pdinfo,
               double *LMF, double *pvF)
{
    int err;

    err = autocorr_test(pmod, pdinfo->pd, pZ, pdinfo,
                  OPT_S | OPT_Q, NULL);

    if (!err) {
      *LMF = get_last_test_statistic(NULL);
      *pvF = get_last_pvalue(NULL);
    } 

    return err;
}

static void autocorr_message (double LMF, double pvF, int order, PRN *prn)
{
    if (!na(LMF) && pvF < 0.05) {
      pputs(prn, "\nConvergence was not reached.  One possible reason "
            "for this is\nautocorrelation in the error term.\n");
      pprintf(prn, "After estimating the model by OLS, the following result "
            "was\nobtained for a test of autocorrelation of order %d:\n",
            order);
      pprintf(prn, "LMF = %g, with p-value %g\n", LMF, pvF);
    }
}

#endif /* GARCH_AUTOCORR_TEST */

#define GARCH_SCALE_SIGMA 1

/* the driver function for the plugin */

MODEL garch_model (const int *cmdlist, double ***pZ, DATAINFO *pdinfo,
               PRN *prn, gretlopt opt) 
{
    MODEL model;
    int *list = NULL;
    int *ols_list = NULL;
    double LMF = NADBL;
    double pvF = NADBL;
    double scale = 1.0;
    int t, err, init_err, yno = 0;

    gretl_model_init(&model);

    list = get_garch_list(cmdlist, (const double **) *pZ,
                    pdinfo, &err);
    if (err) {
      model.errcode = err;
      return model;
    }

    ols_list = make_ols_list(list);
    if (ols_list == NULL) {
      model.errcode = E_ALLOC;
      return model;
    }

    /* run initial OLS */
    model = lsq(ols_list, pZ, pdinfo, OLS, OPT_A | OPT_M);
    if (model.errcode) {
      goto bailout;
    }

#if 0
    printmodel(&model, pdinfo, OPT_NONE, prn);
#endif

#if GARCH_AUTOCORR_TEST
    /* pretest the residuals for autocorrelation */
    if (prn != NULL) {
      garch_pretest(&model, pZ, pdinfo, &LMF, &pvF);
    }
#endif

#if GARCH_SCALE_SIGMA
    yno = ols_list[1];
    scale = model.sigma;
    for (t=0; t<pdinfo->n; t++) {
      if (!na((*pZ)[yno][t])) {
          (*pZ)[yno][t] /= scale;
      }
    }
    for (t=0; t<model.ncoeff; t++) {
      model.coeff[t] /= scale;
    }
    model.ess /= scale * scale;
    model.sigma = 1.0;
#endif 

    /* default variance parameter initialization */
    vparm_init[1] = 0.2;
    vparm_init[list[2]+1] = 0.7;
    vparm_init[0] = model.sigma * model.sigma * 0.1;

    if (opt & OPT_A) {
      /* "--arma-init": try initializing params via ARMA */
      init_err = garch_init_by_arma(&model, list, scale, 
                              pZ, pdinfo);
    }

    do_fcp(list, *pZ, scale, pdinfo, &model, prn, opt); 

    if (scale != 1.0) {
      /* undo scaling of dependent variable */
      for (t=0; t<pdinfo->n; t++) {
          if (!na((*pZ)[yno][t])) {
            (*pZ)[yno][t] *= scale;
          }
      }
    }

 bailout:

    free(ols_list);
    free(list);

#if GARCH_AUTOCORR_TEST
    if (!na(LMF)) {
      if (model.errcode == E_NOCONV) {
          autocorr_message(LMF, pvF, pdinfo->pd, prn);
      } else {
          gretl_model_destroy_tests(&model);
      }
    }
#endif

    return model;
}


Generated by  Doxygen 1.6.0   Back to index