/* gretl - The Gnu Regression, Econometrics and Time-series Library
* Copyright (C) 1999-2006 Allin Cottrell and Riccardo "Jack" Lucchetti
*
* 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.
*/
/*
Modular GARCH routines by Jack Lucchetti, October 2006. For the
moment, meant to replace seamlessly fcp.c, but syntax should evolve in
the future.
*/
#include "libgretl.h"
#include "libset.h" /* unused */
#include "mod_garch.h"
#define GDEBUG 0
enum {
INIT_VAR_THEO,
INIT_VAR_OLS,
INIT_VAR_RESID
};
enum {
DIST_NORM,
DIST_T
};
typedef struct garch_container_ garch_container;
struct garch_container_ {
double *y; /* dependent variable */
const double **X; /* regressors (constant excluded) */
int t1; /* beginning of sample */
int t2; /* end of sample */
int nobs; /* number of observations */
int ncm; /* number of regressors (constant excluded) */
int p; /* GARCH p */
int q; /* GARCH q */
int k; /* total number of parameters */
int init; /* h0 initialisation method */
int distrib; /* innovations distribution (only Gaussian for now) */
double *e; /* residuals */
double *e2; /* squared residuals */
double *h; /* conditional variance */
int ascore; /* 1 for analytical score provided */
double **score_e; /* derivatives of the residuals wrt the parameters */
double **score_h; /* derivatives of the variances wrt the parameters */
double **blockglue; /* derivatives of the loglik wrt residuals and variances */
double **G; /* score matrix */
double *tot_score; /* score vector (sum of G) */
};
#if 0
static void mark (int *n)
{
if (n == NULL) {
fprintf(stderr,"Ha!\n");
} else {
fprintf(stderr,"Ha! (%d)\n", *n);
*n += 1;
}
}
#endif
#define MOD_DEBUG 0
static void free_eh_derivs (garch_container *DH)
{
doubles_array_free(DH->score_e, DH->k);
doubles_array_free(DH->score_h, DH->k);
doubles_array_free(DH->G, DH->k);
doubles_array_free(DH->blockglue, 2);
}
static int allocate_eh_derivs (garch_container *DH)
{
int k = DH->k;
int n = DH->nobs;
int err = 0;
DH->score_e = doubles_array_new(k, n);
DH->score_h = doubles_array_new(k, n);
DH->G = doubles_array_new(k, n);
DH->blockglue = doubles_array_new(2, n);
if (DH->score_e == NULL ||
DH->score_h == NULL ||
DH->G == NULL ||
DH->blockglue == NULL) {
free_eh_derivs(DH);
err = E_ALLOC;
}
return err;
}
static garch_container *
garch_container_new (double *y, const double **X,
int t1, int t2, int nobs, int nx,
int p, int q, int init_method,
double *res, double *h,
int analytical)
{
garch_container *DH = malloc(sizeof *DH);
if (DH == NULL) {
return NULL;
}
DH->e2 = malloc(nobs * sizeof *DH->e2);
if (DH->e2 == NULL) {
free(DH);
return NULL;
}
DH->y = y;
DH->X = X;
DH->t1 = t1;
DH->t2 = t2;
DH->nobs = nobs;
DH->ncm = nx;
DH->p = p;
DH->q = q;
DH->init = init_method;
DH->e = res;
DH->h = h;
DH->ascore = analytical;
DH->k = 1 + nx + 1 + p + q;
DH->score_e = NULL;
DH->score_h = NULL;
DH->G = NULL;
DH->blockglue = NULL;
if (DH->ascore) {
if (allocate_eh_derivs(DH)) {
free(DH->e2);
free(DH);
DH = NULL;
}
}
return DH;
}
static void garch_container_destroy (garch_container *DH)
{
if (DH->ascore) {
free_eh_derivs(DH);
}
free(DH->e2);
free(DH);
}
/* *ARCH log-likelihood for Gaussian innovations */
static double normal_ll (const garch_container *DH)
{
double e2t, ht, ll = 0.0;
int t;
for (t=DH->t1; t<=DH->t2; t++) {
e2t = DH->e2[t];
ht = DH->h[t];
if (na(e2t) || na(ht)) {
return NADBL;
}
ll -= log(ht) + e2t / ht;
}
ll *= 0.5;
ll -= (DH->t2 - DH->t1 + 1) * LN_SQRT_2_PI;
return ll;
}
static void normal_score (const garch_container *DH)
{
double ut;
int t;
for (t=DH->t1; t<=DH->t2; t++) {
DH->blockglue[0][t] = ut = -DH->e[t] / DH->h[t];
DH->blockglue[1][t] = 0.5 * (ut * ut - 1.0 / DH->h[t]);
}
}
/* Compute the GARCH quantities */
static int garch_etht (const double *par, void *ptr)
{
garch_container *DH = (garch_container *) ptr;
int t1 = DH->t1;
int t2 = DH->t2;
int p = DH->p;
int q = DH->q;
int maxlag = (p > q)? p : q;
int i, j, k, ret = 0;
int ncm = DH->ncm;
double **dedq = DH->score_e;
double **dhdq = DH->score_h;
int t, T = t2 - t1 + 1;
double et, ht, tmp, h0 = 0.0;
double u_var = 0.0;
/* compute residuals */
tmp = 0.0;
for (t = t1-maxlag; t <= t2; t++) {
if (t < t1) {
et = 0.0;
} else {
et = DH->y[t] - par[0];
if (DH->X != NULL) {
for (i=0; i<ncm; i++) {
et -= DH->X[i][t]*par[i+1];
}
}
DH->e[t] = et;
DH->e2[t] = et * et;
tmp += DH->e2[t];
}
}
if (DH->ascore) {
for (t=t-maxlag; t<t1; t++) {
for (i=0; i<DH->k; i++) {
dedq[i][t] = 0.0;
}
}
}
/* h0 and derivatives */
switch (DH->init) {
case INIT_VAR_OLS:
h0 = 1.0;
break;
case INIT_VAR_RESID:
h0 = tmp / T;
break;
case INIT_VAR_THEO:
tmp = 1.0;
for (i=ncm+2; i<DH->k; i++) {
tmp -= par[i];
}
u_var = par[ncm+1] / tmp;
h0 = u_var;
break;
}
for (t=t1-maxlag; t<t1; t++) {
DH->h[t] = h0;
DH->e2[t] = h0;
}
if (DH->ascore) {
double dh0;
switch (DH->init) {
case INIT_VAR_OLS:
for (t=t1-maxlag; t<t1; t++) {
for (i=0; i<DH->k; i++) {
dhdq[i][t] = 0.0;
}
}
break;
case INIT_VAR_RESID:
dh0 = 0.0;
for (t=t1; t<=t2; t++) {
dh0 -= DH->e[t];
}
for (t=t1-maxlag; t<t1; t++) {
dhdq[0][t] = dh0 * 2.0 / T;
}
for (i=0; i<ncm; i++) {
dh0 = 0.0;
for (t=t1; t<=t2; t++) {
dh0 -= DH->e[t] * DH->X[i][t];
}
for (t=t1-maxlag; t<t1; t++) {
dhdq[i+1][t] = dh0 * 2.0 / T;
}
}
for (t=t1-maxlag; t<t1; t++) {
for (i=ncm+1; i<DH->k; i++) {
dhdq[i][t] = 0.0;
}
}
break;
case INIT_VAR_THEO:
for (t=t1-maxlag; t<t1; t++) {
for (i=0; i<=ncm; i++) {
dhdq[i][t] = 0.0;
}
}
dh0 = u_var / par[ncm+1];
for (t=t1-maxlag; t<t1; t++) {
dhdq[ncm+1][t] = dh0;
}
dh0 *= u_var;
for (t=t1-maxlag; t<t1; t++) {
for (i=ncm+2; i<DH->k; i++) {
dhdq[i][t] = dh0;
}
}
break;
}
}
/* in-sample loop */
for (t=t1; t<=t2; t++) {
ht = par[ncm+1];
for (i=1; i<=p; i++) {
ht += DH->e2[t-i] * par[ncm+i+1];
}
for (i=1; i<=q; i++) {
ht += DH->h[t-i] * par[ncm+i+p+1];
}
DH->h[t] = ht;
if (DH->ascore) {
/* constant */
dedq[0][t] = -1.0;
k = ncm+1;
dhdq[0][t] = 0.0;
for (i=1; i<=p; i++) {
if (t - p < t1 && DH->init == INIT_VAR_RESID) {
dhdq[0][t] += par[k+i] * dhdq[0][t1-1];
} else {
dhdq[0][t] += 2.0 * par[k+i] * DH->e[t-i] * dedq[0][t-i];
}
}
/* regressors */
for (i=1; i<=ncm; i++) {
dedq[i][t] = -(DH->X[i-1][t]);
k = ncm+1;
dhdq[i][t] = 0.0;
for (j=1; j<=p; j++) {
if (t - p < t1 && DH->init == INIT_VAR_RESID) {
// add INIT_THEO here
dhdq[i][t] += par[k+j] * dhdq[i][t1-1];
} else {
dhdq[i][t] += 2.0 * par[k+j] * DH->e[t-j] * dedq[i][t-j];
}
}
}
/* garch params: omega */
dedq[ncm+1][t] = 0.0;
dhdq[ncm+1][t] = 1.0;
if (t - p < t1 && DH->init == INIT_VAR_THEO) {
for (i=1; i<=p; i++) {
dhdq[ncm+1][t] += par[ncm+1+i] * dhdq[ncm+1][t1-1];
}
}
/* garch params: alphas */
k = ncm + 2;
for (i=1; i<=p; i++) {
dedq[k][t] = 0.0;
dhdq[k][t] = DH->e2[t-i];
if (t - p < t1 && DH->init == INIT_VAR_THEO) {
for (j=0; j<p; j++) {
dhdq[k][t] += par[k+j] * dhdq[k][t1-1];
}
}
k++;
}
/* garch params: betas */
k = ncm + p + 2;
for (i=1; i<=q; i++) {
dedq[k][t] = 0.0;
dhdq[k][t] = DH->h[t-i];
if (t - p < t1 && DH->init == INIT_VAR_THEO) {
for (j=0; j<p; j++) {
dhdq[k][t] += par[k+j-p] * dhdq[k][t1-1];
}
}
k++;
}
/* "real" recursive part */
for (i=0; i<DH->k; i++) {
k = ncm + p + 2;
for (j=1; j<=q; j++) {
dhdq[i][t] += par[k++] * dhdq[i][t-j];
}
}
}
}
#if MOD_DEBUG
fputs("\n\n", stderr);
for (i=0; i<DH->k; i++) {
fprintf(stderr, "garch_etht: par[%d] = %9.6f ", i, par[i]);
}
fputc('\n', stderr);
for (t=t1-maxlag; t<=20; t++) {
if (t < t1) {
fputc('*', stderr);
} else {
fputc(' ', stderr);
}
fprintf(stderr, " t:%4d ", t);
fprintf(stderr, " %8.4f", DH->e[t]);
fprintf(stderr, " %8.4f", DH->e2[t]);
fprintf(stderr, " %8.4f", DH->h[t]);
fprintf(stderr, " %12.8f", dedq[ncm+2][t]);
fprintf(stderr, " %12.8f", dhdq[ncm+2][t]);
fputc('\n', stderr);
}
#endif
return ret;
}
static double loglik (const double *theta, void *ptr)
{
garch_container *DH = (garch_container *) ptr;
double ll = NADBL;
int err;
err = garch_etht(theta, DH);
if (!err) {
ll = normal_ll(DH);
}
return ll;
}
static int score_fill_matrices (const double *theta, void *ptr)
{
garch_container *DH = (garch_container *) ptr;
int i, t, err;
err = garch_etht(theta, DH);
if (err) {
return err;
}
normal_score(DH);
for (t=DH->t1; t<=DH->t2; t++) {
for (i=0; i<DH->k; i++) {
DH->G[i][t] = DH->score_e[i][t] * DH->blockglue[0][t] +
DH->score_h[i][t] * DH->blockglue[1][t];
}
}
return err;
}
static int anal_score (double *theta, double *s, int npar, BFGS_CRIT_FUNC ll,
void *ptr)
{
garch_container *DH = (garch_container *) ptr;
double tmp;
int t, i, err;
err = score_fill_matrices(theta, DH);
if (err) {
return err;
}
for (i=0; i<npar; i++) {
tmp = 0.0;
for (t=DH->t1;t<=DH->t2; t++) {
tmp += DH->G[i][t];
}
s[i] = tmp;
}
return err;
}
static int garch_iinfo (garch_container *DH, gretl_matrix *info)
{
double **tmp_info;
double tmpi, tmpj, tmpx1, tmpx2, x;
int i, j, t;
if (info == NULL) {
return E_ALLOC;
}
tmp_info = doubles_array_new(DH->k, DH->k);
if (tmp_info == NULL) {
return E_ALLOC;
}
for (i=0; i<DH->k; i++) {
for (j=0; j<=i; j++) {
tmp_info[i][j] = 0.0;
}
}
for (t=DH->t1; t<=DH->t2; t++) {
for (i=0; i<=DH->ncm; i++) {
tmpi = DH->score_h[i][t] / DH->h[t];
tmpx1 = (i==0)? 2.0 : 2.0*(DH->X[i-1][t]);
tmpx1 /= DH->h[t];
for (j=0; j<=i; j++) {
tmpj = DH->score_h[j][t] / DH->h[t];
tmpx2 = (j==0)? 1.0 : 1.0*(DH->X[j-1][t]);
x = tmpx1 * tmpx2 + tmpi * tmpj;
tmp_info[i][j] += x;
}
}
for (i=DH->ncm+1; i<DH->k; i++) {
tmpi = DH->score_h[i][t] / DH->h[t];
for (j=DH->ncm+1; j<=i; j++) {
tmpj = DH->score_h[j][t] / DH->h[t];
x = tmpi * tmpj;
tmp_info[i][j] += x;
}
}
}
for (i=0; i<DH->k; i++) {
for (j=0; j<=i; j++) {
gretl_matrix_set(info, i, j, 0.5*tmp_info[i][j]);
if (j < i) {
gretl_matrix_set(info, j, i, 0.5*tmp_info[i][j]);
}
}
}
doubles_array_free(tmp_info, DH->k);
#if GDEBUG
gretl_matrix_print(info, "Information matrix");
#endif
gretl_invert_symmetric_matrix(info);
#if GDEBUG
gretl_matrix_print(info, "Information matrix (inverse)");
#endif
return 0;
}
static int garch_opg (garch_container *DH, gretl_matrix *GG)
{
double **tmp_GG;
double tmpi, x;
int t, i, j;
if (GG == NULL) {
return E_ALLOC;
}
tmp_GG = doubles_array_new(DH->k, DH->k);
if (tmp_GG == NULL) {
return E_ALLOC;
}
for (i=0; i<DH->k; i++) {
for (j=0; j<DH->k; j++) {
tmp_GG[i][j] = 0.0;
}
}
for (t=DH->t1; t<=DH->t2; t++) {
for (i=0; i<DH->k; i++) {
tmpi = DH->G[i][t];
for (j=0; j<=i; j++) {
x = tmpi * DH->G[j][t];
tmp_GG[i][j] += x;
}
}
}
for (i=0; i<DH->k; i++) {
for (j=0; j<=i; j++) {
gretl_matrix_set(GG, i, j, tmp_GG[i][j]);
if (j < i) {
gretl_matrix_set(GG, j, i, tmp_GG[i][j]);
}
}
}
doubles_array_free(tmp_GG, DH->k);
#if GDEBUG
gretl_matrix_print(GG, "OPG matrix");
#endif
return 0;
}
static int garch_ihess (garch_container *DH, double *theta, gretl_matrix *invH)
{
double vij, *V;
int i, j, k, npar = DH->k;
int err = 0;
if (invH == NULL) {
return E_ALLOC;
}
V = numerical_hessian(theta, npar, loglik, DH, &err);
if (V == NULL) {
return err;
}
k = 0;
for (i=0; i<npar; i++) {
for (j=i; j<npar; j++) {
vij = V[k++];
gretl_matrix_set(invH, i, j, vij);
if (i != j) {
gretl_matrix_set(invH, j, i, vij);
}
}
}
#if GDEBUG
gretl_matrix_print(invH, "Hessian (inverse)");
#endif
free(V);
return 0;
}
static int
garch_covariance_matrix (int vopt, double *theta, garch_container *DH,
double *vcv)
{
gretl_matrix *V = NULL;
gretl_matrix *GG = NULL;
gretl_matrix *iinfo = NULL;
gretl_matrix *invhess = NULL;
int npar = DH->k;
int err = 0;
if (vopt == VCV_BW || vopt == VCV_QML) {
/* need a distinct matrix for workspace */
V = gretl_matrix_alloc(npar, npar);
if (V == NULL) {
return E_ALLOC;
}
}
if (vopt == VCV_OP || vopt == VCV_QML || vopt == VCV_BW) {
/* GG' needed */
GG = gretl_matrix_alloc(npar, npar);
err = garch_opg(DH, GG);
}
if (vopt == VCV_IM || vopt == VCV_BW) {
/* information matrix needed */
iinfo = gretl_matrix_alloc(npar, npar);
err = garch_iinfo(DH, iinfo);
}
if (vopt == VCV_QML || vopt == VCV_HESSIAN) {
/* Hessian matrix needed */
invhess = gretl_matrix_alloc(npar, npar);
err = garch_ihess(DH, theta, invhess);
}
if (err) {
goto bailout;
}
switch (vopt) {
case VCV_HESSIAN:
V = invhess;
invhess = NULL;
break;
case VCV_IM:
V = iinfo;
iinfo = NULL;
break;
case VCV_OP:
V = GG;
GG = NULL;
err = gretl_invert_symmetric_matrix(V);
break;
case VCV_BW:
gretl_matrix_qform(iinfo, GRETL_MOD_NONE, GG,
V, GRETL_MOD_NONE);
break;
case VCV_QML:
gretl_matrix_qform(invhess, GRETL_MOD_NONE, GG,
V, GRETL_MOD_NONE);
break;
default:
break;
}
#if GDEBUG
gretl_matrix_print(V, "Variance-covariance matrix");
#endif
if (!err) {
double vij;
int i, j;
for (i=0; i<npar; i++) {
for (j=i; j<npar; j++) {
vij = gretl_matrix_get(V, i, j);
vcv[(i*npar) + j] = vij;
if (i != j) {
vcv[(j*npar) + i] = vij;
}
}
}
}
bailout:
gretl_matrix_free(GG);
gretl_matrix_free(iinfo);
gretl_matrix_free(invhess);
gretl_matrix_free(V);
return err;
}
#if GDEBUG
static void test_score (garch_container *DH, double *theta)
{
double *testa, *testn;
int i, testret;
testa = malloc(DH->k * sizeof *testa);
testn = malloc(DH->k * sizeof *testn);
testret = anal_score(theta, testa, DH->k, loglik, DH);
testret = BFGS_numeric_gradient(theta, testn, DH->k, loglik, DH);
fprintf(stderr, "ret = %d:\n", testret);
for (i=0; i<DH->k; i++) {
fprintf(stderr, "g[%d]: analytical = %14.8f, numerical: = %14.8f, \n",
i, testa[i], testn[i]);
}
fputc('\n', stderr);
free(testa);
free(testn);
}
#endif
/*
Parameters to garch_estimate_mod()
t1: beginning of sample in auxiliary database
t2: end of sample in auxiliary database
nobs: total number of observations in auxiliary database
X: data matrix for auxiliary database (regressors, not needed on
output)
nx: number of columns of X
coeff: vector of coefficient for the conditional mean, normally
initialised by OLS on input (not needed on output) -- does NOT
include the constant
nc: number of elements in coeff
vcv: n^2 vector (0 on input) to store covariance matrix of coeff
res2: vector of 0's on input, squared resids on output (not needed)
e: vector of 0's on input, resids on output
h: null pointer on input, conditional variances on output
y: on input, vector with dep. var., not needed on output
amax: vector; element 0 holds the garch intercept; 1 and 2 the
arch & garch orders; from 3 onwards, the arch & garch
parameters
b: 0 on input, holds vector of coefficient for the conditional
mean on output
scale: double used to scale dep. var.
fncount: 0 on input, holds number of function evaluations on output
fncount: 0 on input, holds number of gradient evaluations on output
prn: print handle for info on iterations and other diagnostic output
*/
int garch_estimate_mod (int t1, int t2, int nobs,
const double **X, int nx, double *coeff, int nc,
double *vcv, double *res2, double *res, double *h,
double *y, double *amax, double *b,
double scale, int *fncount, int *grcount,
PRN *prn, int vopt)
{
garch_container *DH;
int p = amax[1];
int q = amax[2];
int npar = 1 + nx + 1 + p + q;
int analytical = 1;
double *theta = NULL;
int i, j, err = 0;
/* BFGS apparatus */
int maxit = 10000;
double reltol = 1.0e-13;
DH = garch_container_new(y, X, t1, t2, nobs, nx, p, q,
INIT_VAR_RESID, res, h,
analytical);
if (DH == NULL) {
return E_ALLOC;
}
theta = malloc(npar * sizeof *theta);
if (theta == NULL) {
garch_container_destroy(DH);
return E_ALLOC;
}
theta[0] = gretl_mean(t1, t2, y);
for (i=1; i<=nx; i++) {
theta[i] = coeff[i];
}
theta[nx+1] = amax[0];
for (i=nx+2, j=3; i<npar; i++) {
theta[i] = amax[j++];
}
#if 0
if (DH->ascore) {
fputs("\nUsing analytical score\n", stderr);
} else {
fputs("\nUsing numerical score\n", stderr);
}
#endif
err = BFGS_max(theta, npar, maxit, reltol,
fncount, grcount, loglik, C_LOGLIK,
(DH->ascore)? anal_score : NULL,
DH, (prn != NULL)? OPT_V : OPT_NONE,
prn);
amax[0] = loglik(theta, DH) - (t2 - t1 + 1) * log(scale);
#if GDEBUG
test_score(DH, theta);
#endif
err = garch_covariance_matrix(vopt, theta, DH, vcv);
if (!err) {
double sderr;
/* transcribe coefficients and standard errors
note: slightly different from fcp */
for (i=0, j=0; i<npar; i++) {
if (vcv[j] > 0.0) {
sderr = sqrt(vcv[j]);
} else {
sderr = 0.0;
}
j += npar + 1;
amax[i+1] = theta[i];
amax[i+1+npar] = sderr;
}
}
free(theta);
garch_container_destroy(DH);
return err;
}
syntax highlighted by Code2HTML, v. 0.9.1