/*
* Copyright (c) by 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 program 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 program; if not, write to the Free Software
* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*
*/
#include "libgretl.h"
#include "pvalues.h"
#include "gretl_matrix.h"
#include "var.h"
#define JDEBUG 0
#define OLDEIG 0
/*
Critical values for Johansen's likelihood ratio tests
are computed using J. Doornik's gamma approximation
*/
/* Matrices for the trace test */
const double s_mTrace_m_coef[5][6] = {
/* n^2 n 1 n==1 n==2 n^1/2 */
{2, -1.00, 0.07, 0.07, 0, 0},
{2, 2.01, 0, 0.06, 0.05, 0},
{2, 1.05, -1.55, -0.50, -0.23, 0},
{2, 4.05, 0.50, -0.23, -0.07, 0},
{2, 2.85, -5.10, -0.10, -0.06, 1.35}
};
const double s_mTrace_v_coef[5][6] = {
{3, -0.33, -0.55, 0, 0, 0},
{3, 3.6, 0.75, -0.4, -0.3, 0},
{3, 1.8, 0, -2.8, -1.1, 0},
{3, 5.7, 3.2, -1.3, -0.5, 0},
{3, 4.0, 0.8, -5.8, -2.66, 0}
};
const double s_mTrace_m_time[5][7] = {
/* sqrt(n)/T n/T n^2/T^2 n==1/T n==1 n==2 n==3 */
{-0.101, 0.499, 0.896, -0.562, 0.00229, 0.00662, 0},
{ 0, 0.465, 0.984, -0.273,-0.00244, 0, 0},
{ 0.134, 0.422, 1.02, 2.17,-0.00182, 0,-0.00321},
{0.0252, 0.448, 1.09, -0.353, 0, 0, 0},
{-0.819, 0.615, 0.896, 2.43, 0.00149, 0, 0}
};
const double s_mTrace_v_time[5][7] = {
{-0.204, 0.980, 3.11, -2.14, 0.0499, -0.0103,-0.00902},
{ 0.224, 0.863, 3.38, -0.807, 0, 0, -0.0091},
{ 0.422, 0.734, 3.76, 4.32,-0.00606, 0,-0.00718},
{ 0, 0.836, 3.99, -1.33,-0.00298,-0.00139,-0.00268},
{ -1.29, 1.01, 3.92, 4.67, 0.00484,-0.00127, -0.0199}
};
/* Matrices for the lambdamax test */
const double s_mMaxev_m_coef[5][5] = {
/* n 1 n==1 n==2 n^1/2 */
{6.0019, -2.7558, 0.67185, 0.11490, -2.7764},
{5.9498, 0.43402, 0.048360, 0.018198, -2.3669},
{5.8271, -1.6487, -1.6118, -0.25949, -1.5666},
{5.8658, 2.5595, -0.34443, -0.077991, -1.7552},
{5.6364, -0.90531, -3.5166, -0.47966, -0.21447}
};
const double s_mMaxev_v_coef[5][5] = {
{1.8806, -15.499, 1.1136, 0.070508, 14.714},
{2.2231, -7.9064, 0.58592, -0.034324, 12.058},
{2.0785, -9.7846, -3.3680, -0.24528, 13.074},
{1.9955, -5.5428, 1.2425, 0.41949, 12.841},
{2.0899, -5.3303, -7.1523, -0.25260, 12.393}
};
static int
gamma_par_asymp (double tracetest, double lmaxtest, JohansenCode det,
int N, double *pval)
{
/*
Asymptotic critical values for Johansen's LR tests via gamma approximation
params:
tracetest, lmaxtest: trace and lambdamax est. statistics
det: index of setup of deterministic regressors
J_NO_CONST = no constant
J_REST_CONST = restricted constant
J_UNREST_CONST = unrestricted constant
J_REST_TREND = restricted trend
J_UNREST_TREND = unrestricted trend
N: cointegration rank under H0;
pval: on output, array of pvalues, for the two tests;
*/
double mt, vt, ml, vl;
const double *tracem, *tracev, *lmaxm, *lmaxv;
double g, x[7];
int i;
tracem = s_mTrace_m_coef[det];
tracev = s_mTrace_v_coef[det];
lmaxm = s_mMaxev_m_coef[det];
lmaxv = s_mMaxev_v_coef[det];
mt = vt = 0.0;
ml = vl = 0.0;
x[0] = N * N;
x[1] = N;
x[2] = 1.0;
x[3] = (N == 1)? 1.0 : 0.0;
x[4] = (N == 2)? 1.0 : 0.0;
x[5] = sqrt((double) N);
for (i=0; i<6; i++) {
mt += x[i] * tracem[i];
vt += x[i] * tracev[i];
if (i) {
ml += x[i] * lmaxm[i-1];
vl += x[i] * lmaxv[i-1];
}
}
g = gamma_cdf_comp(mt, vt, tracetest, 2);
if (na(g)) {
pval[0] = NADBL;
} else {
pval[0] = 1.0 - g;
if (pval[0] < 0.0) {
pval[0] = 0.0;
}
}
g = gamma_cdf_comp(ml, vl, lmaxtest, 2);
if (na(g)) {
pval[1] = NADBL;
} else {
pval[1] = 1.0 - g;
if (pval[1] < 0.0) {
pval[1] = 0.0;
}
}
return 0;
}
/* Remove a possible excess zero from the end of a floating point
number printed to the given precision p: are we working around
a bug in the C library?
*/
static void fix_xstr (char *s, int p)
{
int n = strlen(s);
if (n > p && strspn(s + n - p, "0") == p) {
s[n-1] = 0;
}
}
enum {
PRINT_ALPHA,
PRINT_BETA
};
/* for cointegration test: print cointegrating vectors or adjustments,
either "raw" or re-scaled */
static void print_beta_or_alpha (JohansenInfo *jv, int k,
const DATAINFO *pdinfo, PRN *prn,
int which, int rescale)
{
gretl_matrix *c = (which == PRINT_BETA)? jv->Beta : jv->Alpha;
int rows = gretl_matrix_rows(c);
char xstr[32];
int i, j;
double x, y;
if (rescale) {
pprintf(prn, "\n%s\n", (which == PRINT_BETA)?
_("renormalized beta") :
_("renormalized alpha"));
} else {
pprintf(prn, "\n%s\n", (which == PRINT_BETA)?
_("beta (cointegrating vectors)") :
_("alpha (adjustment vectors)"));
}
for (i=0; i<rows; i++) {
if (i < jv->list[0]) {
pprintf(prn, "%-10s", pdinfo->varname[jv->list[i+1]]);
} else if (jv->code == J_REST_CONST) {
pprintf(prn, "%-10s", "const");
} else if (jv->code == J_REST_TREND) {
pprintf(prn, "%-10s", "trend");
}
for (j=0; j<k; j++) {
x = gretl_matrix_get(c, i, j);
if (rescale) {
y = gretl_matrix_get(jv->Beta, j, j);
if (which == PRINT_BETA) {
x /= y;
} else {
x *= y;
}
}
if (x == -0.0) {
x = 0.0;
}
sprintf(xstr, "%#.5g", x);
fix_xstr(xstr, 5);
pprintf(prn, "%12s ", xstr);
}
pputc(prn, '\n');
}
}
/* Calculate \alpha (adjustments) matrix as per Johansen, 1991, eqn
2.8, p. 1554. Required for the cointegration test, but not
needed when doing a VECM (in which case we get \alpha via
build_VECM_models() below).
*/
static int compute_alpha (JohansenInfo *jv, int n)
{
gretl_matrix *alpha = NULL;
gretl_matrix *tmp1 = NULL;
gretl_matrix *tmp2 = NULL;
int h, nv = gretl_matrix_rows(jv->Svv);
int err = 0;
if (jv->rank > 0) {
h = jv->rank;
} else {
h = nv;
}
tmp1 = gretl_matrix_alloc(nv, h);
tmp2 = gretl_matrix_alloc(h, h);
alpha = gretl_matrix_alloc(n, h);
if (tmp1 == NULL || tmp2 == NULL || alpha == NULL) {
err = E_ALLOC;
}
if (!err) {
gretl_matrix_qform(jv->Beta, GRETL_MOD_TRANSPOSE, jv->Svv,
tmp2, GRETL_MOD_NONE);
err = gretl_invert_symmetric_matrix(tmp2);
}
if (!err) {
gretl_matrix_multiply(jv->Beta, tmp2, tmp1);
gretl_matrix_multiply(jv->Suv, tmp1, alpha);
}
gretl_matrix_free(tmp1);
gretl_matrix_free(tmp2);
if (!err) {
jv->Alpha = alpha;
} else {
gretl_matrix_free(alpha);
}
return err;
}
/* For cointegration test, print the "long-run matrix," \alpha \beta'
in Johansen's notation, or \Zeta_0 in Hamilton's.
*/
static void print_lr_matrix (JohansenInfo *jv, gretl_matrix *Zeta_0,
int neqns, const DATAINFO *pdinfo, PRN *prn)
{
int cols = gretl_matrix_rows(jv->Svv);
int i, j;
pprintf(prn, "%s\n", _("long-run matrix (alpha * beta')"));
pprintf(prn, "%22s", pdinfo->varname[jv->list[1]]); /* N.B. */
for (j=2; j<=jv->list[0]; j++) {
pprintf(prn, "%13s", pdinfo->varname[jv->list[j]]);
}
if (jv->code == J_REST_CONST) {
pprintf(prn, "%13s", "const");
} else if (jv->code == J_REST_TREND) {
pprintf(prn, "%13s", "trend");
}
pputc(prn, '\n');
for (i=0; i<neqns; i++) {
pprintf(prn, "%-10s", pdinfo->varname[jv->list[i+1]]);
for (j=0; j<cols; j++) {
pprintf(prn, "%#12.5g ", gretl_matrix_get(Zeta_0, i, j));
}
pputc(prn, '\n');
}
pputc(prn, '\n');
}
/* Compute Hamilton's Omega (Johansen 1991 calls it Lambda): the
cross-equation variance matrix.
*/
static int compute_omega (GRETL_VAR *vecm)
{
if (vecm->S == NULL) {
vecm->S = gretl_matrix_alloc(vecm->neqns, vecm->neqns);
}
if (vecm->S == NULL) {
return E_ALLOC;
}
gretl_matrix_multiply_mod(vecm->E, GRETL_MOD_TRANSPOSE,
vecm->E, GRETL_MOD_NONE,
vecm->S, GRETL_MOD_NONE);
gretl_matrix_divide_by_scalar(vecm->S, vecm->T);
return 0;
}
static void gretl_matrix_I (gretl_matrix *A, int n)
{
int i;
gretl_matrix_zero(A);
for (i=0; i<n; i++) {
gretl_matrix_set(A, i, i, 1.0);
}
}
/* compute the EC terms and add them to the dataset, Z, so we
can estimate the VECM by OLS (conditional on \beta) */
static int
add_EC_terms_to_dataset (GRETL_VAR *vecm, double ***pZ, DATAINFO *pdinfo,
int iter)
{
const gretl_matrix *Beta = vecm->jinfo->Beta;
int rank = jrank(vecm);
int *list = vecm->jinfo->list;
double xt, bxt, sb;
int i, j, t, v = pdinfo->v;
int id = gretl_VECM_id(vecm);
int vj, err = 0;
if (iter == 0) {
err = dataset_add_series(rank, pZ, pdinfo);
}
if (!err) {
char vname[VNAMELEN];
for (j=0; j<rank; j++) {
sprintf(vname, "EC%d", j + 1);
if (iter > 0) {
/* series already allocated */
vj = varindex(pdinfo, vname);
} else {
vj = v + j;
strcpy(pdinfo->varname[vj], vname);
make_varname_unique(pdinfo->varname[vj], vj, pdinfo);
sprintf(VARLABEL(pdinfo, vj), "error correction term %d from VECM %d",
j + 1, id);
}
for (t=0; t<pdinfo->n; t++) {
if (t < vecm->t1 || t > vecm->t2) {
(*pZ)[vj][t] = NADBL;
} else {
bxt = 0.0;
/* beta * X(t-1) */
for (i=0; i<vecm->neqns; i++) {
xt = (*pZ)[list[i+1]][t-1];
sb = gretl_matrix_get(Beta, i, j);
sb /= gretl_matrix_get(Beta, j, j);
bxt += sb * xt;
}
/* restricted const or trend */
if (restricted(vecm)) {
sb = gretl_matrix_get(Beta, i, j);
sb /= gretl_matrix_get(Beta, j, j);
if (jcode(vecm) == J_REST_TREND) {
sb *= t;
}
bxt += sb;
}
(*pZ)[vj][t] = bxt;
}
}
}
}
return err;
}
/* After doing OLS estimation of the VECM conditional on \beta: copy
the coefficients on the lagged differences (i.e. form the \Gamma
matrices) so we can compute the VAR representation */
static void copy_coeffs_to_Gamma (MODEL *pmod, int i, gretl_matrix **G,
int maxlag, int nv)
{
int j, k, h;
double x;
for (k=0; k<maxlag; k++) {
h = k + pmod->ifc;
/* successive lags (distinct \Gamma_i matrices) */
for (j=0; j<nv; j++) {
/* successive \Delta x_j */
x = pmod->coeff[h];
gretl_matrix_set(G[k], i, j, x);
h += maxlag;
}
}
}
/* Again, after doing OLS estimation of the VECM conditional on \beta:
copy the coefficients on the EC terms (\beta' X) into the \alpha
matrix. */
static void copy_coeffs_to_Alpha (GRETL_VAR *vecm, int i, gretl_matrix *Alpha,
int maxlag)
{
double x;
const MODEL *pmod = vecm->models[i];
/* position in coeff array of first \alpha term */
int base = vecm->jinfo->nexo + gretl_matrix_rows(Alpha) * maxlag;
int j;
for (j=0; j<vecm->jinfo->rank; j++) {
x = pmod->coeff[base + j];
gretl_matrix_set(Alpha, i, j, x);
}
}
/* Form the matrix \Pi = \alpha \beta': since \beta is augmented
in the case of restricted constant or restricted trend, we
may have to make a reduced copy */
static int form_Pi (GRETL_VAR *vecm, const gretl_matrix *Alpha,
gretl_matrix *Pi)
{
gretl_matrix *Beta = vecm->jinfo->Beta;
int err = 0, freeit = 0;
if (gretl_matrix_rows(Beta) > vecm->neqns) {
Beta = gretl_matrix_alloc(vecm->neqns, vecm->jinfo->rank);
if (Beta == NULL) {
err = E_ALLOC;
} else {
double x;
int i, j;
for (i=0; i<vecm->neqns; i++) {
for (j=0; j<vecm->jinfo->rank; j++) {
x = gretl_matrix_get(vecm->jinfo->Beta, i, j);
gretl_matrix_set(Beta, i, j, x);
}
}
freeit = 1;
}
}
if (!err) {
gretl_matrix_multiply_mod(Alpha, GRETL_MOD_NONE,
Beta, GRETL_MOD_TRANSPOSE,
Pi, GRETL_MOD_NONE);
}
if (freeit) {
gretl_matrix_free(Beta);
}
return err;
}
/* VAR representation: transcribe the coefficient matrix A_i (for lag
i) into its place in the full VAR coefficient matrix, A */
static void add_Ai_to_VAR_A (gretl_matrix *Ai, GRETL_VAR *vecm, int k)
{
int i, j, offset = k * vecm->neqns;
double x;
for (i=0; i<vecm->neqns; i++) {
for (j=0; j<vecm->neqns; j++) {
x = gretl_matrix_get(Ai, i, j);
gretl_matrix_set(vecm->A, i, j + offset, x);
}
}
}
/* Run OLS taking the betas, as calculated via the eigen-analysis,
as given. So obtain estimates and standard errors for the
coefficients on the lagged differences and the unrestricted
deterministic vars. Construct full residuals matrix while
we're at it.
FIXME: when seasonals are included, we're not getting the
same results as JMulTi for the constant (though the results
are the same for the seasonal dummies themselves)?
*/
static int
build_VECM_models (GRETL_VAR *vecm, double ***pZ, DATAINFO *pdinfo, int iter)
{
gretl_matrix *Pi = NULL;
gretl_matrix *A = NULL;
gretl_matrix **G = NULL;
int rv0 = pdinfo->v;
int mt, t, r = vecm->jinfo->rank;
int p = vecm->order;
int nv = vecm->neqns;
int *biglist = vecm->jinfo->biglist;
gretlopt lsqopt = OPT_N | OPT_Z;
int i, j, k;
int err = 0;
/* Note: "vecm->order" is actually the order of the VAR system,
which corresponds to the number of lagged differences on the
RHS of the VAR system. We need that number of G matrices to
hold the coefficients on those lagged differences.
*/
#if JDEBUG
fprintf(stderr, "build_VECM_models: vecm->order = %d\n", vecm->order);
#endif
if (p < 0) {
return E_DATA;
}
/* for computing VAR representation */
Pi = gretl_matrix_alloc(nv, nv);
A = gretl_matrix_alloc(nv, nv);
if (Pi == NULL || A == NULL) {
err = E_ALLOC;
goto bailout;
}
if (p > 0) {
G = gretl_matrix_array_alloc_with_size(p, nv, nv);
if (G == NULL) {
err = E_ALLOC;
goto bailout;
}
}
if (vecm->jinfo->Alpha == NULL) {
vecm->jinfo->Alpha = gretl_matrix_alloc(nv, r);
if (vecm->jinfo->Alpha == NULL) {
err = E_ALLOC;
goto bailout;
}
}
if (iter > 0) {
/* bootstrapping: EC terms already in dataset */
rv0 -= r;
lsqopt |= OPT_A;
}
err = add_EC_terms_to_dataset(vecm, pZ, pdinfo, iter);
for (i=0; i<nv && !err; i++) {
biglist[1] = vecm->jinfo->difflist[i+1];
k = biglist[0] - r + 1;
for (j=0; j<r; j++) {
biglist[k++] = rv0 + j;
}
#if JDEBUG
printlist(biglist, "build_VECM_models: biglist");
#endif
*vecm->models[i] = lsq(biglist, pZ, pdinfo, OLS, lsqopt);
err = vecm->models[i]->errcode;
if (!err) {
vecm->models[i]->ID = i + 1;
vecm->models[i]->aux = AUX_VECM;
vecm->models[i]->adjrsq = NADBL;
if (p > 0) {
copy_coeffs_to_Gamma(vecm->models[i], i, G, p, nv);
}
copy_coeffs_to_Alpha(vecm, i, vecm->jinfo->Alpha, p);
for (t=0; t<vecm->T; t++) {
mt = t + vecm->t1;
gretl_matrix_set(vecm->E, t, i, vecm->models[i]->uhat[mt]);
}
if (i == 0) {
vecm->ncoeff = vecm->models[i]->ncoeff;
}
} else {
fprintf(stderr, "build_VECM_models: error %d from lsq, eqn %d, iter %d\n",
err, i + 1, iter);
}
}
if (!err) {
/* \Pi = \alpha \beta' */
err = form_Pi(vecm, vecm->jinfo->Alpha, Pi);
}
if (err) {
goto bailout;
}
#if JDEBUG
gretl_matrix_print(vecm->jinfo->Alpha, "Alpha from models");
gretl_matrix_print(Pi, "Pi");
for (i=0; i<p; i++) {
fprintf(stderr, "Gamma matrix, lag %d\n\n", i+1);
gretl_matrix_print(G[i], NULL);
}
#endif
if (p == 0) {
gretl_matrix_I(A, nv);
gretl_matrix_add_to(A, Pi);
add_Ai_to_VAR_A(A, vecm, 0);
} else {
for (i=0; i<=p; i++) {
if (i == 0) {
gretl_matrix_I(A, nv);
gretl_matrix_add_to(A, Pi);
gretl_matrix_add_to(A, G[0]);
} else if (i == p) {
gretl_matrix_zero(A);
gretl_matrix_subtract_from(A, G[i-1]);
} else {
gretl_matrix_copy_values(A, G[i]);
gretl_matrix_subtract_from(A, G[i-1]);
}
#if JDEBUG
fprintf(stderr, "A matrix, lag %d\n\n", i+1);
gretl_matrix_print(A, NULL);
#endif
add_Ai_to_VAR_A(A, vecm, i);
}
}
#if JDEBUG
gretl_matrix_print(vecm->A, "vecm->A");
#endif
bailout:
gretl_matrix_free(Pi);
gretl_matrix_array_free(G, p);
gretl_matrix_free(A);
return err;
}
/* Cointegration test: calculate and print the "long-run matrix",
\alpha \beta', (what Hamilton calls \Zeta_0) */
static int
compute_long_run_matrix (JohansenInfo *jv, int n, const DATAINFO *pdinfo, PRN *prn)
{
gretl_matrix *Z0 = NULL;
gretl_matrix *tmp = NULL;
int nv = gretl_matrix_rows(jv->Svv);
int err = 0;
Z0 = gretl_matrix_alloc(n, nv);
tmp = gretl_matrix_alloc(nv, nv);
if (Z0 == NULL || tmp == NULL) {
err = E_ALLOC;
}
if (!err) {
gretl_matrix_multiply_mod(jv->Beta, GRETL_MOD_NONE,
jv->Beta, GRETL_MOD_TRANSPOSE,
tmp, GRETL_MOD_NONE);
/* \Zeta_0 = \Sigma_{uv} A A' */
gretl_matrix_multiply(jv->Suv, tmp, Z0);
print_lr_matrix(jv, Z0, n, pdinfo, prn);
}
gretl_matrix_free(Z0);
gretl_matrix_free(tmp);
return err;
}
/* Print both "raw" and re-scaled versions of the beta and alpha
matrices (cointegrating vectors and vectors of adjustments
respectively).
*/
static int
print_beta_and_alpha (JohansenInfo *jv, gretl_matrix *evals, int h,
const DATAINFO *pdinfo, PRN *prn)
{
int i, err = 0;
pprintf(prn, "\n%s", _("eigenvalue"));
for (i=0; i<h; i++) {
pprintf(prn, "%#12.5g ", evals->val[i]);
}
pputc(prn, '\n');
/* "raw" vectors */
print_beta_or_alpha(jv, h, pdinfo, prn, PRINT_BETA, 0);
print_beta_or_alpha(jv, h, pdinfo, prn, PRINT_ALPHA, 0);
/* re-scaled versions */
print_beta_or_alpha(jv, h, pdinfo, prn, PRINT_BETA, 1);
print_beta_or_alpha(jv, h, pdinfo, prn, PRINT_ALPHA, 1);
pputc(prn, '\n');
return err;
}
/*
renormalize \beta such that its uppermost submatrix of
size rank * rank is the identity matrix:
\beta' = [ I | *free elements* ]
*/
static int phillips_normalize_beta (GRETL_VAR *vecm)
{
gretl_matrix *c = NULL;
gretl_matrix *beta_c = NULL;
int r = jrank(vecm);
int n = gretl_matrix_rows(vecm->jinfo->Beta);
int i, j, err = 0;
double x;
c = gretl_matrix_alloc(r, r);
beta_c = gretl_matrix_alloc(n, r);
if (c == NULL || beta_c == NULL) {
err = E_ALLOC;
goto bailout;
}
for (i=0; i<r; i++) {
for (j=0; j<r; j++) {
x = gretl_matrix_get(vecm->jinfo->Beta, i, j);
gretl_matrix_set(c, i, j, x);
}
}
/* form \beta_c = \beta c^{-1} */
gretl_invert_general_matrix(c);
gretl_matrix_multiply(vecm->jinfo->Beta, c, beta_c);
/* correct rounding error: set true zeros in \beta_c */
for (i=0; i<r; i++) {
for (j=0; j<r; j++) {
if (i == j) {
gretl_matrix_set(beta_c, i, j, 1.0);
} else {
gretl_matrix_set(beta_c, i, j, 0.0);
}
}
}
#if JDEBUG
gretl_matrix_print(vecm->jinfo->Beta, "original beta");
gretl_matrix_print(beta_c, "beta_c = beta * c^{-1}");
#endif
gretl_matrix_copy_values(vecm->jinfo->Beta, beta_c);
bailout:
gretl_matrix_free(c);
gretl_matrix_free(beta_c);
return err;
}
/* VECM: compute the variance of the estimator of \beta, after doing
Phillips normalization */
static int beta_variance (GRETL_VAR *vecm)
{
gretl_matrix *O = NULL;
gretl_matrix *aOa = NULL;
gretl_matrix *HSH = NULL;
int r = jrank(vecm);
int m = gretl_matrix_cols(vecm->jinfo->Alpha);
int n = gretl_matrix_rows(vecm->jinfo->Beta);
int i, j, k, err = 0;
double x;
O = gretl_matrix_copy(vecm->S);
aOa = gretl_matrix_alloc(m, m);
HSH = gretl_matrix_alloc(n - r, n - r);
if (O == NULL || aOa == NULL || HSH == NULL) {
err = E_ALLOC;
goto bailout;
}
/* compute \alpha' \Omega^{-1} \alpha */
err = gretl_invert_symmetric_matrix(O);
if (err) {
goto bailout;
}
gretl_matrix_qform(vecm->jinfo->Alpha, GRETL_MOD_TRANSPOSE, O,
aOa, GRETL_MOD_NONE);
#if JDEBUG
gretl_matrix_print(vecm->S, "vecm->S");
gretl_matrix_print(O, "O = inverse(vecm->S)");
gretl_matrix_print(vecm->jinfo->Alpha, "alpha_c");
gretl_matrix_print(aOa, "aOa = alpha_c' * O * alpha_c");
#endif
/* compute H'SH (just keep the south-east corner) */
for (i=r; i<n; i++) {
for (j=r; j<n; j++) {
x = gretl_matrix_get(vecm->jinfo->Svv, i, j);
gretl_matrix_set(HSH, i - r, j - r, x);
}
}
#if JDEBUG
gretl_matrix_print(vecm->jinfo->Svv, "full Svv");
gretl_matrix_print(HSH, "HSH = subset(Svv)");
#endif
vecm->jinfo->Bvar = gretl_matrix_kronecker_product_new(aOa, HSH);
if (vecm->jinfo->Bvar == NULL) {
err = E_ALLOC;
goto bailout;
}
err = gretl_invert_symmetric_matrix(vecm->jinfo->Bvar);
if (err) {
goto bailout;
}
gretl_matrix_divide_by_scalar(vecm->jinfo->Bvar, vecm->T);
vecm->jinfo->Bse = gretl_matrix_alloc(n - r, r);
if (vecm->jinfo->Bse == NULL) {
err = E_ALLOC;
goto bailout;
}
k = 0;
for (j=0; j<r; j++) {
/* cointegrating vector j */
for (i=0; i<n-r; i++) {
x = gretl_matrix_get(vecm->jinfo->Bvar, k, k);
gretl_matrix_set(vecm->jinfo->Bse, i, j, sqrt(x));
k++;
}
}
#if JDEBUG
gretl_matrix_print(vecm->jinfo->Bvar, "var(beta)");
gretl_matrix_print(vecm->jinfo->Bse, "se(beta)");
#endif
bailout:
gretl_matrix_free(O);
gretl_matrix_free(aOa);
gretl_matrix_free(HSH);
return err;
}
static int johansen_ll_calc (GRETL_VAR *jvar, const gretl_matrix *eigvals)
{
gretl_matrix *Suu;
double ldet, T_2 = (double) jvar->T / 2.0;
int n = jvar->neqns;
int h, i, err = 0;
h = (jrank(jvar) > 0)? jrank(jvar) : n;
Suu = gretl_matrix_copy(jvar->jinfo->Suu);
if (Suu == NULL) {
err = E_ALLOC;
} else {
ldet = gretl_matrix_log_determinant(Suu, &err);
jvar->ll = - T_2 * n * (1.0 + LN_2_PI) - T_2 * ldet;
for (i=0; i<h; i++) {
jvar->ll -= T_2 * log(1.0 - eigvals->val[i]);
}
gretl_matrix_free(Suu);
}
return err;
}
static int vecm_ll_stats (GRETL_VAR *vecm)
{
gretl_matrix *S;
int T = vecm->T;
int g = vecm->neqns;
int k = g * (vecm->order + 1);
S = gretl_matrix_copy(vecm->S);
if (S == NULL) {
return E_ALLOC;
}
vecm->ldet = gretl_vcv_log_determinant(S);
gretl_matrix_free(S);
/* FIXME: is k right (in all cases)? */
k += vecm->jinfo->nexo;
vecm->AIC = (-2.0 * vecm->ll + 2.0 * k * g) / T;
vecm->BIC = (-2.0 * vecm->ll + log(T) * k * g) / T;
vecm->HQC = (-2.0 * vecm->ll + 2.0 * log(log(T)) * k * g) / T;
return 0;
}
static int
compute_coint_test (GRETL_VAR *jvar, const gretl_matrix *evals, PRN *prn)
{
int T = jvar->T;
int n = jvar->neqns;
double cumeig = 0.0;
double *lmax = NULL;
double *trace = NULL;
double pvals[2];
int i;
trace = malloc(n * sizeof *trace);
lmax = malloc(n * sizeof *lmax);
if (trace == NULL || lmax == NULL) {
free(trace);
free(lmax);
return E_ALLOC;
}
for (i=n-1; i>=0; i--){
lmax[i] = -T * log(1.0 - evals->val[i]);
cumeig += lmax[i];
trace[i] = cumeig;
}
pputc(prn, '\n');
print_Johansen_test_case(jcode(jvar), prn);
pprintf(prn, "\n%s %s %s %s %s %s\n", _("Rank"), _("Eigenvalue"),
_("Trace test"), _("p-value"),
_("Lmax test"), _("p-value"));
for (i=0; i<n; i++) {
gamma_par_asymp(trace[i], lmax[i], jcode(jvar), n - i, pvals);
pprintf(prn, "%4d%#11.5g%#11.5g [%6.4f]%#11.5g [%6.4f]\n",
i, evals->val[i], trace[i], pvals[0], lmax[i], pvals[1]);
}
pputc(prn, '\n');
free(lmax);
free(trace);
return 0;
}
static int johansen_get_eigenvalues (gretl_matrix *Suu,
const gretl_matrix *Suv,
const gretl_matrix *Svv,
gretl_matrix *M,
gretl_matrix **evals,
int rank)
{
gretl_matrix *Tmp = NULL;
int n = Svv->cols;
int err;
err = gretl_invert_symmetric_matrix(Suu);
if (err) {
return err;
}
Tmp = gretl_matrix_alloc(n, n);
if (Tmp == NULL) {
return E_ALLOC;
}
gretl_matrix_qform(Suv, GRETL_MOD_TRANSPOSE,
Suu, Tmp, GRETL_MOD_NONE);
*evals = gretl_gensymm_eigenvals(Tmp, Svv, M, &err);
if (!err) {
err = gretl_symmetric_eigen_sort(*evals, M, rank);
}
gretl_matrix_free(Tmp);
return err;
}
/* Public entry point for both cointegration test and VECM
estimation */
int johansen_analysis (GRETL_VAR *jvar, double ***pZ, DATAINFO *pdinfo,
gretlopt opt, PRN *prn)
{
gretl_matrix *M = NULL;
gretl_matrix *Suu = NULL;
gretl_matrix *evals = NULL;
int n = jvar->neqns;
int nv = gretl_matrix_cols(jvar->jinfo->Svv);
int rank = jrank(jvar);
int err = 0;
#if JDEBUG
fprintf(stderr, "\n*** starting johansen_analysis()\n\n");
#endif
M = gretl_matrix_alloc(nv, nv);
/* We want to preserve Suu, so copy before processing */
Suu = gretl_matrix_copy(jvar->jinfo->Suu);
if (Suu == NULL || M == NULL) {
err = E_ALLOC;
goto bailout;
}
err = johansen_get_eigenvalues(Suu, jvar->jinfo->Suv, jvar->jinfo->Svv,
M, &evals, rank);
if (err) {
pputs(prn, _("Failed to find eigenvalues\n"));
goto bailout;
}
#if JDEBUG
gretl_matrix_print(M, "raw eigenvector(s)");
#endif
if (!err) {
johansen_ll_calc(jvar, evals);
if (rank == 0) {
compute_coint_test(jvar, evals, prn);
/* just running cointegration test */
if (opt & OPT_Q) {
goto bailout;
}
jvar->jinfo->Beta = M;
M = NULL;
err = compute_alpha(jvar->jinfo, n);
if (!err) {
print_beta_and_alpha(jvar->jinfo, evals, n, pdinfo, prn);
compute_long_run_matrix(jvar->jinfo, n, pdinfo, prn);
}
} else {
/* estimating VECM */
int do_stderrs = rank < jvar->neqns;
jvar->jinfo->Beta = gretl_matrix_copy(M);
if (jvar->jinfo->Beta == NULL) {
err = E_ALLOC;
}
if (!err) {
err = phillips_normalize_beta(jvar);
}
if (!err) {
err = build_VECM_models(jvar, pZ, pdinfo, 0);
}
if (!err) {
err = compute_omega(jvar);
}
if (!err && do_stderrs) {
err = beta_variance(jvar);
}
if (!err) {
err = gretl_VAR_do_error_decomp(jvar->S, jvar->C);
}
if (!err) {
err = vecm_ll_stats(jvar);
}
}
}
bailout:
gretl_matrix_free(M);
gretl_matrix_free(Suu);
gretl_matrix_free(evals);
return err;
}
/* Simplified version of the Johansen procedure, to be called in
the process of computing bootstrap confidence intervals for
impulse response functions. We just have to do enough to
generate the VAR representation.
*/
int
johansen_bootstrap_round (GRETL_VAR *jvar, double ***pZ, DATAINFO *pdinfo,
int iter)
{
gretl_matrix *M = NULL;
gretl_matrix *evals = NULL;
int nv = gretl_matrix_cols(jvar->jinfo->Svv);
int err = 0;
#if JDEBUG
fprintf(stderr, "\n*** starting johansen_bootstrap_round()\n\n");
#endif
M = gretl_matrix_alloc(nv, nv);
if (M == NULL) {
err = E_ALLOC;
goto bailout;
}
err = johansen_get_eigenvalues(jvar->jinfo->Suu, jvar->jinfo->Suv,
jvar->jinfo->Svv, M, &evals,
jrank(jvar));
if (err) {
goto bailout;
}
#if JDEBUG
gretl_matrix_print(M, "raw eigenvector(s)");
#endif
if (!err) {
if (jvar->jinfo->Beta == NULL) {
jvar->jinfo->Beta = gretl_matrix_copy(M);
} else {
gretl_matrix_copy_values(jvar->jinfo->Beta, M);
}
if (jvar->jinfo->Beta == NULL) {
err = E_ALLOC;
}
if (!err) {
err = phillips_normalize_beta(jvar);
}
if (!err) {
err = build_VECM_models(jvar, pZ, pdinfo, iter);
}
if (!err) {
err = compute_omega(jvar);
}
}
bailout:
gretl_matrix_free(M);
gretl_matrix_free(evals);
return err;
}
static int
johansen_LR_calc (GRETL_VAR *jvar, const gretl_matrix *evals, PRN *prn)
{
gretl_matrix *Suu;
double llr = 0.0;
double ldet = 0.0;
double T_2 = (double) jvar->T / 2.0;
int n = jvar->neqns;
int h, i, err = 0;
h = (jrank(jvar) > 0)? jrank(jvar) : n;
Suu = gretl_matrix_copy(jvar->jinfo->Suu);
if (Suu == NULL) {
err = E_ALLOC;
} else {
ldet = gretl_matrix_log_determinant(Suu, &err);
}
if (!err) {
llr = - T_2 * n * (1.0 + LN_2_PI) - T_2 * ldet;
for (i=0; i<h; i++) {
pprintf(prn, "eigenvalue %d = %g\n", i+1, evals->val[i]);
llr -= T_2 * log(1.0 - evals->val[i]);
}
pputc(prn, '\n');
}
if (Suu != NULL) {
gretl_matrix_free(Suu);
}
if (!err) {
double x = 2.0 * (jvar->ll - llr);
int nb = gretl_matrix_rows(jvar->jinfo->Beta);
int df = h * (nb - gretl_matrix_cols(jvar->jinfo->D));
pprintf(prn, "Unrestricted loglikelihood (lu) = %g\n", jvar->ll);
pprintf(prn, "Restricted loglikelihood (lr) = %g\n", llr);
pprintf(prn, "2 * (lu - lr) = %g\n", x);
pprintf(prn, "P(Chi-Square(%d) > %g = %g\n", df, x,
chisq_cdf_comp(x, df));
}
return err;
}
/*
Test of (homogeneous) linear restrictions on the cointegrating
relations in a VECM. This all needs verification and possibly
fixing in the case where the "unrestricted" VECM includes a
restricted constant or trend
*/
int vecm_beta_test (GRETL_VAR *jvar, const DATAINFO *pdinfo, PRN *prn)
{
gretl_matrix *M = NULL;
gretl_matrix *Svv = NULL;
gretl_matrix *Suv = NULL;
gretl_matrix *Suu = NULL;
gretl_matrix *evals = NULL;
int n = jvar->neqns;
int m = gretl_matrix_cols(jvar->jinfo->D);
int rank = jrank(jvar);
int err = 0;
M = gretl_matrix_alloc(m, m);
Svv = gretl_matrix_alloc(m, m);
Suv = gretl_matrix_alloc(n, m);
Suu = gretl_matrix_copy(jvar->jinfo->Suu);
if (M == NULL || Svv == NULL || Suv == NULL || Suu == NULL) {
err = E_ALLOC;
goto bailout;
}
pputs(prn, "\nTest of restrictions on cointegrating relations\n\n");
/* calculate Svv <- D' Svv D */
gretl_matrix_qform(jvar->jinfo->D, GRETL_MOD_TRANSPOSE,
jvar->jinfo->Svv, Svv, GRETL_MOD_NONE);
if (!err) {
/* Suv <- SuvD */
err = gretl_matrix_multiply(jvar->jinfo->Suv,
jvar->jinfo->D,
Suv);
}
err = johansen_get_eigenvalues(Suu, Suv, Svv, M, &evals, rank);
if (!err) {
johansen_LR_calc(jvar, evals, prn);
}
if (!err) {
gretl_matrix_multiply_mod(jvar->jinfo->D, GRETL_MOD_NONE,
M, GRETL_MOD_NONE,
jvar->jinfo->Beta, GRETL_MOD_NONE);
err = compute_alpha(jvar->jinfo, n);
}
if (!err) {
err = phillips_normalize_beta(jvar);
}
if (!err) {
print_beta_or_alpha(jvar->jinfo, rank, pdinfo, prn, PRINT_BETA, 0);
print_beta_or_alpha(jvar->jinfo, rank, pdinfo, prn, PRINT_ALPHA, 0);
pputc(prn, '\n');
compute_long_run_matrix(jvar->jinfo, n, pdinfo, prn);
}
bailout:
gretl_matrix_free(M);
gretl_matrix_free(Svv);
gretl_matrix_free(Suv);
gretl_matrix_free(Suu);
gretl_matrix_free(evals);
return err;
}
syntax highlighted by Code2HTML, v. 0.9.1