/* 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.
*/
#include "libgretl.h"
#define ADEBUG 0
typedef struct arbond_ arbond;
struct unit_info {
int t1; /* first usable obs for unit */
int t2; /* last usable obs */
char *skip; /* mask for obs to be skipped, if any */
};
struct diag_info {
int v; /* ID number of variable */
int minlag; /* minimum lag order */
int maxlag; /* maximum lag order */
};
struct arbond_ {
gretlopt opt; /* option flags */
int doZX; /* calculation flag */
int step; /* what step are we on? */
int yno; /* ID number of dependent var */
int p; /* lag order for dependent variable */
int qmax; /* longest lag of y used as instrument */
int qmin; /* shortest lag of y used as instrument */
int nx; /* number of independent variables */
int nz; /* number of regular instruments */
int nzb; /* number of block-diagonal instruments (other than y) */
int m; /* number of columns in instrument matrix, Z */
int pc0; /* column in Z where predet vars start */
int xc0; /* column in Z where exog vars start */
int N; /* total number of units in sample */
int effN; /* number of units with usable observations */
int T; /* total number of observations per unit */
int maxTi; /* maximum equations for any given unit */
int k; /* number of parameters estimated */
int nobs; /* total observations used */
int t1; /* initial offset into dataset */
int t1min; /* first usable observation, any unit */
int ndum; /* number of time dummies to use */
double SSR; /* sum of squared residuals */
double s2; /* residual variance */
double AR1; /* z statistic for AR(1) errors */
double AR2; /* z statistic for AR(2) errors */
double sargan; /* overidentification test statistic */
double wald; /* Wald test statistic */
int wdf; /* degrees of freedom for Wald test */
int *xlist; /* list of independent variables */
int *ilist; /* list of instruments */
gretl_matrix *beta; /* parameter estimates */
gretl_matrix *vbeta; /* parameter variance matrix */
gretl_matrix *uhat; /* residuals, differenced version */
gretl_matrix *H; /* step 1 error covariance matrix */
gretl_matrix *A;
gretl_matrix *Acpy;
gretl_matrix *V;
gretl_matrix *ZT; /* transpose of full instrument matrix */
gretl_matrix *Zi; /* per-unit instrument matrix */
gretl_matrix *dy; /* first difference of dependent var */
gretl_matrix *dX; /* lagged differences of y and indep vars */
gretl_matrix *tmp1;
gretl_matrix *kmtmp;
gretl_matrix *kktmp;
gretl_matrix *den;
gretl_matrix *L1;
gretl_matrix *XZA;
gretl_matrix *R1;
gretl_matrix *ZX;
struct diag_info *d; /* info on block-diagonal instruments */
struct unit_info *ui; /* info on panel units */
};
#define data_index(ab,i) (i * ab->T + ab->t1)
static void arbond_free (arbond *ab)
{
int i;
gretl_matrix_free(ab->beta);
gretl_matrix_free(ab->vbeta);
gretl_matrix_free(ab->uhat);
gretl_matrix_free(ab->H);
gretl_matrix_free(ab->A);
gretl_matrix_free(ab->Acpy);
gretl_matrix_free(ab->V);
gretl_matrix_free(ab->Zi);
gretl_matrix_free(ab->ZT);
gretl_matrix_free(ab->dy);
gretl_matrix_free(ab->dX);
gretl_matrix_free(ab->tmp1);
gretl_matrix_free(ab->kmtmp);
gretl_matrix_free(ab->kktmp);
gretl_matrix_free(ab->den);
gretl_matrix_free(ab->L1);
gretl_matrix_free(ab->XZA);
gretl_matrix_free(ab->R1);
gretl_matrix_free(ab->ZX);
free(ab->xlist);
free(ab->ilist);
free(ab->d);
for (i=0; i<ab->N; i++) {
free(ab->ui[i].skip);
}
free(ab->ui);
}
static int arbond_allocate (arbond *ab)
{
int T2 = ab->maxTi;
ab->beta = gretl_matrix_alloc(ab->k, 1);
ab->vbeta = gretl_matrix_alloc(ab->k, ab->k);
ab->uhat = gretl_matrix_alloc(ab->nobs, 1);
ab->ZT = gretl_matrix_alloc(ab->m, ab->nobs);
ab->H = gretl_matrix_alloc(T2, T2);
ab->A = gretl_matrix_alloc(ab->m, ab->m);
ab->Acpy = gretl_matrix_alloc(ab->m, ab->m);
ab->Zi = gretl_matrix_alloc(T2, ab->m);
ab->dy = gretl_column_vector_alloc(ab->nobs);
ab->dX = gretl_matrix_alloc(ab->nobs, ab->k);
ab->tmp1 = gretl_matrix_alloc(ab->m, ab->m);
ab->kmtmp = gretl_matrix_alloc(ab->k, ab->m);
ab->kktmp = gretl_matrix_alloc(ab->k, ab->k);
ab->den = gretl_matrix_alloc(ab->k, ab->k);
ab->L1 = gretl_matrix_alloc(1, ab->m);
ab->XZA = gretl_matrix_alloc(ab->k, ab->m);
ab->R1 = gretl_matrix_alloc(ab->m, 1);
ab->ZX = gretl_matrix_alloc(ab->m, ab->k);
if (ab->ZT == NULL || ab->H == NULL || ab->A == NULL ||
ab->Zi == NULL || ab->dy == NULL || ab->dX == NULL ||
ab->tmp1 == NULL || ab->kmtmp == NULL ||
ab->kktmp == NULL || ab->L1 == NULL || ab->XZA == NULL ||
ab->R1 == NULL || ab->ZX == NULL || ab->den == NULL ||
ab->beta == NULL || ab->vbeta == NULL ||
ab->uhat == NULL || ab->Acpy == NULL) {
return E_ALLOC;
} else {
return 0;
}
}
/* if the const has been included among the regressors but not the
instruments, add it to the instruments */
static int maybe_add_const_to_ilist (arbond *ab)
{
int i, addc = 0;
int err = 0;
if (ab->xlist == NULL || (ab->nz == 0 && ab->nzb == 0)) {
/* no x's, or all x's treated as exogenous already */
return 0;
}
for (i=1; i<=ab->xlist[0]; i++) {
if (ab->xlist[i] == 0) {
addc = 1;
break;
}
}
if (addc && ab->ilist != NULL) {
for (i=1; i<=ab->ilist[0]; i++) {
if (ab->ilist[i] == 0) {
addc = 0;
break;
}
}
}
if (addc) {
ab->ilist = gretl_list_append_term(&ab->ilist, 0);
if (ab->ilist == NULL) {
err = E_ALLOC;
} else {
ab->nz += 1;
}
}
return err;
}
static int arbond_make_lists (arbond *ab, const int *list, int xpos,
const DATAINFO *pdinfo)
{
int i, nz = 0, spos = 0;
int err = 0;
for (i=xpos; i<=list[0]; i++) {
if (list[i] == LISTSEP) {
spos = i;
break;
}
}
#if ADEBUG
printlist(list, "incoming list in arbond_make_lists");
fprintf(stderr, "separator pos = %d\n", spos);
#endif
if (spos > 0) {
ab->nx = spos - xpos;
nz = list[0] - spos;
} else {
ab->nx = list[0] - (xpos - 1);
}
#if ADEBUG
fprintf(stderr, "got intial ab->nx = %d, nz = %d\n", ab->nx, nz);
#endif
if (ab->nx > 0) {
/* compose indep vars list */
ab->xlist = gretl_list_new(ab->nx);
if (ab->xlist == NULL) {
return E_ALLOC;
}
for (i=0; i<ab->nx; i++) {
ab->xlist[i+1] = list[xpos + i];
}
#if ADEBUG
printlist(ab->xlist, "ab->xlist");
#endif
if (nz == 0 && ab->nzb == 0) {
/* implicitly all x vars are exogenous */
ab->ilist = gretl_list_copy(ab->xlist);
if (ab->ilist == NULL) {
return E_ALLOC;
}
ab->nz = ab->ilist[0];
}
}
if (nz > 0) {
/* compose regular instruments list */
ab->ilist = gretl_list_new(nz);
if (ab->ilist == NULL) {
return E_ALLOC;
}
for (i=0; i<nz; i++) {
ab->ilist[i+1] = list[spos + i + 1];
}
ab->nz = nz;
}
err = maybe_add_const_to_ilist(ab);
#if ADEBUG
printlist(ab->ilist, "ab->ilist");
#endif
return err;
}
static void arbond_zero_matrices (arbond *ab)
{
ab->beta = NULL;
ab->vbeta = NULL;
ab->uhat = NULL;
ab->ZT = NULL;
ab->H = NULL;
ab->A = NULL;
ab->Acpy = NULL;
ab->V = NULL;
ab->Zi = NULL;
ab->dy = NULL;
ab->dX = NULL;
ab->tmp1 = NULL;
ab->kmtmp = NULL;
ab->kktmp = NULL;
ab->den = NULL;
ab->L1 = NULL;
ab->XZA = NULL;
ab->R1 = NULL;
ab->ZX = NULL;
}
static void arbond_free_lists (arbond *ab)
{
free(ab->xlist);
free(ab->ilist);
}
static int
arbond_init (arbond *ab, const int *list, const DATAINFO *pdinfo,
gretlopt opt, struct diag_info *d, int nzb)
{
int i, xpos = 0;
int err = 0;
ab->d = d;
ab->nzb = nzb;
if (list[0] < 3) {
return E_PARSE;
}
ab->p = list[1];
ab->qmin = 2;
if (list[2] == LISTSEP) {
ab->qmax = 0;
ab->yno = list[3];
xpos = 4;
} else if (list[3] == LISTSEP && list[0] >= 4) {
ab->qmax = list[2];
ab->yno = list[4];
xpos = 5;
} else if (list[4] == LISTSEP && list[0] >= 5) {
ab->qmax = list[2];
ab->qmin = list[3];
ab->yno = list[5];
xpos = 6;
} else {
return E_PARSE;
}
if (ab->p < 1 || (ab->qmax != 0 && ab->qmax < ab->p + 1) ||
(ab->qmax != 0 && ab->qmin > ab->qmax)) {
/* is this all right? */
return E_INVARG;
}
/* FIXME: qmin doesn't actually do anything below, yet */
ab->xlist = NULL;
ab->ilist = NULL;
ab->opt = opt;
ab->step = 1;
ab->doZX = 1;
ab->nx = 0;
ab->nz = 0;
ab->t1min = 0;
ab->ndum = 0;
if (list[0] >= xpos) {
err = arbond_make_lists(ab, list, xpos, pdinfo);
if (err) {
return err;
}
}
ab->t1 = pdinfo->t1;
ab->T = pdinfo->pd;
ab->effN = ab->N = (pdinfo->t2 - ab->t1 + 1) / ab->T;
ab->k = ab->p + ab->nx;
ab->maxTi = 0;
ab->m = 0;
ab->pc0 = 0;
ab->xc0 = 0;
ab->nobs = 0;
ab->SSR = NADBL;
ab->s2 = NADBL;
ab->AR1 = NADBL;
ab->AR2 = NADBL;
ab->sargan = NADBL;
ab->wald = NADBL;
ab->wdf = 0;
#if ADEBUG
fprintf(stderr, "yno = %d, p = %d, qmax = %d, qmin = %d, nx = %d, k = %d\n",
ab->yno, ab->p, ab->qmax, ab->qmin, ab->nx, ab->k);
fprintf(stderr, "t1 = %d, T = %d, N = %d\n", ab->t1, ab->T, ab->N);
#endif
arbond_zero_matrices(ab);
ab->ui = malloc(ab->N * sizeof *ab->ui);
if (ab->ui == NULL) {
arbond_free_lists(ab);
err = E_ALLOC;
} else {
for (i=0; i<ab->N; i++) {
ab->ui[i].skip = NULL;
}
}
return err;
}
/* see if there are NAs for the dependent variable, or
for any of the independent variables, at a given
observation */
static int anymiss (arbond *ab, const double **Z, int s)
{
int i;
for (i=0; i<=ab->p+1; i++) {
if (na(Z[ab->yno][s-i])) {
return 1;
}
}
if (ab->xlist != NULL) {
for (i=1; i<=ab->xlist[0]; i++) {
if (na(Z[ab->xlist[i]][s])) {
return 1;
}
}
}
return 0;
}
static int bzcols (arbond *ab, int i)
{
int j, k, nc = 0;
int t = i + ab->p + 1; /* ?? */
for (j=0; j<ab->nzb; j++) {
for (k=ab->d[j].minlag; k<=ab->d[j].maxlag; k++) {
if (t - k >= 0) {
nc++;
}
}
}
return nc;
}
static int
arbond_sample_check (arbond *ab, const int *list,
const double **Z, const DATAINFO *pdinfo)
{
const double *y = Z[ab->yno];
char *mask = NULL;
int t1min = ab->T - 1;
int t1imin = ab->T - 1;
int t2max = 0;
int i, s, t;
int err = 0;
for (i=0; i<ab->N; i++) {
/* find the first y observation, all units */
if (t1min > 0) {
s = data_index(ab, i);
for (t=0; t<ab->T; t++) {
if (!na(y[s+t])) {
if (t < t1min) {
t1min = t;
}
break;
}
}
}
}
#if ADEBUG
fprintf(stderr, "arbond_sample_check, initial scan: ab->T = %d, t1min = %d\n",
ab->T, t1min);
#endif
mask = malloc(ab->T);
if (mask == NULL) {
return E_ALLOC;
}
if (ab->qmax == 0) {
ab->qmax = ab->T;
}
for (i=0; i<ab->N; i++) {
int t1i = ab->T - 1, t2i = 0;
int Ti = 0, usable = 0;
#if ADEBUG
fprintf(stderr, "Checking unit %d\n", i);
#endif
for (t=0; t<ab->T; t++) {
mask[t] = 0;
}
/* identify the observations at which we can form the required
Delta y terms, have the requisite independent variables,
and can construct at least one orthogonality condition
using a lagged level of y
*/
s = data_index(ab, i);
for (t=ab->p+1; t<ab->T; t++) {
if (anymiss(ab, Z, s + t)) {
mask[t] = 1;
} else {
usable++;
if (t < t1i) t1i = t;
if (t > t2i) t2i = t;
}
}
if (usable == 0) {
ab->effN -= 1;
ab->ui[i].t1 = -1;
ab->ui[i].t2 = -1;
#if ADEBUG
fprintf(stderr, "not usable\n");
#endif
continue;
}
Ti = t2i - t1i + 1;
if (usable < Ti) {
/* there were gaps: steal and replace the skip mask */
ab->ui[i].skip = mask;
mask = NULL;
if (i < ab->N - 1) {
mask = malloc(ab->T);
if (mask == NULL) {
return E_ALLOC;
}
}
}
if (usable > ab->maxTi) {
ab->maxTi = usable;
}
ab->nobs += usable;
if (t1i < t1imin) {
t1imin = t1i;
}
if (t2i > t2max) {
t2max = t2i;
}
ab->ui[i].t1 = t1i;
ab->ui[i].t2 = t2i;
#if ADEBUG
fprintf(stderr, "t1 = %d, t2 = %d, Ti = %d, usable obs = %d\n",
t1i, t2i, Ti, usable);
#endif
}
/* record first usable obs, any unit */
ab->t1min = t1imin;
/* figure number of time dummies, if wanted */
if (ab->opt & OPT_D) {
ab->ndum = t2max - t1imin;
ab->k += ab->ndum;
}
/* is t1min actually "reachable"? */
if (t1min < t1imin - ab->qmax) {
t1min = t1imin - ab->qmax;
}
#if ADEBUG
fprintf(stderr, "Number of units with usable observations = %d\n", ab->effN);
fprintf(stderr, "Total usable observations = %d\n", ab->nobs);
fprintf(stderr, "Max equations for any given unit = %d\n", ab->maxTi);
fprintf(stderr, "Maximal relevant time-series span: %d to %d = %d\n",
t1min, t2max, t2max - t1min + 1);
#endif
if (ab->effN == 0) {
err = E_MISSDATA;
} else {
/* compute the number of columns in Zi */
int tau = t2max - t1min + 1;
int nblocks = tau - ab->p - 1; /* was tau - 2 */
int bcols, cols = 0;
#if ADEBUG
fprintf(stderr, "\ntau = %d (ab->p = %d)\n", tau, ab->p);
#endif
for (i=0; i<nblocks; i++) {
/* lagged y values */
cols = (ab->p + i > ab->qmax - 1)? ab->qmax - 1 : ab->p + i;
ab->m += cols;
#if ADEBUG
fprintf(stderr, "i=%d: adding %d cols for y-lags\n", i, cols);
#endif
if (ab->nzb > 0) {
/* other block-diagonal instruments */
bcols = bzcols(ab, i);
ab->m += bcols;
#if ADEBUG
fprintf(stderr, " plus %d cols for z-lags\n", bcols);
#endif
}
}
#if ADEBUG
fprintf(stderr, "'basic' m = %d\n", ab->m);
#endif
ab->qmax = cols + 1;
/* record the column where the exogenous vars start */
ab->xc0 = ab->m;
ab->m += ab->nz;
ab->m += ab->ndum;
#if ADEBUG
fprintf(stderr, "total m = %d (dummies=%d, exog=%d)\n",
ab->m, ab->ndum, ab->nz);
#endif
}
free(mask);
return err;
}
static int unit_nobs (arbond *ab, int i)
{
int n;
if (ab->ui[i].t1 < 0) return 0;
n = ab->ui[i].t2 - ab->ui[i].t1 + 1;
if (ab->ui[i].skip != NULL) {
int t;
for (t=ab->ui[i].t1; t<=ab->ui[i].t2; t++) {
if (ab->ui[i].skip[t]) {
n--;
}
}
}
return n;
}
static int skip_obs (arbond *ab, int i, int t)
{
return (ab->ui[i].skip == NULL)? 0 : ab->ui[i].skip[t];
}
static int skip_unit (arbond *ab, int i)
{
return ab->ui[i].t1 < 0;
}
static int ar_skip_unit (arbond *ab, int i, int k)
{
int t;
for (t=ab->ui[i].t1+k; t<=ab->ui[i].t2; t++) {
if (!skip_obs(ab, i, t) && !skip_obs(ab, i, t-k)) {
return 0;
}
}
return 1;
}
/* see if we have sufficient data to calculate A & B's z statistic for
AR(k) errors: return the length of the needed arrays, or 0 if we
can't do it
*/
static int ar_data_check (arbond *ab, int k)
{
int i, t, T = 0;
for (i=0; i<ab->N; i++) {
if (skip_unit(ab, i)) {
continue;
}
for (t=ab->ui[i].t1+k; t<=ab->ui[i].t2; t++) {
if (!skip_obs(ab, i, t) && !skip_obs(ab, i, t-k)) {
T++;
}
}
}
return T;
}
static int arbond_const_pos (arbond *ab)
{
int i;
if (ab->xlist == NULL) {
return 0;
}
for (i=1; i<=ab->xlist[0]; i++) {
if (ab->xlist[i] == 0) {
return i;
}
}
return 0;
}
/* FIXME try to detect and omit periodic dummies? */
static int arbond_wald_test (arbond *ab)
{
gretl_matrix *vcv = NULL;
gretl_vector *b = NULL;
double x;
int cpos = arbond_const_pos(ab);
int i, j, k, kc = ab->p + ab->nx;
int ri, rj;
int err;
/* position of const in coeff vector? */
if (cpos != 0) {
k = kc - 1;
cpos += ab->p - 1;
} else {
k = kc;
cpos = -1;
}
b = gretl_column_vector_alloc(k);
vcv = gretl_matrix_alloc(k, k);
if (b == NULL || vcv == NULL) {
err = E_ALLOC;
goto bailout;
}
ri = 0;
for (i=0; i<kc; i++) {
if (i != cpos) {
b->val[ri++] = ab->beta->val[i];
}
}
ri = 0;
for (i=0; i<kc; i++) {
if (i != cpos) {
rj = 0;
for (j=0; j<kc; j++) {
if (j != cpos) {
x = gretl_matrix_get(ab->vbeta, i, j);
gretl_matrix_set(vcv, ri, rj++, x);
}
}
ri++;
}
}
err = gretl_invert_symmetric_matrix(vcv);
if (err) {
fprintf(stderr, "arbond_wald_test, error inverting vcv\n");
goto bailout;
}
x = gretl_scalar_qform(b, vcv, &err);
if (err) {
fprintf(stderr, _("Failed to compute test statistic\n"));
goto bailout;
}
if (!err) {
ab->wald = x;
ab->wdf = k;
}
#if ADEBUG
fprintf(stderr, "Wald chi^2(%d) = %g\n", k, x);
#endif
bailout:
gretl_matrix_free(vcv);
gretl_vector_free(b);
return err;
}
static int sargan_test (arbond *ab)
{
gretl_matrix *Zu = NULL;
gretl_matrix *m1 = NULL;
int err = 0;
Zu = gretl_matrix_alloc(ab->m, 1);
m1 = gretl_matrix_alloc(ab->m, 1);
if (Zu == NULL || m1 == NULL) {
gretl_matrix_free(Zu);
gretl_matrix_free(m1);
return E_ALLOC;
}
gretl_matrix_multiply(ab->ZT, ab->uhat, Zu);
gretl_matrix_divide_by_scalar(ab->A, ab->effN);
gretl_matrix_multiply(ab->A, Zu, m1);
ab->sargan = gretl_matrix_dot_product(Zu, GRETL_MOD_TRANSPOSE,
m1, GRETL_MOD_NONE,
&err);
if (ab->step == 1) {
ab->sargan *= 2.0 / ab->s2; /* allow for scale factor in H */
}
#if ADEBUG
fprintf(stderr, "Sargan df = m - k = %d - %d\n",
ab->m, ab->k);
fprintf(stderr, "Sargan test: Chi-square(%d) = %g\n",
ab->m - ab->k, ab->sargan);
#endif
gretl_matrix_free(Zu);
gretl_matrix_free(m1);
return err;
}
/* Compute the z test for AR errors, if possible. This should
perhaps be rolled into arbond_variance() for the sake of
efficiency
*/
static int ar_test (arbond *ab, const gretl_matrix *C)
{
gretl_matrix *v = NULL;
gretl_matrix *vk = NULL;
gretl_matrix *X = NULL; /* this is the trimmed "X_{*}" */
gretl_matrix *vkX = NULL;
gretl_matrix *tmpk = NULL;
gretl_matrix *ui = NULL;
gretl_matrix *m1 = NULL;
gretl_matrix *SZv = NULL;
double x, num, den;
double den2, den3;
int s, t, q, Q;
int sbig, k = 1;
int i, j, err = 0;
restart:
Q = ar_data_check(ab, k);
if (Q == 0) {
if (k == 1) {
return 0;
} else {
/* don't leak */
goto bailout;
}
}
#if ADEBUG
fprintf(stderr, "AR(%d) test: number of usable obs = %d\n", k, Q);
#endif
if (k == 1) {
v = gretl_column_vector_alloc(Q);
vk = gretl_column_vector_alloc(Q);
X = gretl_matrix_alloc(Q, ab->k);
vkX = gretl_matrix_alloc(1, ab->k);
tmpk = gretl_matrix_alloc(1, ab->k);
ui = gretl_column_vector_alloc(ab->maxTi);
m1 = gretl_matrix_alloc(ab->m, 1);
SZv = gretl_zero_matrix_new(ab->m, 1);
if (v == NULL || vk == NULL || X == NULL ||
vkX == NULL || tmpk == NULL || ui == NULL ||
m1 == NULL || SZv == NULL) {
err = E_ALLOC;
goto bailout;
}
} else {
gretl_matrix_reuse(v, Q, 1);
gretl_matrix_reuse(vk, Q, 1);
gretl_matrix_reuse(X, Q, -1);
gretl_matrix_reuse(m1, ab->m, 1);
gretl_matrix_zero(SZv);
}
den = 0.0;
q = s = sbig = 0;
for (i=0; i<ab->N; i++) {
int Ti = unit_nobs(ab, i);
double den1i = 0.0;
int si = 0;
if (Ti == 0) {
continue;
}
if (ar_skip_unit(ab, i, k)) {
sbig += Ti;
s += Ti;
continue;
}
gretl_matrix_reuse(ui, Ti, -1);
gretl_matrix_reuse(ab->Zi, ab->m, Ti);
/* extract full-length Z'_i and u_i */
for (t=ab->ui[i].t1; t<=ab->ui[i].t2; t++) {
if (!skip_obs(ab, i, t)) {
for (j=0; j<ab->m; j++) {
x = gretl_matrix_get(ab->ZT, j, sbig);
gretl_matrix_set(ab->Zi, j, si, x);
}
x = ab->uhat->val[sbig];
gretl_vector_set(ui, si, x);
sbig++;
si++;
}
}
for (t=ab->ui[i].t1; t<ab->ui[i].t1+k; t++) {
/* skip any obs prior to t1 + k */
if (!skip_obs(ab, i, t)) s++;
}
/* extract lagged residuals vk along with v_{*} and X_{*} */
for (t=ab->ui[i].t1+k; t<=ab->ui[i].t2; t++) {
if (!skip_obs(ab, i, t)) {
if (!skip_obs(ab, i, t-k)) {
v->val[q] = ab->uhat->val[s];
vk->val[q] = ab->uhat->val[s-k];
den1i += v->val[q] * vk->val[q];
for (j=0; j<ab->k; j++) {
x = gretl_matrix_get(ab->dX, s, j);
gretl_matrix_set(X, q, j, x);
}
q++;
}
s++;
}
}
/* cumulate Z'_i u_i u_i'_* u_{i,-k} */
gretl_matrix_multiply_by_scalar(ui, den1i);
gretl_matrix_multiply_mod(ab->Zi, GRETL_MOD_NONE,
ui, GRETL_MOD_NONE,
SZv, GRETL_MOD_CUMULATE);
den += den1i * den1i;
}
#if ADEBUG > 1
gretl_matrix_print(vk, "lagged residuals");
gretl_matrix_print(v, "current residuals");
#endif
/* the numerator */
num = 0.0;
for (i=0; i<Q; i++) {
num += vk->val[i] * v->val[i];
}
/* the denominator has three components, the first of which
is already handled */
/* required component "vkX" = \hat{v}'_{-k} X_{*} */
gretl_matrix_multiply_mod(vk, GRETL_MOD_TRANSPOSE,
X, GRETL_MOD_NONE,
vkX, GRETL_MOD_NONE);
/* vk' X_* (X'ZAZ'X)^{-1} X'ZA(sum stuff) */
gretl_matrix_multiply(vkX, C, tmpk);
gretl_matrix_reuse(m1, 1, ab->m);
gretl_matrix_multiply(tmpk, ab->XZA, m1);
den2 = gretl_matrix_dot_product(m1, GRETL_MOD_NONE,
SZv, GRETL_MOD_NONE,
&err);
den -= 2.0 * den2;
/* additive term vk' X_* vbeta X_*' vk */
gretl_matrix_multiply(vkX, ab->vbeta, tmpk);
den3 = gretl_matrix_dot_product(tmpk, GRETL_MOD_NONE,
vkX, GRETL_MOD_TRANSPOSE,
&err);
den += den3;
if (den < 0) {
err = E_NAN;
goto bailout;
}
/* now "x" = m1 or m2 */
x = num / sqrt(den);
#if ADEBUG
fprintf(stderr, "AR(%d) test: z = %.4g [%.3f]\n", k, x,
normal_pvalue_2(x));
#endif
if (k == 1) {
ab->AR1 = x;
k = 2;
goto restart;
} else {
ab->AR2 = x;
}
bailout:
gretl_matrix_free(v);
gretl_matrix_free(vk);
gretl_matrix_free(X);
gretl_matrix_free(vkX);
gretl_matrix_free(tmpk);
gretl_matrix_free(ui);
gretl_matrix_free(m1);
gretl_matrix_free(SZv);
return err;
}
/* Windmeijer, Journal of Econometrics, 126 (2005), page 33:
finite-sample correction for step-2 variance matrix. Note: from a
computational point of view this calculation is expressed more
clearly in the working paper by Bond and Windmeijer, "Finite Sample
Inference for GMM Estimators in Linear Panel Data Models" (2002),
in particular the fact that the matrix named "D" below must be
built column by column, taking the derivative of the "W" (or "A")
matrix with respect to the successive independent variables.
*/
static int windmeijer_correct (arbond *ab, const gretl_matrix *uhat1,
const gretl_matrix *varb1)
{
gretl_matrix *aV = NULL; /* standard asymptotic variance */
gretl_matrix *D = NULL; /* finite-sample factor */
gretl_matrix *dWj = NULL; /* one component of the above */
gretl_matrix *ui = NULL; /* per-unit residuals */
gretl_matrix *xij = NULL; /* per-unit X_j values */
gretl_matrix *TT = NULL; /* workspace follows */
gretl_matrix *mT = NULL;
gretl_matrix *km = NULL;
gretl_matrix *k1 = NULL;
int i, j, t;
int err = 0;
aV = gretl_matrix_copy(ab->vbeta);
D = gretl_matrix_alloc(ab->k, ab->k);
dWj = gretl_matrix_alloc(ab->m, ab->m);
ui = gretl_column_vector_alloc(ab->maxTi);
xij = gretl_column_vector_alloc(ab->maxTi);
TT = gretl_matrix_alloc(ab->maxTi, ab->maxTi);
mT = gretl_matrix_alloc(ab->m, ab->nobs);
km = gretl_matrix_alloc(ab->k, ab->m);
k1 = gretl_matrix_alloc(ab->k, 1);
if (aV == NULL || D == NULL || dWj == NULL ||
TT == NULL || mT == NULL || km == NULL ||
k1 == NULL || ui == NULL || xij == NULL) {
err = E_ALLOC;
goto bailout;
}
/* form -(1/N) * asyV * XZW^{-1} */
gretl_matrix_multiply(aV, ab->XZA, ab->kmtmp);
gretl_matrix_multiply_by_scalar(ab->kmtmp, -1.0 / ab->effN);
/* form W^{-1}Z'v_2 */
gretl_matrix_multiply(ab->A, ab->ZT, mT);
gretl_matrix_multiply(mT, ab->uhat, ab->R1);
for (j=0; j<ab->k; j++) { /* loop across the X's */
int s = 0;
gretl_matrix_zero(dWj);
/* form dWj = -(1/N) \sum Z_i'(X_j u' + u X_j')Z_i */
for (i=0; i<ab->N; i++) {
int Ti = unit_nobs(ab, i);
if (Ti == 0) {
continue;
}
gretl_matrix_reuse(ui, Ti, 1);
gretl_matrix_reuse(xij, Ti, 1);
gretl_matrix_reuse(ab->Zi, Ti, ab->m);
gretl_matrix_reuse(TT, Ti, Ti);
/* extract ui */
for (t=0; t<Ti; t++) {
ui->val[t] = uhat1->val[s++];
}
/* extract xij */
gretl_matrix_extract_matrix(xij, ab->dX, s - Ti, j,
GRETL_MOD_NONE);
gretl_matrix_multiply_mod(ui, GRETL_MOD_NONE,
xij, GRETL_MOD_TRANSPOSE,
TT, GRETL_MOD_NONE);
gretl_matrix_add_self_transpose(TT);
/* extract Zi */
gretl_matrix_extract_matrix(ab->Zi, ab->ZT, 0, s - Ti,
GRETL_MOD_TRANSPOSE);
gretl_matrix_qform(ab->Zi, GRETL_MOD_TRANSPOSE,
TT, dWj, GRETL_MOD_CUMULATE);
}
gretl_matrix_multiply_by_scalar(dWj, -1.0 / ab->effN);
/* D[.,j] = -aV * XZW^{-1} * dWj * W^{-1}Z'v_2 */
gretl_matrix_multiply(ab->kmtmp, dWj, km);
gretl_matrix_multiply(km, ab->R1, k1);
/* write into appropriate column of D (k x k) */
for (i=0; i<ab->k; i++) {
gretl_matrix_set(D, i, j, k1->val[i]);
}
}
/* add to AsyV: D * AsyV */
gretl_matrix_multiply_mod(D, GRETL_MOD_NONE,
aV, GRETL_MOD_NONE,
ab->vbeta, GRETL_MOD_CUMULATE);
/* add to AsyV: AsyV * D' */
gretl_matrix_multiply_mod(aV, GRETL_MOD_NONE,
D, GRETL_MOD_TRANSPOSE,
ab->vbeta, GRETL_MOD_CUMULATE);
/* add to AsyV: D * var(\hat{\beta}_1) * D' */
gretl_matrix_qform(D, GRETL_MOD_NONE, varb1,
ab->vbeta, GRETL_MOD_CUMULATE);
bailout:
gretl_matrix_free(D);
gretl_matrix_free(dWj);
gretl_matrix_free(ui);
gretl_matrix_free(xij);
gretl_matrix_free(TT);
gretl_matrix_free(mT);
gretl_matrix_free(km);
gretl_matrix_free(k1);
gretl_matrix_free(aV);
return err;
}
/*
Compute the variance matrix:
N * C^{-1} * (X'*Z*A_N*\hat{V}_N*A_N*Z'*X) * C^{-1}
where C = X'*Z*A_N*Z'*X and
\hat{V}_N = N^{-1} \sum Z_i'*v_i*v_i'*Z_i,
(v_i being the step-1 residuals)
we compute the residuals while we're at it
*/
static int arbond_variance (arbond *ab)
{
gretl_matrix *kk = NULL;
gretl_matrix *V = NULL;
gretl_matrix *ui = NULL;
gretl_matrix *u1 = NULL;
gretl_matrix *C = ab->den;
double x, ut, SSR;
int i, j, t, k, c;
int err = 0;
kk = gretl_matrix_alloc(ab->k, ab->k);
if (kk == NULL) {
return E_ALLOC;
}
if (ab->step == 1) {
V = gretl_zero_matrix_new(ab->m, ab->m);
ui = gretl_column_vector_alloc(ab->maxTi);
if (V == NULL || ui == NULL) {
err = E_ALLOC;
goto bailout;
}
} else {
u1 = gretl_matrix_copy(ab->uhat);
if (u1 == NULL) {
err = E_ALLOC;
goto bailout;
}
/* copy 1-step coefficient variance */
gretl_matrix_copy_values(kk, ab->vbeta);
}
SSR = 0.0;
c = k = 0;
for (i=0; i<ab->N; i++) {
int Ti = unit_nobs(ab, i);
if (Ti == 0) {
continue;
}
if (ab->step == 1) {
gretl_matrix_reuse(ab->Zi, Ti, ab->m);
gretl_matrix_reuse(ui, Ti, 1);
gretl_matrix_extract_matrix(ab->Zi, ab->ZT, 0, c,
GRETL_MOD_TRANSPOSE);
c += Ti;
}
for (t=0; t<Ti; t++) {
ut = ab->dy->val[k];
for (j=0; j<ab->k; j++) {
x = gretl_matrix_get(ab->dX, k, j);
ut -= ab->beta->val[j] * x;
}
SSR += ut * ut;
ab->uhat->val[k] = ut;
if (ui != NULL) {
ui->val[t] = ut;
}
k++;
}
if (ab->step == 1) {
gretl_matrix_multiply_mod(ui, GRETL_MOD_TRANSPOSE,
ab->Zi, GRETL_MOD_NONE,
ab->L1, GRETL_MOD_NONE);
gretl_matrix_multiply_mod(ab->L1, GRETL_MOD_TRANSPOSE,
ab->L1, GRETL_MOD_NONE,
V, GRETL_MOD_CUMULATE);
}
}
if (ab->step == 1) {
gretl_matrix_divide_by_scalar(V, ab->effN);
/* form X'Z A_N Z'X */
gretl_matrix_multiply_mod(ab->ZX, GRETL_MOD_TRANSPOSE,
ab->A, GRETL_MOD_NONE,
ab->kmtmp, GRETL_MOD_NONE);
gretl_matrix_qform(ab->kmtmp, GRETL_MOD_NONE, V,
kk, GRETL_MOD_NONE);
/* pre- and post-multiply by C^{-1} */
err = gretl_invert_symmetric_matrix(C);
gretl_matrix_qform(C, GRETL_MOD_NONE, kk, ab->vbeta,
GRETL_MOD_NONE);
if (ab->opt & OPT_T) {
/* preserve V for second stage */
ab->V = V;
V = NULL;
}
} else {
/* second step */
gretl_matrix_qform(ab->ZX, GRETL_MOD_TRANSPOSE, ab->V,
ab->vbeta, GRETL_MOD_NONE);
err = gretl_invert_symmetric_matrix(ab->vbeta);
if (!err) {
err = gretl_invert_symmetric_matrix(C); /* for AR test */
}
}
if (err) {
goto bailout;
}
gretl_matrix_multiply_by_scalar(ab->vbeta, ab->effN);
ab->SSR = SSR;
ab->s2 = SSR / (ab->nobs - ab->k);
if (ab->step == 2 && !(ab->opt & OPT_A)) {
windmeijer_correct(ab, u1, kk);
}
#if ADEBUG
gretl_matrix_print(ab->vbeta, "Var(beta)");
for (i=0; i<ab->k; i++) {
x = gretl_matrix_get(ab->vbeta, i, i);
fprintf(stderr, "se(beta[%d]) = %g\n", i, sqrt(x));
}
fprintf(stderr, "\nSSR = %.11g\n", ab->SSR);
fprintf(stderr, "sigma^2 = %.7g\n", ab->s2);
fprintf(stderr, "sigma = %.7g\n", sqrt(ab->s2));
#endif
/* while we're at it... */
ar_test(ab, C);
sargan_test(ab);
arbond_wald_test(ab);
bailout:
gretl_matrix_free(V);
gretl_matrix_free(ui);
gretl_matrix_free(u1);
gretl_matrix_free(kk);
return err;
}
static int next_obs (arbond *ab, int i, int j0, int n)
{
int j;
for (j=j0; j<n; j++) {
if (ab->ui[i].skip[j+ab->ui[i].t1] == 0) {
return j;
}
}
return 0;
}
/* construct the H matrix for first-differencing
as applied to unit i */
static int make_first_diff_matrix (arbond *ab, int i)
{
static int *rc;
int n, m;
double x;
int k, j, adjacent, skip = 0;
if (ab == NULL) {
/* clean-up signal */
free(rc);
rc = NULL;
return 0;
}
if (rc == NULL) {
rc = malloc((ab->T) * sizeof *rc);
if (rc == NULL) {
return E_ALLOC;
}
}
n = ab->ui[i].t2 - ab->ui[i].t1 + 1;
m = unit_nobs(ab, i);
if (m < n) {
skip = 1;
j = next_obs(ab, i, 0, n);
for (k=0; k<m; k++) {
rc[k] = j;
j = next_obs(ab, i, j+1, n);
}
}
gretl_matrix_reuse(ab->H, m, m);
for (j=0; j<m; j++) {
for (k=j; k<m; k++) {
if (skip) {
adjacent = (abs(rc[k] - rc[j]) == 1);
} else {
adjacent = (abs(k-j) == 1);
}
x = (k==j)? 2 : (adjacent)? -1 : 0;
gretl_matrix_set(ab->H, j, k, x);
gretl_matrix_set(ab->H, k, j, x);
}
}
return 0;
}
static int arbond_prepare_model (MODEL *pmod, arbond *ab,
const int *list, const char *istr,
const double **X, const DATAINFO *pdinfo)
{
const double *y = X[ab->yno];
double x;
int i, j;
int err = 0;
pmod->t1 = pdinfo->t1;
pmod->t2 = pdinfo->t2;
pmod->dfn = ab->k;
pmod->dfd = ab->nobs - ab->k;
pmod->list = gretl_list_copy(list);
if (pmod->list == NULL) {
pmod->errcode = E_ALLOC;
return pmod->errcode;
}
gretl_model_set_int(pmod, "yno", ab->yno);
gretl_model_set_int(pmod, "n_included_units", ab->effN);
pmod->ci = ARBOND;
pmod->ncoeff = ab->k;
pmod->nobs = ab->nobs;
pmod->full_n = pdinfo->n;
pmod->ess = ab->SSR;
if (ab->s2 >= 0) {
pmod->sigma = sqrt(ab->s2);
}
pmod->rsq = pmod->adjrsq = NADBL;
pmod->fstt = NADBL;
pmod->lnL = NADBL;
gretl_model_allocate_params(pmod, ab->k);
if (pmod->errcode) {
return pmod->errcode;
}
j = 0;
for (i=0; i<ab->p; i++) {
/* FIXME possible varname overflow */
sprintf(pmod->params[j++], "D%.10s(-%d)", pdinfo->varname[ab->yno], i+1);
}
for (i=0; i<ab->nx; i++) {
strcpy(pmod->params[j++], pdinfo->varname[ab->xlist[i+1]]);
}
for (i=0; i<ab->ndum; i++) {
sprintf(pmod->params[j++], "T%d", i + 2);
}
err = gretl_model_allocate_storage(pmod);
if (!err) {
gretl_model_new_vcv(pmod, NULL);
}
if (!err) {
for (i=0; i<ab->k; i++) {
pmod->coeff[i] = ab->beta->val[i];
for (j=0; j<=i; j++) {
x = gretl_matrix_get(ab->vbeta, i, j);
pmod->vcv[ijton(i, j, ab->k)] = x;
if (i == j) {
pmod->sderr[i] = sqrt(x);
}
}
}
}
if (!err) {
int s, t, k = 0;
/* add uhat, yhat */
for (i=0; i<ab->N; i++) {
if (skip_unit(ab, i)) {
continue;
}
for (t=0; t<ab->T; t++) {
if (t >= ab->ui[i].t1 && t <= ab->ui[i].t2) {
if (!skip_obs(ab, i, t)) {
s = data_index(ab, i) + t;
pmod->uhat[s] = ab->uhat->val[k];
pmod->yhat[s] = y[s] - pmod->uhat[s];
k++;
}
}
}
}
}
/* additional arbond-specific data */
if (!err) {
gretl_model_set_int(pmod, "step", ab->step);
if (!na(ab->AR1)) {
gretl_model_set_double(pmod, "AR1", ab->AR1);
}
if (!na(ab->AR2)) {
gretl_model_set_double(pmod, "AR2", ab->AR2);
}
if (!na(ab->sargan)) {
gretl_model_set_int(pmod, "sargan_df", ab->m - ab->k);
gretl_model_set_double(pmod, "sargan", ab->sargan);
}
if (!na(ab->wald)) {
gretl_model_set_int(pmod, "wald_df", ab->wdf);
gretl_model_set_double(pmod, "wald", ab->wald);
}
if (istr != NULL && *istr != '\0') {
gretl_model_set_string_as_data(pmod, "istr", gretl_strdup(istr));
}
if (ab->opt & OPT_D) {
gretl_model_set_int(pmod, "time-dummies", 1);
}
if ((ab->opt & OPT_T) && (ab->opt & OPT_A)) {
gretl_model_set_int(pmod, "asy", 1);
}
}
return err;
}
/* exclude any independent variables that are zero at all
relevant observations */
static int arbond_zero_check (arbond *ab, const double **Z)
{
const double *x;
int drop = 0;
int i, j, k, t;
for (j=0; j<ab->nx; j++) {
int all0 = 1;
x = Z[ab->xlist[j+1]];
for (i=0; i<ab->N && all0; i++) {
if (skip_unit(ab, i)) {
continue;
}
for (t=ab->ui[i].t1; t<=ab->ui[i].t2 && all0; t++) {
k = data_index(ab, i) + t;
if (x[k] != 0.0 && !na(x[k])) {
all0 = 0;
}
}
}
if (all0) {
ab->xlist[j+1] = -1;
drop++;
}
}
if (drop > 0) {
for (i=1; i<=ab->xlist[0]; i++) {
if (ab->xlist[i] < 0) {
gretl_list_delete_at_pos(ab->xlist, i);
i--;
}
}
ab->nx = ab->xlist[0];
ab->k -= drop;
}
return 0;
}
static int count_mask (const char *mask, int n)
{
int i, c = 0;
for (i=0; i<n; i++) {
c += mask[i] == 0;
}
return c;
}
static gretl_matrix *
matrix_copy_masked (const gretl_matrix *m, const char *mask)
{
int n = count_mask(mask, m->rows);
gretl_matrix *a;
double x;
int i, j, k, l;
a = gretl_matrix_alloc(n, n);
if (a == NULL) {
return NULL;
}
k = 0;
for (i=0; i<m->rows; i++) {
if (!mask[i]) {
l = 0;
for (j=0; j<m->cols; j++) {
if (!mask[j]) {
x = gretl_matrix_get(m, i, j);
gretl_matrix_set(a, k, l++, x);
}
}
k++;
}
}
return a;
}
static char *zero_row_mask (const gretl_matrix *m, int *err)
{
char *mask = NULL;
int any0 = 0, row0, i, j;
mask = calloc(m->rows, 1);
if (mask == NULL) {
*err = E_ALLOC;
return NULL;
}
for (i=0; i<m->rows; i++) {
row0 = 1;
for (j=0; j<m->cols; j++) {
if (gretl_matrix_get(m, i, j) != 0.0) {
row0 = 0;
break;
}
}
if (row0) {
mask[i] = 1;
any0 = 1;
}
}
if (!any0) {
free(mask);
mask = NULL;
}
return mask;
}
static void cut_rows (gretl_matrix *a, const char *mask)
{
int n = count_mask(mask, a->rows);
int i, j, k;
double x;
for (j=0; j<a->cols; j++) {
k = 0;
for (i=0; i<a->rows; i++) {
if (!mask[i]) {
x = gretl_matrix_get(a, i, j);
a->val[j * n + k] = x;
k++;
}
}
}
a->rows = n;
}
static char *
make_rank_mask (const gretl_matrix *A, int *err)
{
gretl_matrix *Q = NULL;
gretl_matrix *R = NULL;
char *mask = NULL;
double test;
int i, n = A->cols;
Q = gretl_matrix_copy(A);
if (Q == NULL) {
*err = E_ALLOC;
return NULL;
}
R = gretl_matrix_alloc(n, n);
if (R == NULL) {
gretl_matrix_free(Q);
*err = E_ALLOC;
return NULL;
}
*err = gretl_matrix_QR_decomp(Q, R);
if (!*err) {
mask = calloc(n, 1);
if (mask == NULL) {
*err = E_ALLOC;
}
}
if (!*err) {
for (i=0; i<n; i++) {
test = gretl_matrix_get(R, i, i);
if (fabs(test) < R_DIAG_MIN) {
mask[i] = 1;
}
}
}
if (*err) {
free(mask);
mask = NULL;
}
gretl_matrix_free(Q);
gretl_matrix_free(R);
return mask;
}
/* Based on reduction of the A matrix, trim ZT to match and
adjust the sizes of all workspace matrices that have a
dimension involving ab->m.
*/
static void real_shrink_matrices (arbond *ab, const char *mask)
{
fprintf(stderr, "A matrix: shrinking m from %d to %d\n",
ab->m, ab->A->rows);
cut_rows(ab->ZT, mask);
ab->m = ab->A->rows;
gretl_matrix_reuse(ab->Acpy, ab->m, ab->m);
gretl_matrix_reuse(ab->tmp1, ab->m, ab->m);
gretl_matrix_reuse(ab->kmtmp, -1, ab->m);
gretl_matrix_reuse(ab->L1, -1, ab->m);
gretl_matrix_reuse(ab->XZA, -1, ab->m);
gretl_matrix_reuse(ab->R1, ab->m, -1);
gretl_matrix_reuse(ab->ZX, ab->m, -1);
}
/* Remove zero rows/cols from A, as indicated by mask, and delete the
corresponding columns from Z. At this point we leave open the
question of whether the reduced A matrix is positive definite.
*/
static int
reduce_Z_and_A (arbond *ab, const char *mask)
{
gretl_matrix *A;
int err = 0;
A = matrix_copy_masked(ab->A, mask);
if (A == NULL) {
err = E_ALLOC;
} else {
gretl_matrix_free(ab->A);
ab->A = A;
real_shrink_matrices(ab, mask);
}
return err;
}
/* We already removed any zero rows/columns from A, but it couldn't
be inverted. Now we try reducing the dimension of A to its rank.
We need to read from the backup, Acpy, since ab->A will have been
mangled in the failed inversion attempt. We proceed to check that
the reduced A is invertible, and if not we flag an error.
*/
static int try_alt_inverse (arbond *ab)
{
char *mask = NULL;
int err = 0;
mask = make_rank_mask(ab->Acpy, &err);
if (err) {
return err;
}
gretl_matrix_free(ab->A);
ab->A = matrix_copy_masked(ab->Acpy, mask);
if (ab->A == NULL) {
free(mask);
return E_ALLOC;
}
err = gretl_invert_symmetric_matrix(ab->A);
if (!err) {
real_shrink_matrices(ab, mask);
} else {
fprintf(stderr, "try_alt_inverse: error inverting\n");
}
free(mask);
return err;
}
static int arbond_calculate (arbond *ab)
{
int err = 0;
/* find Z'X, if this job is not already done */
if (ab->doZX) {
gretl_matrix_multiply(ab->ZT, ab->dX, ab->ZX);
ab->doZX = 0;
}
/* find X'Z A */
gretl_matrix_multiply_mod(ab->ZX, GRETL_MOD_TRANSPOSE,
ab->A, GRETL_MOD_NONE,
ab->XZA, GRETL_MOD_NONE);
/* calculate "numerator", X'ZAZ'y */
gretl_matrix_multiply(ab->ZT, ab->dy, ab->R1);
gretl_matrix_multiply(ab->XZA, ab->R1, ab->beta);
/* calculate "denominator", X'ZAZ'X */
gretl_matrix_multiply(ab->XZA, ab->ZX, ab->den);
gretl_matrix_xtr_symmetric(ab->den);
gretl_matrix_copy_values(ab->kktmp, ab->den);
err = gretl_cholesky_decomp_solve(ab->kktmp, ab->beta);
#if ADEBUG
if (!err) {
int i;
fprintf(stderr, "%d-step estimates:\n\n", ab->step);
for (i=0; i<ab->k; i++) {
fprintf(stderr, "beta[%d] = %g\n", i, ab->beta->val[i]);
}
fputc('\n', stderr);
}
#endif
if (!err) {
err = arbond_variance(ab);
}
return err;
}
static int arbond_step_2 (arbond *ab, PRN *prn)
{
int err;
#if ADEBUG
gretl_matrix_print(ab->V, "V, in arbond_step_2");
#endif
if (gretl_matrix_rows(ab->V) > ab->effN) {
err = 1; /* we know this case won't work */
} else {
gretl_matrix_copy_values(ab->Acpy, ab->V);
err = gretl_invert_symmetric_matrix(ab->V);
if (err) {
gretl_matrix_copy_values(ab->V, ab->Acpy);
}
}
if (err) {
err = gretl_SVD_invert_matrix(ab->V);
if (err) {
return err;
}
gretl_matrix_xtr_symmetric(ab->V);
}
gretl_matrix_copy_values(ab->A, ab->V);
ab->step = 2;
err = arbond_calculate(ab);
if (err) {
fprintf(stderr, "step 2: arbond_calculate returned %d\n", err);
}
return err;
}
static void arbond_make_dy_and_X (arbond *ab, const double *y,
const double **Z)
{
int i, j, s, t, k = 0;
double x;
for (i=0; i<ab->N; i++) {
if (skip_unit(ab, i)) {
continue;
}
for (t=ab->ui[i].t1; t<=ab->ui[i].t2; t++) {
if (skip_obs(ab, i, t)) {
continue;
}
s = data_index(ab, i) + t;
/* current difference of dependent var */
gretl_vector_set(ab->dy, k, y[s] - y[s-1]);
for (j=0; j<ab->p; j++) {
/* lagged difference of dependent var */
gretl_matrix_set(ab->dX, k, j, y[s-j-1] - y[s-j-2]);
}
for (j=0; j<ab->nx; j++) {
/* independent vars */
x = Z[ab->xlist[j+1]][s];
gretl_matrix_set(ab->dX, k, j + ab->p, x);
}
for (j=0; j<ab->ndum; j++) {
/* time dummies */
x = (t - ab->t1min - 1 == j)? 1 : 0;
gretl_matrix_set(ab->dX, k, j + ab->p + ab->nx, x);
}
k++;
}
}
#if ADEBUG
gretl_matrix_print(ab->dy, "dy");
gretl_matrix_print(ab->dX, "X");
#endif
}
static int arbond_make_Z_and_A (arbond *ab, const double *y,
const double **Z)
{
int i, j, k, s, t, c = 0;
int zi, zj, zk;
double x;
char *zmask;
#if ADEBUG
char zstr[8];
#endif
int err = 0;
gretl_matrix_zero(ab->A);
for (i=0; i<ab->N && !err; i++) {
int ycols = ab->p; /* intial y block width */
int offj = 0; /* initialize column offset */
int Ti = unit_nobs(ab, i);
if (Ti == 0) {
continue;
}
gretl_matrix_reuse(ab->Zi, Ti, ab->m);
gretl_matrix_zero(ab->Zi);
for (t=ab->p+1; t<ab->ui[i].t1; t++) {
/* unbalanced case: compute initial column offset */
offj += ycols;
if (ycols < ab->qmax - 1) {
ycols++;
}
zj = 0;
for (zi=0; zi<ab->nzb; zi++) {
for (zk=ab->d[zi].minlag; zk<=ab->d[zi].maxlag; zk++) {
if (t - zk >= 0) {
zj++;
}
}
}
offj += zj;
}
k = 0;
for (t=ab->ui[i].t1; t<=ab->ui[i].t2; t++) {
int skip = skip_obs(ab, i, t);
int offincr = ycols;
if (!skip) {
/* lagged y (GMM instr) columns */
for (j=0; j<ycols; j++) {
s = data_index(ab, i) + t - (ycols + 1) + j;
if (!na(y[s])) {
gretl_matrix_set(ab->Zi, k, j + offj, y[s]);
}
}
}
/* additional block-diagonal columns, if required --
needs for checking for the unbalanced case
*/
zj = 0;
for (zi=0; zi<ab->nzb; zi++) {
for (zk=ab->d[zi].minlag; zk<=ab->d[zi].maxlag; zk++) {
if (t - zk >= 0) { /* ?? */
if (!skip) {
s = data_index(ab, i) + t - zk;
x = Z[ab->d[zi].v][s];
if (!na(x)) {
gretl_matrix_set(ab->Zi, k, zj + offj + ycols, x);
}
}
zj++;
}
}
}
offincr += zj;
if (!skip) {
/* additional full-length instrument columns */
s = data_index(ab, i) + t;
for (j=0; j<ab->nz; j++) {
x = Z[ab->ilist[j+1]][s];
gretl_matrix_set(ab->Zi, k, ab->xc0 + j, x);
}
/* plus time dummies, if wanted */
for (j=0; j<ab->ndum; j++) {
x = (t - ab->t1min - 1 == j)? 1 : 0;
gretl_matrix_set(ab->Zi, k, ab->xc0 + ab->nz + j, x);
}
k++; /* increment target row */
}
offj += offincr; /* starting column for next block */
if (ycols < ab->qmax - 1) {
/* increment y block width, if we're not already at max */
ycols++;
}
}
#if ADEBUG
sprintf(zstr, "Z_%d", i + 1);
gretl_matrix_print(ab->Zi, zstr);
#endif
err = make_first_diff_matrix(ab, i);
/* Cumulate Z_i' H Z_i into A_N */
gretl_matrix_qform(ab->Zi, GRETL_MOD_TRANSPOSE,
ab->H, ab->A, GRETL_MOD_CUMULATE);
/* Write Zi into ZT at offset 0, c */
gretl_matrix_inscribe_matrix(ab->ZT, ab->Zi, 0, c, GRETL_MOD_TRANSPOSE);
c += Ti;
}
/* clean up */
make_first_diff_matrix(NULL, 0);
if (!err) {
/* mask zero rows of ZT, if required */
zmask = zero_row_mask(ab->ZT, &err);
if (zmask != NULL) {
err = reduce_Z_and_A(ab, zmask);
free(zmask);
}
}
if (!err) {
gretl_matrix_divide_by_scalar(ab->A, ab->N);
}
#if ADEBUG
gretl_matrix_print(ab->ZT, "ZT");
gretl_matrix_print(ab->A, "N^{-1} * \\sum Z_i' H Z_i");
#endif
return err;
}
static int parse_diag_info (const char *s, struct diag_info *d,
const DATAINFO *pdinfo)
{
char vname[VNAMELEN];
int v, m1, m2;
int err = 0;
if (s == NULL) {
err = E_ALLOC;
} else if (sscanf(s, "GMM(%15[^, ],%d,%d)", vname, &m1, &m2) != 3) {
err = E_PARSE;
} else {
if (m2 == 0) {
m2 = 99;
}
v = varindex(pdinfo, vname);
if (v == pdinfo->v) {
err = E_UNKVAR;
} else if (m1 < 0 || m2 < m1) {
err = E_DATA;
} else {
d->v = v;
d->minlag = m1;
d->maxlag = m2;
}
}
return err;
}
static int arbond_parse_istr (const char *istr, const DATAINFO *pdinfo,
struct diag_info **pd, int *ns)
{
struct diag_info *d = NULL;
char *s0 = gretl_strdup(istr);
char *s, *p, *spec;
int i, err = 0;
if (s0 == NULL) {
return E_ALLOC;
}
/* count ')'-terminated fields */
s = s0;
while (*s) {
if (*s == ')') {
*ns += 1;
}
s++;
}
if (*ns == 0) {
err = E_PARSE;
} else {
d = malloc(*ns * sizeof *d);
if (d == NULL) {
err = E_ALLOC;
}
}
/* parse and record individual GMM instrument specs */
s = s0;
i = 0;
while (*s && !err) {
while (*s == ' ') s++;
p = s;
while (*p && *p != ')') p++;
spec = gretl_strndup(s, p - s + 1);
err = parse_diag_info(spec, &d[i++], pdinfo);
free(spec);
s = p + 1;
}
free(s0);
if (err) {
free(d);
*ns = 0;
} else {
*pd = d;
}
return err;
}
MODEL
arbond_estimate (const int *list, const char *istr, const double **X,
const DATAINFO *pdinfo, gretlopt opt,
PRN *prn)
{
struct diag_info *d = NULL;
int nzb = 0;
MODEL mod;
arbond ab;
int err = 0;
gretl_model_init(&mod);
gretl_model_smpl_init(&mod, pdinfo);
if (istr != NULL && *istr != 0) {
mod.errcode = arbond_parse_istr(istr, pdinfo, &d, &nzb);
if (mod.errcode) {
fprintf(stderr, "Error %d in arbond_parse_istr\n", mod.errcode);
return mod;
}
}
mod.errcode = arbond_init(&ab, list, pdinfo, opt, d, nzb);
if (mod.errcode) {
fprintf(stderr, "Error %d in arbond_init\n", mod.errcode);
free(d);
return mod;
}
mod.errcode = arbond_sample_check(&ab, list, X, pdinfo);
if (mod.errcode) {
fprintf(stderr, "Error %d in arbond_sample_check\n", mod.errcode);
goto bailout;
}
if (ab.nx > 0) {
arbond_zero_check(&ab, X);
}
mod.errcode = arbond_allocate(&ab);
if (mod.errcode) {
fprintf(stderr, "Error %d in arbond_allocate\n", mod.errcode);
goto bailout;
}
/* build the differenced vector \bar{y} plus the
X matrix
*/
arbond_make_dy_and_X(&ab, X[ab.yno], X);
/* build instrument matrix blocks, Z_i, and insert into
big Z' matrix; cumulate first-stage A_N as we go */
err = arbond_make_Z_and_A(&ab, X[ab.yno], X);
/* invert A_N: we back up ab.A in case inversion fails */
gretl_matrix_copy_values(ab.Acpy, ab.A);
err = gretl_invert_symmetric_matrix(ab.A);
if (err) {
fprintf(stderr, "inverting ab.A failed on first pass\n");
/* failed: try again, reducing A based on its rank */
err = try_alt_inverse(&ab);
}
#if ADEBUG
gretl_matrix_print(ab.A, "A_N");
#endif
if (!err) {
err = arbond_calculate(&ab);
}
if (!err && (opt & OPT_T)) {
err = arbond_step_2(&ab, prn);
}
bailout:
if (err && !mod.errcode) {
mod.errcode = err;
}
if (!mod.errcode) {
mod.errcode = arbond_prepare_model(&mod, &ab, list, istr,
X, pdinfo);
}
arbond_free(&ab);
return mod;
}
syntax highlighted by Code2HTML, v. 0.9.1