/*
The functions here are derived from the Fortran in VSGARCMX.FOR,
from the Fiorentini, Calzolari and Panattoni GARCH implementation
(Journal of Applied Econometrics, 1996), the original notice for
which is reproduced infra. The relevant code was extracted from
its original Monte Carlo context, translated to C using f2c, and
then edited extensively to turn it into more idiomatic C, and
to convert from fixed-size arrays to dynamic memory allocation.
I have also modified the comments at certain places in the code,
where reference is made to the equations in the FCP paper. The
comments in the Fortran apparently pertained to a draft of the
paper; I have updated them relative to the paper as published in
JAE, 1996, pp. 399-417.
Allin Cottrell, Wake Forest University, March 2004.
*/
/* Gabriele FIORENTINI, Giorgio CALZOLARI, Lorenzo PANATTONI
Journal of APPLIED ECONOMETRICS, 1996
mixed gradient algorithm
garch(p,q) estimates of a linear equation
SEE BOLLERSLEV, JOE 31(1986), 307-327.
*/
#include "libgretl.h"
#include "libset.h"
#include "fcp.h"
#include "f2c.h"
#include "clapack_double.h"
#undef FDEBUG
#define ABNUM 3 /* max number of GARCH lags */
#define SMALL_HT 1.0e-7
int global_np;
double gscale;
int gncoeff;
#define vix(i,j) ((i) + global_np * (j))
/* private functions */
static int invert (double *g, int dim);
static
double garch_ll (double *c, int nc, double *res2,
double *res, double *yhat,
const double *y, const double **X,
int nexo, int t1, int t2, const double *param,
double *b, double *a0,
int q, int p, double *h);
static int vcv_setup (int t1, int t2, double *c, int nc,
double *res2, double *res, int *count,
const double **g, double *grad,
double *param, int nparam, double *a0,
int q, int p, double *h,
double **dhdp, double *zt,
double ***H, double *vco, int code);
static int
garch_info_matrix (int t1, int t2,
const double **X, int nexo,
double *yhat, double *c, int nc, double *res2,
double *res, const double *y, double toler,
int *count, double *vcv, const double **g,
double *grad, double *param, int nparam,
double *b, double *a0,
int q, int p, double *h,
double **dhdp, double *zt);
static int
garch_full_hessian (int t1, int t2,
const double **X, int nexo,
double *yhat, double *c, int nc, double *res2,
double *res, const double *y, double toler,
int *count, double *vcv, const double **g,
double *grad, double *param, int nparam,
double *b, double *a0,
int q, int p, double *h,
double **dhdp, double *zt);
static int vs_allocate (double ***pdhdp, double ***pg,
double **pparam, double **pgrad,
double **pc, double **paux, double **pvch,
double **pparpre, double **pyhat,
int np, int nrc, int T)
{
double *param = NULL, *grad = NULL;
double *c = NULL, *aux = NULL, *vch = NULL, *parpre = NULL;
double **D = NULL, **G = NULL;
double *yhat = NULL;
int t;
param = malloc(np * sizeof *param);
grad = malloc(np * sizeof *grad);
parpre = malloc(np * sizeof *parpre);
if (param == NULL || grad == NULL || parpre == NULL) {
goto bailout;
}
c = malloc(nrc * sizeof *c);
aux = malloc(nrc * sizeof *aux);
yhat = malloc(T * sizeof *yhat);
if (c == NULL || aux == NULL || yhat == NULL) {
goto bailout;
}
vch = malloc(np * np * sizeof *vch);
if (vch == NULL) {
goto bailout;
}
D = doubles_array_new(np, T);
if (D == NULL) {
goto bailout;
}
G = doubles_array_new(nrc, T);
if (G == NULL) {
goto bailout;
}
for (t=0; t<T; t++) {
yhat[t] = 0.0;
}
*pdhdp = D;
*pg = G;
*pparam = param;
*pgrad = grad;
*pc = c;
*paux = aux;
*pvch = vch;
*pparpre = parpre;
*pyhat = yhat;
return 0;
bailout:
free(param);
free(grad);
free(c);
free(aux);
free(vch);
free(parpre);
free(yhat);
doubles_array_free(D, np);
doubles_array_free(G, nrc);
return 1;
}
static void vs_free (double **dhdp, int np, double **g, int nrc,
double *param, double *grad,
double *c, double *aux, double *vcv,
double *parpre, double *yhat)
{
doubles_array_free(dhdp, np);
doubles_array_free(g, nrc);
free(param);
free(grad);
free(c);
free(aux);
free(vcv);
free(parpre);
free(yhat);
}
static void garch_iter_info (int iter, double *theta, int m, double ll,
int hess, PRN *prn)
{
int i;
double x;
pprintf(prn, "\n*** %s %d%s\n", ("iteration"), iter + 1,
(hess)? _(" (using Hessian)") : " (using Information Matrix)");
pputs(prn, "Parameters:\n");
for (i=0; i<m; i++) {
if (i && i % 5 == 0) {
pputc(prn, '\n');
}
x = theta[i];
if (i < gncoeff - 1) {
x *= gscale;
} else if (i == gncoeff - 1) {
x *= gscale * gscale;
}
pprintf(prn, "%#12.5g ", x);
}
pprintf(prn, "\nll = %f\n", ll);
}
static double
get_yhat (const double **X, int nx, int t, const double *a)
{
int j;
double yhat = a[0];
for (j = 0; j < nx; j++) {
yhat += a[j+1] * X[j][t];
}
return yhat;
}
static void copy_coeff (const double *c, int nc, double *b)
{
int i;
for (i=0; i<nc; i++) {
b[i] = c[i];
}
}
int ols_ (int t1, int t2, const double **X, int nx,
double *c, int nc, const double *y,
double *amax, double *aux, double *b, double **g)
{
int i, j, k, t, err;
double deltc, deriv, relinc = 0.5;
double oldc, yy;
double *vc;
vc = malloc(nc * nc * sizeof *vc);
if (vc == NULL) return 1;
copy_coeff(c, nc, b);
for (t=t1; t<=t2; t++) {
amax[t] = get_yhat(X, nx, t, b);
}
for (i=0; i<nc; i++) {
aux[i] = 0.0;
for (j=0; j<nc; j++) {
vc[i + j * nc] = 0.0;
}
}
for (t=t1; t<=t2; t++) {
for (k=0; k<nc; k++) {
oldc = c[k];
deltc = relinc;
if (oldc != 0.0) {
deltc = oldc * relinc;
}
c[k] = oldc + deltc;
copy_coeff(c, nc, b);
yy = get_yhat(X, nx, t, b);
deltc = c[k] - oldc;
deriv = (yy - amax[t]) / deltc;
c[k] = oldc;
g[k][t] = deriv;
}
copy_coeff(c, nc, b);
/* cumulates all the w'z into diagonal blocks of vc,
and w'y into elements of aux */
for (i=0; i<nc; i++) {
aux[i] += g[i][t] * y[t];
for (j=0; j<nc; j++) {
vc[i + j * nc] += g[i][t] * g[j][t];
}
}
}
err = invert(vc, nc);
if (err == 0) {
/* compute coefficients */
for (i=0; i<nc; i++) {
c[i] = 0.0;
for (j=0; j<nc; j++) {
c[i] += vc[i + j * nc] * aux[j];
}
}
copy_coeff(c, nc, b);
} else {
fputs("OLS: matrix is singular, initial coefficients are unchanged\n",
stderr);
for (i=0; i<nc; i++) {
for (j=0; j<nc; j++) {
vc[i + j * nc] = 0.0;
}
}
}
free(vc);
return 0;
}
/* Calculate robust VCV. If the first input (vci) is the inverse of
the Hessian, the output is the QML (White) estimator; if vci is the
inverse of the information matrix, the output is the
Bollerslev-Wooldridge estimator. The second input matrix (vco) is
the (uninverted) outer product (OP) matrix.
*/
static double *robust_vcv (const double *vci, const double *vco,
int nparam)
{
int i, j, k;
int np2 = nparam * nparam;
double *vv = NULL, *vcr = NULL;
vv = malloc(np2 * sizeof *vv); /* temporary workspace */
vcr = malloc(np2 * sizeof *vcr);
if (vv == NULL || vcr == NULL) {
free(vv);
free(vcr);
return NULL;
}
/* multiply H^{-1} (or I^{-1}) into OP */
for (i=0; i<nparam; i++) {
for (j=0; j<nparam; j++) {
vv[vix(i,j)] = 0.0;
for (k=0; k<nparam; k++) {
vv[vix(i,j)] += vci[vix(i,k)] * vco[vix(k,j)];
}
}
}
/* post-multiply by H^{-1} (or I^{-1}) */
for (i=0; i<nparam; i++) {
for (j=0; j<nparam; j++) {
vcr[vix(i,j)] = 0.0;
for (k=0; k<nparam; k++) {
vcr[vix(i,j)] += vv[vix(i,k)] * vci[vix(k,j)];
}
}
}
free(vv);
return vcr;
}
static int
make_garch_vcv (int t1, int t2,
const double **X, int nx,
double *yhat, const double *y,
double *c, int nc,
double *res2, double *res,
const double **g, double *grad,
double *param, int nparam,
double *b, double *a0,
int q, int p,
double *h, double **dhdp, double *zt,
const double *vch, double *vcv,
int vopt)
{
double *vco = NULL, *vcr = NULL, *vci = NULL;
int np2 = nparam * nparam;
int k, err = 0;
/* OP and robust variants need OP matrix */
if (vopt == VCV_OP || vopt == VCV_QML || vopt == VCV_BW) {
vco = malloc(nparam * nparam * sizeof *vco);
if (vco == NULL) {
err = 1;
goto bailout;
}
vcv_setup(t1, t2, c, nc,
res2, res,
NULL, g, grad,
param, nparam,
a0, q, p, h,
dhdp, zt, NULL, vco,
VCV_OP);
if (vopt == VCV_OP) {
err = invert(vco, nparam);
for (k=0; k<np2; k++) {
vcv[k] = vco[k];
}
}
}
/* IM and BW variants need the info matrix */
if (vopt == VCV_IM || vopt == VCV_BW) {
vci = malloc(nparam * nparam * sizeof *vci);
if (vci == NULL) {
err = 1;
goto bailout;
}
garch_info_matrix(t1, t2, X, nx, yhat, c,
nc, res2, res, y,
0.0, NULL, vci, g,
grad, param, nparam, b,
a0, q, p, h, dhdp, zt);
if (vopt == VCV_IM) {
for (k=0; k<np2; k++) {
vcv[k] = vci[k];
}
} else {
/* Bollerslev-Wooldridge */
vcr = robust_vcv(vci, vco, nparam);
if (vcr == NULL) {
err = 1;
goto bailout;
}
for (k=0; k<np2; k++) {
vcv[k] = vcr[k];
}
}
} else if (vopt == VCV_QML) {
vcr = robust_vcv(vch, vco, nparam);
if (vcr == NULL) {
err = 1;
goto bailout;
}
for (k=0; k<np2; k++) {
vcv[k] = vcr[k];
}
} else if (vopt == VCV_HESSIAN) {
for (k=0; k<np2; k++) {
vcv[k] = vch[k];
}
}
bailout:
free(vco);
free(vci);
free(vcr);
return err;
}
/*
Parameters to garch_estimate()
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)
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)
res: 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 dependent variable
iters: int, 0 on input, holds number of iterations on output
prn: print handle for info on iterations and other diagnostic output
*/
int garch_estimate (int t1, int t2, int nobs,
const double **X, int nx, double *coeff, int nc,
double *vcv, double *res2, double *res, double *h,
const double *y, double *amax, double *b,
double scale, int *iters, PRN *prn, int vopt)
{
int it1, it2, ittot;
int count, conv = 0;
int i, q, p, nparam;
double zt[6]; /* 6 == max value of (1 + q + p) */
double pappo, tol1, tol2;
double sumgra;
double ll, a0, s1, s2;
double *param = NULL, *grad = NULL;
double **dhdp = NULL, **g = NULL;
double *c = NULL, *aux = NULL;
double *parpre, *yhat = NULL;
double *vch = NULL;
int err = 0;
q = (int) amax[1];
p = (int) amax[2];
/* "export" scale as file-scope global */
gscale = scale;
/* number of parameters of unconcentrated likelihood */
nparam = nc + 1 + q + p;
global_np = nparam;
gncoeff = nc + 1;
if (vs_allocate(&dhdp, &g, ¶m, &grad,
&c, &aux, &vch, &parpre, &yhat,
nparam, nc, nobs)) {
pprintf(prn, "Out of memory\n");
return E_ALLOC;
}
for (i=0; i<q; i++) {
param[nc + i + 1] = amax[3 + i];
}
for (i=0; i<p; i++) {
param[nc + q + i + 1] = amax[3 + q + i];
}
param[nc] = amax[0];
tol1 = .05; /* tolerance before switching to Hessian */
tol2 = 1e-8;
for (i=0; i<nc; i++) {
c[i] = param[i] = coeff[i];
}
/* this is only to calculate matrix of regressors (g) */
ols_(t1, t2, X, nx, c, nc, y, amax, aux, b, g);
#ifdef FFDEBUG
for (i=0; i<nc; i++) {
fprintf(stderr, "after ols g[%d] = %.9g\n", i, g[i]);
}
#endif
/* iterative estimation */
count = 0;
for (it1=0; it1<100; it1++) {
ll = garch_ll(c, nc, res2, res, yhat, y,
X, nx, t1, t2, param, b, &a0,
q, p, h);
garch_iter_info(it1, param, nparam, ll, 0, prn);
/* store previous coefficients */
for (i=0; i<nparam; i++) {
parpre[i] = param[i];
}
#ifdef FDEBUG
fprintf(stderr, "*** Calling garch_info_matrix, round %d\n", it1);
#endif
err = garch_info_matrix(t1, t2, X, nx, yhat, c,
nc, res2, res, y,
tol1, &count, vch, (const double **) g,
grad, param, nparam, b,
&a0, q, p, h, dhdp, zt);
if (err) {
return E_NOCONV;
}
s1 = s2 = 0.0;
for (i=0; i<nparam; i++) {
s1 += parpre[i] * parpre[i];
pappo = param[i] - parpre[i];
s2 += pappo * pappo;
}
if (s1 == 0.0) {
s1 = 1e-10;
}
if (s2 / s1 <= tol1 * tol1) {
break;
}
}
#ifdef FDEBUG
fprintf(stderr, "\n\n*** Now going to Hessian ***\n\n");
#endif
ittot = it1 + 1;
for (it2=0; it2<100 && !conv; it2++) {
/* compute residuals for covariance matrix */
ll = garch_ll(c, nc, res2, res, yhat, y,
X, nx, t1, t2, param, b, &a0,
q, p, h);
/* store previous coefficients */
for (i=0; i<nparam; i++) {
parpre[i] = param[i];
}
#ifdef FDEBUG
fprintf(stderr, "*** Calling garch_full_hessian, round %d\n", it2);
#endif
garch_full_hessian(t1, t2, X, nx, yhat, c, nc,
res2, res, y, tol2, &it2, vch,
(const double **) g,
grad, param, nparam, b, &a0,
q, p, h, dhdp, zt);
garch_iter_info(ittot++, param, nparam, ll, 1, prn);
s1 = 0.0;
s2 = 0.0;
for (i=0; i<nparam; i++) {
s1 += parpre[i] * parpre[i];
pappo = param[i] - parpre[i];
s2 += pappo * pappo;
}
if (s1 == 0.0) {
s1 = 1e-10;
}
if (s2 / s1 > tol2 * tol2) {
/* not converged yet */
continue;
}
sumgra = 0.0;
for (i=0; i<nparam; i++) {
sumgra += grad[i] * grad[i];
}
if (sumgra >= 1.0e-4) {
pprintf(prn, "\nParameters and gradients at iteration %d:\n\n",
ittot);
for (i=0; i<nparam; i++) {
pprintf(prn, "%12.6f (%9.6f)\n", param[i], grad[i]);
}
pprintf(prn, "\nSum of squared gradients = %.9g (should be less "
"than %g)\n", sumgra, 1.0e-4);
err = E_NOCONV;
goto garch_exit;
} else {
pprintf(prn, "\nFull Hessian convergence at iteration %d, "
"tol = %.9g\n\n", ittot, tol2);
conv = 1;
*iters = ittot;
amax[0] = ll;
}
}
if (!conv && !err) {
err = E_NOCONV;
}
if (!err) {
double sderr;
/* build the desired VCV variant */
err = make_garch_vcv(t1, t2, X, nx,
yhat, y, c, nc,
res2, res,
(const double **) g, grad,
param, nparam,
b, &a0, q, p,
h, dhdp, zt,
(const double *) vch, vcv,
vopt);
if (!err) {
/* transcribe coefficients and standard errors */
for (i=0; i<nparam; i++) {
if (vcv[vix(i,i)] > 0.0) {
sderr = sqrt(vcv[vix(i,i)]);
} else {
sderr = 0.0;
}
amax[i+1] = param[i];
amax[i+1+nparam] = sderr;
}
}
}
garch_exit:
vs_free(dhdp, nparam, g, nc, param, grad, c, aux, vch,
parpre, yhat);
return err;
}
/* Compute the GARCH log-likelihood. Parameters are passed in the
"param" array. res, res2 and ht are be computed here (res2
holds squared residuals).
*/
static double
garch_ll (double *c, int nc, double *res2,
double *res, double *yhat,
const double *y, const double **X, int nx,
int t1, int t2, const double *param, double *b,
double *a0, int q, int p, double *h)
{
int i, t, lag;
int n = t2 - t1 + 1;
double uncvar, ll;
const double *alpha = param + nc + 1;
const double *beta = param + nc + q + 1;
for (i=0; i<nc; i++) {
c[i] = param[i];
}
*a0 = param[nc];
#ifdef FDEBUG
fprintf(stderr, "garch_ll: a0 = %.9g\n", *a0);
for (i=0; i<q; i++) {
fprintf(stderr, " alpha[%d] = %.9g\n", i, alpha[i]);
}
for (i=0; i<p; i++) {
fprintf(stderr, " beta[%d] = %.9g\n", i, beta[i]);
}
#endif
/* Compute residuals, squared residuals, and unconditional
variance over the real estimation period */
copy_coeff(c, nc, b);
uncvar = 0.0;
for (t = t1; t <= t2; t++) {
yhat[t] = get_yhat(X, nx, t, b);
res[t] = y[t] - yhat[t];
res2[t] = res[t] * res[t];
uncvar += res2[t];
}
uncvar /= n;
#ifdef FDEBUG
fprintf(stderr, "uncvar = %.9g (T=%d)\n", uncvar, n);
#endif
/*
We use sample unconditional variance as the starting value
(at time 0, -1, -2, etc.) for the squared residuals and ht;
we use 0 as the starting value for residuals.
*/
lag = (p > q)? p : q;
for (t = t1-lag; t < t1; ++t) {
res[t] = 0.0;
res2[t] = h[t] = uncvar;
}
for (t=t1; t<=t2; t++) {
h[t] = *a0;
for (i=1; i<=q; i++) {
h[t] += res2[t-i] * alpha[i-1];
}
for (i=1; i<=p; i++) {
h[t] += h[t-i] * beta[i-1];
}
/* arbitrary */
if (h[t] <= 0.0) {
h[t] = SMALL_HT;
}
}
ll = 0.0;
for (t=t1; t<=t2; t++) {
double hts = h[t] * gscale * gscale;
ll -= 0.5 * log(hts) + 0.5 * res2[t] / h[t] + LN_SQRT_2_PI;
}
return ll;
}
/*
Check that the values of the parameters of the conditional
variance, ht, are in the set of the admissible values. If a0 is
less or equal than zero it is set to SMALL_HT. If alpha and beta
are less than zero they are set to zero; also if the sum of alpha
and beta is > 1.0, then alpha and beta are normalized (divided by
sum).
*/
static void check_ht (double *b, int np)
{
double sum = 0.0;
int i;
if (b[0] <= 0.0) {
b[0] = SMALL_HT;
}
for (i=1; i<np; i++) {
if (b[i] < 0.0) {
b[i] = 0.0;
}
sum += b[i];
}
if (sum > 1.0) {
for (i=1; i<np; i++) {
b[i] /= sum;
}
}
}
/* combined setup for OP matrix, Information Matrix and Hessian */
static int vcv_setup (int t1, int t2, double *c, int nc,
double *res2, double *res, int *count,
const double **g, double *grad,
double *param, int nparam, double *a0,
int q, int p, double *h,
double **dhdp, double *zt,
double ***H, double *vcv, int code)
{
int i, j, k, t, n, lag, nvparm;
double *asum2;
const double *alpha = param + nc + 1;
const double *beta = param + nc + q + 1;
*a0 = param[nc];
#ifdef FDEBUG
fprintf(stderr, "vcv_setup: a0 = %.9g\n", *a0);
for (i=0; i<q; i++) {
fprintf(stderr, " alpha[%d] = %.9g\n", i, alpha[i]);
}
for (i=0; i<p; i++) {
fprintf(stderr, " beta[%d] = %.9g\n", i, beta[i]);
}
#endif
/* some useful abbreviations */
nvparm = 1 + q + p;
lag = (p > q)? p : q;
n = t2 - t1 + 1;
#ifdef FDEBUG
fprintf(stderr, "make vcv: lag=%d, nc=%d, nparam=%d\n",
lag, nc, nparam);
#endif
asum2 = malloc(nc * sizeof *asum2);
if (asum2 == NULL) {
return 1;
}
if (count != NULL) {
++(*count);
}
for (i=0; i<nc; i++) {
c[i] = param[i];
}
/* Begin computation of dhtdp wrt the variance parameters; we
start computing derivatives of starting values; for ht starting
values are obtained from the unconditional variance of the
residuals.
*/
for (k=1; k<=p; k++) {
for (i=0; i<nvparm; i++) {
dhdp[nc+i][t1-k] = 0.0;
if (H != NULL) { /* hessian only */
for (j=0; j<nvparm; j++) {
H[nc+i][nc+j][k] = 0.0;
}
}
}
}
for (t=t1; t<=t2; t++) {
/* fill in zt at time t (see p. 401) */
zt[0] = 1.0;
for (i=1; i<=q; i++) {
zt[i] = res2[t-i];
}
for (i=1; i<=p; i++) {
zt[q+i] = h[t-i];
}
/* Fill in dhtdp at time t, part relative to variance parameters
(eq. 7, p. 402)
*/
for (i=0; i<nvparm; i++) {
dhdp[nc+i][t] = zt[i];
for (j=1; j<=p; j++) {
dhdp[nc+i][t] += dhdp[nc+i][t-j] * beta[j-1];
}
}
}
/* Build matrix dhtdp, block for regression coefficients
(eq. 13). We use 0 as starting value (time 0, -1, etc.) for the
derivatives of ht wrt the coefficients; we also use 0 as
starting values for the residuals.
*/
/* building blocks for pre-sample terms */
for (i=0; i<nc; i++) {
asum2[i] = 0.0;
for (t=t1; t<=t2; t++) {
asum2[i] -= res[t] * 2.0 * g[i][t];
}
asum2[i] /= n;
}
/* pre-sample range */
for (t=t1-lag; t<t1; t++) {
for (i=0; i<nc; i++) {
dhdp[i][t] = asum2[i];
}
}
/* actual sample range */
for (t=t1; t<=t2; t++) {
for (i=0; i<nc; i++) {
dhdp[i][t] = 0.0;
for (j=1; j<=q; j++) {
if (t - q < t1) {
dhdp[i][t] += alpha[j-1] * asum2[i];
} else {
dhdp[i][t] -= alpha[j-1] * 2.0 * g[i][t-j] * res[t-j];
}
}
for (j=1; j<=p; j++) {
dhdp[i][t] += dhdp[i][t-j] * beta[j-1];
}
}
}
/* Initialize gradient and vcv */
for (i=0; i<nparam; i++) {
grad[i] = 0.0;
for (j=0; j<nparam; j++) {
vcv[vix(i,j)] = 0.0;
}
}
for (t=t1; t<=t2; t++) {
double r_h = res[t] / h[t];
double r2_h = res[t] * r_h;
double aa, bb;
/*
First part, relative to regression coefficients (eq. 10, p. 402)
*/
for (i=0; i<nc; i++) {
aa = r_h * g[i][t] + .5 / h[t] * dhdp[i][t] * (r2_h - 1.0);
grad[i] += aa;
if (code == VCV_OP) {
for (j=0; j<=i; j++) {
bb = r_h * g[j][t] + .5 / h[t] * dhdp[j][t] * (r2_h - 1.0);
vcv[vix(i,j)] += aa * bb;
vcv[vix(j,i)] = vcv[vix(i,j)];
}
for (j=0; j<nvparm; j++) {
int ncj = nc + j;
vcv[vix(i,ncj)] += aa * 0.5 / h[t] * dhdp[ncj][t] * (r2_h - 1.0);
vcv[vix(ncj,i)] = vcv[vix(i,ncj)];
}
}
}
/*
Second part, relative to variance parameters (eq. 6, p. 401)
*/
for (i=0; i<nvparm; i++) {
int nci = nc + i;
aa = .5 / h[t] * dhdp[nci][t] * (r2_h - 1.0);
grad[nci] += aa;
if (code == VCV_OP) {
for (j=0; j<=i; j++) {
int ncj = nc + j;
vcv[vix(nci,ncj)] += aa * 0.5 / h[t] * dhdp[ncj][t] * (r2_h - 1.0);
vcv[vix(ncj,nci)] = vcv[vix(nci,ncj)];
}
}
}
}
if (code == VCV_IM) {
for (t=t1; t<=t2; t++) {
double ht2 = h[t] * h[t];
/* Part relative to the coefficients (eq. 30, p. 406).
Since we take the expected value, only the first two terms
remain.
*/
for (i=0; i<nc; i++) {
for (j=0; j<nc; j++) {
vcv[vix(i,j)] = vcv[vix(i,j)]
- g[i][t] * g[j][t] / h[t]
- .5 * dhdp[i][t] * dhdp[j][t] / ht2;
}
}
/* Part relative to the variance parameters (eq. 29, p. 406).
Since we take the expected value, only the second term
remains.
*/
for (i=nc; i<nparam; i++) {
for (j=nc; j<nparam; j++) {
vcv[vix(i,j)] -= .5 * dhdp[i][t] * dhdp[j][t] / ht2;
}
}
}
}
if (code != VCV_HESSIAN) {
free(asum2);
return 0;
}
/* Initial values in dhdpdp (here, "H") are 2/t x'x, and zero for
off-diagonal (mixed) blocks. */
for (k=0; k<lag; k++) {
for (i=0; i<nc; i++) {
for (j=0; j<nc; j++) {
H[i][j][k+1] = 0.;
}
}
for (t=t1; t<=t2; t++) {
for (i=0; i<nc; i++) {
for (j=0; j<nc; j++) {
H[i][j][k+1] += 2.0 * g[i][t] * g[j][t] / n;
}
}
}
for (i=0; i<nc; i++) {
for (j=0; j<nvparm; j++) {
/* mod. by AC: zero _all_ mixed entries */
H[i][nc+j][k+1] = H[nc+j][i][k+1] = 0.0;
}
}
}
/* Now we fill out the full Hessian */
for (t=t1; t<=t2; ++t) {
double r_h = res[t] / h[t];
double r2_h = r_h * res[t];
double r2_h3 = r2_h / (h[t] * h[t]);
double u_h2 = 1.0 / (h[t] * h[t]);
for (i=0; i<nparam; i++) {
for (j=0; j<nparam; j++) {
H[i][j][0] = 0.0;
}
}
if (lag <= 0) {
goto L90;
}
for (k=1; k<=q; k++) {
for (i=0; i<nc; i++) {
for (j=0; j<nc; j++) {
if (t - q < t1) {
H[i][j][0] += H[i][j][q] * alpha[k-1];
} else {
H[i][j][0] += 2.0 *
g[i][t-k] * g[j][t-k] * alpha[k-1];
}
}
}
}
for (k=1; k<=p; k++) {
for (i=0; i<nc; i++) {
for (j=0; j<nc; j++) {
H[i][j][0] += H[i][j][k] * beta[k-1];
}
}
}
for (i=0; i<nc; i++) {
for (k=1; k<=q; k++) {
if (t - q < t1) {
H[i][nc+k][0] += asum2[i];
} else {
H[i][nc+k][0] -= 2.0 * g[i][t-k] * res[t-k];
}
}
for (k=1; k<=p; k++) {
H[i][nc+q+k][0] += dhdp[i][t-k];
}
}
for (k=1; k<=p; k++) {
for (i=0; i<nc; i++) {
for (j=0; j<nvparm; j++) {
H[i][nc+j][0] += H[i][nc+j][k] * beta[k-1];
}
}
}
L90:
/* Part relative to the coefficients (eq. 15, p. 403).
Since we take the expected value, only the first two terms
remain.
*/
for (i=0; i<nc; i++) {
for (j=0; j<nc; j++) {
vcv[vix(i,j)] = vcv[vix(i,j)]
- g[i][t] * g[j][t] / h[t]
- .5 * r2_h3 * dhdp[i][t] * dhdp[j][t]
- (r_h * g[j][t] * dhdp[i][t]) / h[t]
- (r_h * g[i][t] * dhdp[j][t]) / h[t]
+ 0.5 * (r2_h - 1.0) *
(H[i][j][0] / h[t] - dhdp[i][t]
* dhdp[j][t] / (h[t] * h[t]));
}
}
/* Part relative to the variance parameters (eq. 14, p. 403).
Since we take the expected value, only the second term
remains.
*/
if (p > 0) {
for (i=0; i<nvparm; i++) {
for (j=1; j<=p; j++) {
H[nc+i][nc+q+j][0] += dhdp[nc+i][t-j];
}
}
for (i=1; i<=p; i++) {
for (j=0; j<nvparm; j++) {
H[nc+q+i][nc+j][0] += dhdp[nc+j][t-i];
}
}
for (k=1; k<=p; k++) {
for (i=0; i<nvparm; i++) {
for (j=0; j<nvparm; j++) {
H[nc+i][nc+j][0] += H[nc+i][nc+j][k] * beta[k-1];
#ifdef FDEBUG
if ((t==t1 || t==t2)) {
fprintf(stderr, "Set H[%d][%d][0] = %.9g\n",
nc+i, nc+j, H[nc+i][nc+j][0]);
fprintf(stderr, "incr = %.9g * %.9g\n",
H[nc+i][nc+j][k], beta[k-1]);
}
#endif
}
}
}
}
for (i=nc; i<nparam; i++) {
for (j=nc; j<nparam; j++) {
vcv[vix(i,j)] = vcv[vix(i,j)]
+ .5 * u_h2 * dhdp[i][t] * dhdp[j][t]
- r2_h3 * dhdp[i][t] * dhdp[j][t]
+ .5 * (r2_h - 1.0) / h[t] * H[i][j][0];
#ifdef FDEBUG
if ((t==t1 || t==t2)) {
fprintf(stderr, "Set vcv[%d] = %.9g (using Hval %.9g)\n",
vix(i,j), vcv[vix(i,j)], H[i][j][0]);
}
#endif
}
}
/* top-right mixed part (eq. 17, p. 403) */
for (i=0; i<nc; i++) {
for (j=0; j<nvparm; j++) {
vcv[vix(i,nc+j)] = vcv[vix(i,nc+j)]
- g[i][t] * r_h * dhdp[nc+j][t] / h[t]
- .5 * (r2_h - 1.0) * dhdp[nc+j][t] * dhdp[i][t] / (h[t] * h[t])
+ .5 * (r2_h - 1.0) * H[i][nc+j][0] / h[t]
- .5 * r2_h * u_h2 * dhdp[i][t] * dhdp[nc+j][t];
/* and bottom left too */
vcv[vix(nc+j, i)] = vcv[vix(i, nc+j)];
}
}
/* before quitting time t, tidy up dhdpdp */
for (k=0; k<lag; k++) {
for (i=0; i<nparam; i++) {
for (j=0; j<nparam; j++) {
#ifdef FDEBUG
if (t < 5) {
fprintf(stderr, "t=%d: setting H(%d,%d,%d) to H(%d,%d,%d)=%g\n",
t,i,j,lag-k,i,j,lag-k-1, H[i][j][lag-k-1]);
}
#endif
H[i][j][lag-k] = H[i][j][lag-k-1];
}
}
}
}
free(asum2);
return 0;
} /* vcv_setup */
/* block-diagonal information matrix.
parameters are passed in the array "c"
*/
static int
garch_info_matrix (int t1, int t2,
const double **X, int nx, double *yhat, double *c,
int nc, double *res2, double *res, const double *y,
double toler, int *count, double *vcv,
const double **g, double *grad,
double *param, int nparam, double *b,
double *a0, int q, int p,
double *h, double **dhdp, double *zt)
{
int i, j;
double d, d0, d1, d2, ll2, d3;
double ll3, d12, d31, d23, dd;
double di, ff, dm, ds;
int iv, nexp;
double a1s, a2s, a3s;
int it1, it2, it3, it4, it5;
double d12s, dac, dub, d23s, d31s;
int err;
double s1, s2, stre, bigd;
double cappa;
int ncall;
int nvparm;
double oldstp = 9.0e+39;
static double ll1 = 0.0, fs = 0.0;
double *gg = NULL, *step = NULL;
iv = 0;
it1 = it2 = it3 = it4 = it5 = 0;
nvparm = q + p + 1;
/* calculate information matrix */
vcv_setup(t1, t2, c, nc, res2, res,
count, g, grad, param, nparam,
a0, q, p, h,
dhdp, zt, NULL, vcv,
VCV_IM);
/* invert the information matrix */
err = invert(vcv, nparam);
if (err) {
fprintf(stderr, "garch_info_matrix: matrix inversion failed\n");
return 1;
}
if (count == NULL) {
/* just calculating vcv at convergence */
goto vcv_exit;
}
gg = malloc(nparam * sizeof *gg);
if (gg == NULL) {
return 1;
}
step = malloc(nparam * sizeof *step);
if (step == NULL) {
free(gg);
return 1;
}
/* here we start iteration */
/* calculate the step for the new coefficients */
s1 = s2 = 0.0;
for (i=0; i<nparam; i++) {
s1 += param[i] * param[i];
gg[i] = param[i];
step[i] = 0.0;
for (j=0; j<nparam; j++) {
step[i] -= vcv[vix(i,j)] * grad[j];
}
s2 += step[i] * step[i];
}
if (s1 == 0.0) {
s1 = 1e-10;
}
stre = s2 / s1;
s2 = sqrt(s2);
stre = sqrt(stre);
oldstp = s2;
#ifdef FDEBUG
fprintf(stderr, "s2 = %.9g\n", s2);
#endif
for (i=0; i<nparam; i++) {
step[i] /= s2;
}
it4 += iv;
ds = s2;
if (stre <= toler) {
goto L496;
}
ncall = 0;
nexp = *count / 5;
if (nexp > 5) {
nexp = 5;
}
cappa = pow(2.0, nexp);
d0 = s2;
d0 /= cappa;
dac = d0 * .001;
dub = d0 * 4.0;
if (*count == 1) {
ll1 = -garch_ll(c, nc, res2, res, yhat, y,
X, nx, t1, t2, param, b, a0,
q, p, h);
#ifdef FFDEBUG
fprintf(stderr, "count=1, ll1=%.9g\n", ll1);
#endif
}
for (i=0; i<nparam; i++) {
param[i] = gg[i] + step[i] * d0;
}
check_ht(param + nc, nvparm);
#ifdef FDEBUG
for (i=0; i<nparam; i++) {
fprintf(stderr, "param[%d] in matinf(1) = %.9g\n", i, param[i]);
}
#endif
ll2 = -garch_ll(c, nc, res2, res, yhat, y,
X, nx, t1, t2, param, b, a0,
q, p, h);
#ifdef FDEBUG
fprintf(stderr, "ll2=%.9g, ll1=%.9g\n", ll2, ll1);
#endif
if (ll2 > ll1) {
goto L307;
}
d1 = 0.0;
d2 = d0;
d3 = d0 + d0;
for (i=0; i<nparam; i++) {
param[i] = gg[i] + step[i] * d3;
}
check_ht(param + nc, nvparm);
#ifdef FDEBUG
for (i=0; i<nparam; i++) {
fprintf(stderr, "param[%d] in matinf(2) = %.9g\n", i, param[i]);
}
#endif
ll3 = -garch_ll(c, nc, res2, res, yhat, y,
X, nx, t1, t2, param, b, a0,
q, p, h);
goto L325;
L307:
d1 = -d0;
d2 = 0.0;
d3 = d0;
ll3 = ll2;
ll2 = ll1;
for (i=0; i<nparam; i++) {
param[i] = gg[i] + step[i] * d1;
}
check_ht(param + nc, nvparm);
#ifdef FDEBUG
for (i=0; i<nparam; i++) {
fprintf(stderr, "param[%d] in matinf(3) = %.9g\n", i, param[i]);
}
#endif
ll1 = -garch_ll(c, nc, res2, res, yhat, y,
X, nx, t1, t2, param, b, a0,
q, p, h);
L325:
d23 = d2 - d3;
d31 = d3 - d1;
d12 = d1 - d2;
di = d23 * ll1 + d31 * ll2 + d12 * ll3;
bigd = di * -2.0 / (d23 * d31 * d12);
if (bigd > 0.0) {
goto L400;
}
if (ll3 <= ll1) {
goto L341;
}
L329:
d3 = d2;
ll3 = ll2;
d2 = d1;
ll2 = ll1;
d1 -= dub;
for (i=0; i<nparam; ++i) {
param[i] = gg[i] + d1 * step[i];
}
check_ht(param + nc, nvparm);
#ifdef FDEBUG
for (i=0; i<nparam; i++) {
fprintf(stderr, "param[%d] in matinf(4) = %.9g\n", i, param[i]);
}
#endif
ll1 = -garch_ll(c, nc, res2, res, yhat, y,
X, nx, t1, t2, param, b, a0,
q, p, h);
if (++ncall > 100) {
goto L490;
}
goto L325;
L341:
d1 = d2;
ll1 = ll2;
d2 = d3;
ll2 = ll3;
d3 += dub;
for (i=0; i<nparam; i++) {
param[i] = gg[i] + d3 * step[i];
}
check_ht(param + nc, nvparm);
#ifdef FDEBUG
for (i=0; i<nparam; i++) {
fprintf(stderr, "param[%d] in matinf(5) = %.9g\n", i, param[i]);
}
#endif
ll3 = -garch_ll(c, nc, res2, res, yhat, y,
X, nx, t1, t2, param, b, a0,
q, p, h);
if (++ncall > 100) {
goto L490;
}
goto L325;
L400:
d23s = d23 * (d2 + d3);
d31s = d31 * (d3 + d1);
d12s = d12 * (d1 + d2);
ds = (d23s * ll1 + d31s * ll2 + d12s * ll3) * .5 / di;
for (i=0; i<nparam; i++) {
param[i] = gg[i] + step[i] * ds;
}
check_ht(param + nc, nvparm);
#ifdef FDEBUG
for (i=0; i<nparam; i++) {
fprintf(stderr, "param[%d] in matinf(6) = %.9g\n", i, param[i]);
}
#endif
fs = -garch_ll(c, nc, res2, res, yhat, y,
X, nx, t1, t2, param, b, a0,
q, p, h);
if (++ncall > 100) {
goto L490;
}
a1s = (d = d1 - ds, fabs(d));
a2s = (d = d2 - ds, fabs(d));
a3s = (d = d3 - ds, fabs(d));
dm = a1s;
if (a3s < dm) {
dm = a3s;
}
if (dub >= dm) {
goto L422;
}
if (ds < d1 - dub) {
goto L329;
}
if (ds > d3 + dub) {
goto L341;
}
L422:
if (a1s < dac || a2s < dac || a3s < dac) {
goto L490;
}
if (ll1 < ll2 || ll1 < ll3) {
goto L434;
}
d1 = ds;
ll1 = fs;
goto L459;
L434:
if (ll2 < ll3 || ll2 < ll1) {
goto L447;
}
d2 = ds;
ll2 = fs;
goto L459;
L447:
d3 = ds;
ll3 = fs;
L459:
if (d2 <= d3) {
goto L463;
}
dd = d2;
ff = ll2;
d2 = d3;
ll2 = ll3;
d3 = dd;
ll3 = ff;
L463:
if (d1 <= d2) {
goto L325;
}
dd = d1;
ff = ll1;
d1 = d2;
ll1 = ll2;
d2 = dd;
ll2 = ff;
goto L459;
L490:
if (fs <= ll1) {
goto L491;
}
fs = ll1;
ds = d1;
L491:
if (fs <= ll2) {
goto L492;
}
fs = ll2;
ds = d2;
L492:
if (fs <= ll3) {
goto L496;
}
fs = ll3;
ds = d3;
L496:
for (i=0; i<nparam; i++) {
param[i] = gg[i] + ds * step[i];
}
check_ht(param + nc, nvparm);
#ifdef FDEBUG
for (i=0; i<nparam; i++) {
fprintf(stderr, "param[%d] in matinf(7) = %.9g\n", i, param[i]);
}
#endif
ll1 = fs;
fs = -fs;
it5 += iv;
vcv_exit:
/* change the sign of the matrix */
for (i=0; i<nparam; i++) {
for (j=0; j<nparam; j++) {
vcv[vix(i,j)] *= -1.0;
}
}
free(gg);
free(step);
return 0;
} /* garch_info_matrix */
static void free_dhdpdp (double ***H, int np)
{
int i, j;
for (i=0; i<np; i++) {
for (j=0; j<np; j++) {
free(H[i][j]);
}
free(H[i]);
}
free(H);
}
static double ***allocate_dhdpdp (int np, int p, int q)
{
double ***H;
int i, j, lag = (p > q)? p : q;
H = malloc(np * sizeof *H);
if (H == NULL) {
return NULL;
}
for (i=0; i<np; i++) {
H[i] = malloc(np * sizeof **H);
if (H[i] == NULL) {
return NULL;
}
for (j=0; j<np; j++) {
H[i][j] = malloc((lag + 1) * sizeof ***H);
if (H[i][j] == NULL) {
return NULL;
}
}
}
return H;
}
static int
garch_full_hessian (int t1, int t2,
const double **X, int nx, double *yhat, double *c,
int nc, double *res2, double *res, const double *y,
double toler, int *count, double *vcv,
const double **g, double *grad,
double *param, int nparam, double *b, double *a0,
int q, int p, double *h, double **dhdp, double *zt)
{
int i, j;
double d;
double d0, d1, d2, ll2, d3;
double ll3, d12, d31, d23, dd, di, ff;
int iv, ncall, nexp, nvparm;
double dm, ds;
double a1s, a2s, a3s;
int it1, it2, it3, it4, it5;
double dac, d12s, dub, d23s, d31s;
int err;
double s1, s2, stre, bigd;
double cappa;
int lag;
static double ll1 = 0.0, fs = 0.0;
double oldstp = 9.0e+39;
double *gg, *step;
double ***H;
gg = malloc(nparam * sizeof *gg);
if (gg == NULL) {
return 1;
}
step = malloc(nparam * sizeof *step);
if (step == NULL) {
free(gg);
return 1;
}
/* 3rd dimension of dhdpdp is max(p,q) + 1 */
H = allocate_dhdpdp(nparam, p, q);
if (H == NULL) {
free(gg);
free(step);
return 1;
}
iv = 0;
it1 = it2 = it3 = it4 = it5 = 0;
nvparm = q + p + 1;
lag = (p > q)? p : q;
/* calculate the full Hessian */
vcv_setup(t1, t2, c, nc, res2, res,
count, g, grad, param, nparam,
a0, q, p, h,
dhdp, zt, H, vcv,
VCV_HESSIAN);
/* invert the Hessian */
err = invert(vcv, nparam);
if (err) {
fprintf(stderr, "garch_full_hessian: matrix inversion failed\n");
}
/* Start iteration here */
/* calculate the step for the new coefficients */
s2 = 0.0;
for (i=0; i<nparam; i++) {
gg[i] = param[i];
step[i] = 0.0;
for (j=0; j<nparam; j++) {
step[i] -= vcv[vix(i,j)] * grad[j];
}
s2 += step[i] * step[i];
}
s1 = 0.0;
for (i=0; i<nparam; i++) {
s1 += param[i] * param[i];
}
if (s1 == 0.0) {
s1 = 1e-10;
}
stre = s2 / s1;
s2 = sqrt(s2);
stre = sqrt(stre);
oldstp = s2;
for (i=0; i<nparam; i++) {
step[i] /= s2;
}
it4 += iv;
ds = s2;
if (stre <= toler) {
goto L496;
}
ncall = 0;
nexp = *count / 5;
if (nexp > 5) {
nexp = 5;
}
cappa = pow(2.0, nexp);
d0 = s2;
d0 /= cappa;
dac = d0 * .001;
dub = d0 * 4.0;
if (*count == 1) {
ll1 = -garch_ll(c, nc, res2, res, yhat, y,
X, nx, t1, t2, param, b, a0,
q, p, h);
#ifdef FDEBUG
fprintf(stderr, "hess: ll1 = %.9g\n", ll1);
#endif
}
for (i=0; i<nparam; i++) {
param[i] = gg[i] + step[i] * d0;
}
check_ht(param + nc, nvparm);
ll2 = -garch_ll(c, nc, res2, res, yhat, y,
X, nx, t1, t2, param, b, a0,
q, p, h);
#ifdef FDEBUG
fprintf(stderr, "hess: ll2 = %.9g\n", ll2);
#endif
if (ll2 > ll1) {
goto L307;
}
d1 = 0.0;
d2 = d0;
d3 = d0 + d0;
for (i=0; i<nparam; ++i) {
param[i] = gg[i] + step[i] * d3;
}
check_ht(param + nc, nvparm);
ll3 = -garch_ll(c, nc, res2, res, yhat, y,
X, nx, t1, t2, param, b, a0,
q, p, h);
#ifdef FDEBUG
fprintf(stderr, "hess: ll3 = %.9g\n", ll3);
#endif
goto L325;
L307:
d1 = -d0;
d2 = 0.0;
d3 = d0;
ll3 = ll2;
ll2 = ll1;
for (i=0; i<nparam; i++) {
param[i] = gg[i] + step[i] * d1;
}
check_ht(param + nc, nvparm);
ll1 = -garch_ll(c, nc, res2, res, yhat, y,
X, nx, t1, t2, param, b, a0,
q, p, h);
#ifdef FDEBUG
fprintf(stderr, "hess: ll1(2) = %.9g\n", ll1);
#endif
L325:
d23 = d2 - d3;
d31 = d3 - d1;
d12 = d1 - d2;
di = d23 * ll1 + d31 * ll2 + d12 * ll3;
bigd = di * -2.0 / (d23 * d31 * d12);
#ifdef FDEBUG
fprintf(stderr, "hess: bigd = %.9g\n", bigd);
#endif
if (bigd > 0.0) {
goto L400;
}
if (ll3 <= ll1) {
goto L341;
}
L329:
d3 = d2;
ll3 = ll2;
d2 = d1;
ll2 = ll1;
d1 -= dub;
for (i=0; i<nparam; i++) {
param[i] = gg[i] + step[i] * d1;
}
check_ht(param + nc, nvparm);
ll1 = -garch_ll(c, nc, res2, res, yhat, y,
X, nx, t1, t2, param, b, a0,
q, p, h);
#ifdef FDEBUG
fprintf(stderr, "hess: ll1(3) = %.9g\n", ll1);
#endif
if (++ncall > 100) {
goto L490;
}
goto L325;
L341:
d1 = d2;
ll1 = ll2;
d2 = d3;
ll2 = ll3;
d3 += dub;
for (i=0; i<nparam; i++) {
param[i] = gg[i] + step[i] * d3;
}
check_ht(param + nc, nvparm);
ll3 = -garch_ll(c, nc, res2, res, yhat, y,
X, nx, t1, t2, param, b, a0,
q, p, h);
#ifdef FDEBUG
fprintf(stderr, "hess: ll3(2) = %.9g\n", ll3);
#endif
if (++ncall > 100) {
goto L490;
}
goto L325;
L400:
d23s = d23 * (d2 + d3);
d31s = d31 * (d3 + d1);
d12s = d12 * (d1 + d2);
ds = (d23s * ll1 + d31s * ll2 + d12s * ll3) * .5 / di;
for (i=0; i<nparam; i++) {
param[i] = gg[i] + step[i] * ds;
}
check_ht(param + nc, nvparm);
fs = -garch_ll(c, nc, res2, res, yhat, y,
X, nx, t1, t2, param, b, a0,
q, p, h);
#ifdef FDEBUG
fprintf(stderr, "hess: fs = %.9g, ds = %.9g\n", fs, ds);
#endif
if (++ncall > 100) {
goto L490;
}
a1s = (d = d1 - ds, fabs(d));
a2s = (d = d2 - ds, fabs(d));
a3s = (d = d3 - ds, fabs(d));
dm = a1s;
if (a3s < dm) {
dm = a3s;
}
if (dub >= dm) {
goto L422;
}
if (ds < d1 - dub) {
goto L329;
}
if (ds > d3 + dub) {
goto L341;
}
L422:
if (a1s < dac || a2s < dac || a3s < dac) {
goto L490;
}
if (ll1 < ll2 || ll1 < ll3) {
goto L434;
}
d1 = ds;
ll1 = fs;
goto L459;
L434:
if (ll2 < ll3 || ll2 < ll1) {
goto L447;
}
d2 = ds;
ll2 = fs;
goto L459;
L447:
d3 = ds;
ll3 = fs;
L459:
if (d2 <= d3) {
goto L463;
}
dd = d2;
ff = ll2;
d2 = d3;
ll2 = ll3;
d3 = dd;
ll3 = ff;
L463:
if (d1 <= d2) {
goto L325;
}
dd = d1;
ff = ll1;
d1 = d2;
ll1 = ll2;
d2 = dd;
ll2 = ff;
goto L459;
L490:
if (fs <= ll1) {
goto L491;
}
fs = ll1;
ds = d1;
L491:
if (fs <= ll2) {
goto L492;
}
fs = ll2;
ds = d2;
L492:
if (fs <= ll3) {
goto L496;
}
fs = ll3;
ds = d3;
L496:
for (i=0; i<nparam; i++) {
param[i] = gg[i] + step[i] * ds;
}
check_ht(param + nc, nvparm);
ll1 = fs;
fs = -fs;
it5 += iv;
/* change the sign of the inverse */
for (i=0; i<nparam; i++) {
for (j=0; j<nparam; j++) {
vcv[vix(i,j)] *= -1.0;
}
}
free(gg);
free(step);
free_dhdpdp(H, nparam);
return 0;
} /* garch_full_hessian */
static int invert (double *g, int dim)
{
integer m = dim;
integer info, lwork;
integer *ipiv;
double *work;
#ifdef FDEBUG
int i, j;
static int k;
fprintf(stderr, "Matrix inversion %d, dim=%d\n", k, dim);
for (i=0; i<dim; i++) {
for (j=0; j<dim; j++) {
fprintf(stderr, "%14.9g ", g[i + dim * j]);
}
fputc('\n', stderr);
}
k++;
#endif
ipiv = malloc(m * sizeof *ipiv);
if (ipiv == NULL) {
return 1;
}
work = malloc(sizeof *work);
if (work == NULL) {
free(ipiv);
return 1;
}
dgetrf_(&m, &m, g, &m, ipiv, &info);
if (info != 0) {
free(ipiv);
free(work);
return 1;
}
lwork = -1;
dgetri_(&m, g, &m, ipiv, work, &lwork, &info);
if (info != 0 || work[0] <= 0.0) {
free(ipiv);
free(work);
return 1;
}
lwork = (integer) work[0];
work = realloc(work, lwork * sizeof *work);
if (work == NULL) {
free(ipiv);
return 1;
}
dgetri_(&m, g, &m, ipiv, work, &lwork, &info);
free(work);
free(ipiv);
return (int) info;
}
syntax highlighted by Code2HTML, v. 0.9.1