1 /*
2 * gretl -- Gnu Regression, Econometrics and Time-series Library
3 * Copyright (C) 2001 Allin Cottrell and Riccardo "Jack" Lucchetti
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, either version 3 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program. If not, see <http://www.gnu.org/licenses/>.
17 *
18 */
19
20 #include "libgretl.h"
21 #include "qr_estimate.h"
22 #include "gretl_matrix.h"
23 #include "matrix_extra.h"
24 #include "libset.h"
25 #include "gretl_panel.h"
26 #include "estim_private.h"
27
28 #include "gretl_f2c.h"
29 #include "clapack_double.h"
30
31 #define ESSZERO 1e-22 /* SSR less than this counts as zero */
32
33 enum {
34 VCV_SIMPLE,
35 VCV_ROBUST,
36 VCV_XPX
37 };
38
39 static int qr_make_cluster_vcv (MODEL *pmod, int ci,
40 const DATASET *dset,
41 gretl_matrix *XX,
42 gretlopt opt);
43
44 #define gls_rho(p) gretl_model_get_double_default(p, "rho_gls", 0.0)
45
46 /* General note: in fortran arrays, column entries are contiguous.
47 Columns of data matrix X hold variables, rows hold observations.
48 So in a fortran array, entries for a given variable are contiguous.
49 */
50
qr_get_tss(MODEL * pmod,const DATASET * dset,int * ifc,int * yconst)51 static double qr_get_tss (MODEL *pmod, const DATASET *dset,
52 int *ifc, int *yconst)
53 {
54 int yno = pmod->list[1];
55 int pwe = (pmod->opt & OPT_P);
56 double rho = gls_rho(pmod);
57 double yt, y0 = 0.0, ymean = 0.0;
58 double x, tss = 0.0;
59 double ctss = 0.0;
60 int t;
61
62 if (*ifc == 0) {
63 *ifc = check_for_effective_const(pmod, dset);
64 }
65
66 *yconst = 1;
67
68 if (rho != 0.0) {
69 double ry, d;
70
71 for (t=pmod->t1; t<=pmod->t2; t++) {
72 ry = dset->Z[yno][t];
73 if (t == pmod->t1 && pwe) {
74 ry *= sqrt(1.0 - rho * rho);
75 } else {
76 ry -= rho * dset->Z[yno][t-1];
77 }
78 ymean += ry;
79 }
80 ymean /= pmod->nobs;
81
82 for (t=pmod->t1; t<=pmod->t2; t++) {
83 ry = dset->Z[yno][t];
84 if (t == pmod->t1 && pwe) {
85 ry *= sqrt(1.0 - rho * rho);
86 } else {
87 ry -= rho * dset->Z[yno][t-1];
88 }
89 if (t == pmod->t1) {
90 y0 = ry;
91 } else if (ry != y0) {
92 *yconst = 0;
93 }
94 d = ry - ymean;
95 if (*ifc) {
96 tss += d * d;
97 } else {
98 tss += ry * ry;
99 ctss += d * d;
100 }
101 }
102 } else {
103 for (t=pmod->t1; t<=pmod->t2; t++) {
104 if (!na(pmod->yhat[t])) {
105 ymean += dset->Z[yno][t];
106 }
107 }
108 ymean /= pmod->nobs;
109
110 y0 = dset->Z[yno][pmod->t1];
111
112 for (t=pmod->t1; t<=pmod->t2; t++) {
113 if (!na(pmod->yhat[t])) {
114 if (dset->Z[yno][t] != y0) {
115 *yconst = 0;
116 }
117 x = dset->Z[yno][t] - ymean;
118 if (*ifc) {
119 tss += x * x;
120 } else {
121 yt = dset->Z[yno][t];
122 tss += yt * yt;
123 ctss += x * x;
124 }
125 }
126 }
127 }
128
129 if (!*ifc && ctss > 0) {
130 double cR2 = 1 - (pmod->ess / ctss);
131
132 gretl_model_set_double(pmod, "centered-R2", cR2);
133 }
134
135 return tss;
136 }
137
qr_compute_stats(MODEL * pmod,const DATASET * dset,int n,gretlopt opt)138 static void qr_compute_stats (MODEL *pmod, const DATASET *dset,
139 int n, gretlopt opt)
140 {
141 int yconst, ifc = pmod->ifc;
142 int yno = pmod->list[1];
143
144 pmod->tss = qr_get_tss(pmod, dset, &ifc, &yconst);
145 pmod->chisq = NADBL;
146
147 if (yconst && pmod->dfd > 0) {
148 double y0 = dset->Z[yno][pmod->t1];
149
150 if (y0 > 0) {
151 double tss = pmod->nobs * y0 * y0;
152
153 pmod->rsq = 1 - (pmod->ess / tss);
154 gretl_model_set_int(pmod, "uncentered", 1);
155 }
156 } else if (pmod->dfd > 0) {
157 double den = pmod->tss * pmod->dfd;
158
159 pmod->rsq = 1.0 - (pmod->ess / pmod->tss);
160 pmod->adjrsq = 1.0 - (pmod->ess * (n - 1) / den);
161 } else {
162 pmod->rsq = 1.0;
163 }
164
165 if (pmod->ncoeff == 1 && pmod->ifc) {
166 pmod->fstt = NADBL;
167 return;
168 }
169
170 if (pmod->dfd > 0 && pmod->dfn > 0) {
171 if (opt & OPT_R) {
172 pmod->fstt = wald_omit_F(NULL, pmod);
173 } else if (pmod->rsq == 1.0) {
174 pmod->fstt = NADBL;
175 } else {
176 pmod->fstt = (pmod->tss - pmod->ess) * pmod->dfd /
177 (pmod->ess * pmod->dfn);
178 }
179 } else {
180 pmod->fstt = NADBL;
181 }
182 }
183
qr_make_vcv(MODEL * pmod,gretl_matrix * V,int flag)184 static int qr_make_vcv (MODEL *pmod, gretl_matrix *V, int flag)
185 {
186 int k = pmod->ncoeff;
187 int m = k * (k + 1) / 2;
188 double x;
189 int i, j, idx;
190
191 pmod->vcv = malloc(m * sizeof *pmod->vcv);
192 if (pmod->vcv == NULL) {
193 return E_ALLOC;
194 }
195
196 if (flag == VCV_XPX) {
197 gretl_model_set_int(pmod, "vcv_xpx", 1);
198 }
199
200 for (i=0; i<k; i++) {
201 for (j=0; j<=i; j++) {
202 idx = ijton(i, j, k);
203 x = gretl_matrix_get(V, i, j);
204 if (flag == VCV_SIMPLE) {
205 x *= pmod->sigma * pmod->sigma;
206 }
207 pmod->vcv[idx] = x;
208 if (i == j) {
209 pmod->sderr[i] = sqrt(x);
210 if (flag == VCV_XPX) {
211 pmod->sderr[i] *= pmod->sigma;
212 }
213 }
214 }
215 }
216
217 return 0;
218 }
219
get_resids_and_SSR(MODEL * pmod,const DATASET * dset,gretl_matrix * yhat,int fulln)220 static void get_resids_and_SSR (MODEL *pmod, const DATASET *dset,
221 gretl_matrix *yhat, int fulln)
222 {
223 double rho = gls_rho(pmod);
224 int qdiff = (rho != 0.0);
225 int pwe = (pmod->opt & OPT_P);
226 int yvar = pmod->list[1];
227 double *u = pmod->uhat;
228 double yt;
229 int t, i = 0;
230
231 pmod->ess = 0.0;
232
233 if (qdiff) {
234 for (t=0; t<fulln; t++) {
235 if (t < pmod->t1 || t > pmod->t2) {
236 pmod->yhat[t] = u[t] = NADBL;
237 } else {
238 yt = dset->Z[yvar][t];
239 if (t == pmod->t1 && pwe) {
240 yt *= sqrt(1.0 - rho * rho);
241 } else {
242 yt -= rho * dset->Z[yvar][t-1];
243 }
244 pmod->yhat[t] = yhat->val[i];
245 u[t] = yt - yhat->val[i];
246 pmod->ess += u[t] * u[t];
247 i++;
248 }
249 }
250 } else if (pmod->nwt) {
251 double wt;
252
253 for (t=0; t<fulln; t++) {
254 if (t < pmod->t1 || t > pmod->t2 || model_missing(pmod, t)) {
255 pmod->yhat[t] = u[t] = NADBL;
256 } else {
257 wt = dset->Z[pmod->nwt][t];
258 yt = dset->Z[yvar][t];
259 if (wt == 0.0 || na(wt)) {
260 pmod->yhat[t] = NADBL;
261 } else {
262 yt *= sqrt(dset->Z[pmod->nwt][t]);
263 pmod->yhat[t] = yhat->val[i];
264 u[t] = yt - yhat->val[i];
265 pmod->ess += u[t] * u[t];
266 i++;
267 }
268 }
269 }
270 } else {
271 for (t=0; t<fulln; t++) {
272 if (t < pmod->t1 || t > pmod->t2 || model_missing(pmod, t)) {
273 pmod->yhat[t] = u[t] = NADBL;
274 } else {
275 pmod->yhat[t] = yhat->val[i];
276 u[t] = dset->Z[yvar][t] - yhat->val[i];
277 pmod->ess += u[t] * u[t];
278 i++;
279 }
280 }
281 }
282
283 /* if SSR is small enough, treat it as zero */
284 if (fabs(pmod->ess) < ESSZERO) {
285 pmod->ess = 0.0;
286 }
287 }
288
289 static void
get_data_X(gretl_matrix * X,const MODEL * pmod,const DATASET * dset)290 get_data_X (gretl_matrix *X, const MODEL *pmod, const DATASET *dset)
291 {
292 int wt = pmod->nwt;
293 int i, t, vi, s = 0;
294
295 /* copy independent vars into matrix X */
296
297 if (!wt && pmod->missmask == NULL) {
298 /* use faster method */
299 int T = pmod->t2 - pmod->t1 + 1;
300 size_t sz = T * sizeof(double);
301
302 for (i=2; i<=pmod->list[0]; i++) {
303 vi = pmod->list[i];
304 memcpy(X->val + s, dset->Z[vi] + pmod->t1, sz);
305 s += T;
306 }
307 } else {
308 for (i=2; i<=pmod->list[0]; i++) {
309 vi = pmod->list[i];
310 for (t=pmod->t1; t<=pmod->t2; t++) {
311 if (!model_missing(pmod, t)) {
312 if (wt) {
313 X->val[s++] = sqrt(dset->Z[wt][t]) * dset->Z[vi][t];
314 } else {
315 X->val[s++] = dset->Z[vi][t];
316 }
317 }
318 }
319 }
320 }
321 }
322
make_data_X(const MODEL * pmod,const DATASET * dset)323 static gretl_matrix *make_data_X (const MODEL *pmod,
324 const DATASET *dset)
325 {
326 gretl_matrix *X;
327
328 X = gretl_matrix_alloc(pmod->nobs, pmod->ncoeff);
329 if (X != NULL) {
330 get_data_X(X, pmod, dset);
331 }
332
333 return X;
334 }
335
336 /* Calculate W(t)-transpose * W(t-lag) */
337
wtw(gretl_matrix * wt,const gretl_matrix * H,int t,int lag)338 static void wtw (gretl_matrix *wt, const gretl_matrix *H,
339 int t, int lag)
340 {
341 int i, j, s = t - lag;
342 double hti, htj;
343
344 for (i=0; i<H->cols; i++) {
345 hti = gretl_matrix_get(H, t, i);
346 for (j=0; j<H->cols; j++) {
347 htj = gretl_matrix_get(H, s, j);
348 gretl_matrix_set(wt, i, j, hti * htj);
349 }
350 }
351 }
352
353 /* special handling for quadratic spectral kernel */
354
qs_hac_weight(double bt,int i)355 double qs_hac_weight (double bt, int i)
356 {
357 double di = i / bt;
358 double mi = 6 * M_PI * di / 5;
359 double w;
360
361 w = 25 / (12 * M_PI * M_PI * di * di);
362 w *= sin(mi) / mi - cos(mi);
363
364 return w;
365 }
366
hac_weight(int kern,int h,int i)367 double hac_weight (int kern, int h, int i)
368 {
369 double ai = fabs((double) i) / (h + 1.0);
370 double w;
371
372 if (kern == KERNEL_PARZEN) {
373 if (ai <= 0.5) {
374 w = 1.0 - 6*ai*ai + 6*pow(ai, 3.0);
375 } else {
376 w = 2.0 * pow(1.0 - ai, 3.0);
377 }
378 } else {
379 /* Bartlett kernel */
380 w = 1.0 - ai;
381 }
382
383 return w;
384 }
385
386 #define NW_DEBUG 0
387
388 /* Newey and West's parameter 'n' for truncation when
389 calculating s^(0) and s^(1) in the course of doing
390 the data-determined bandwidth thing.
391 */
392
lag_trunc_param(int kern,int prewhitened,int T)393 static int lag_trunc_param (int kern, int prewhitened, int T)
394 {
395 if (kern == KERNEL_BARTLETT) {
396 /* Newey and West (1994) */
397 if (prewhitened) {
398 return (int) (3.0 * pow(T/100.0, 2.0 / 9));
399 } else {
400 return (int) (4.0 * pow(T/100.0, 2.0 / 9));
401 }
402 /* Hall just gives "O(T^{2/9})" */
403 /* return (int) pow((double) T, 2.0 / 9); */
404 } else if (kern == KERNEL_PARZEN) {
405 return (int) pow((double) T, 4.0 / 25);
406 } else {
407 return (int) pow((double) T, 2.0 / 25);
408 }
409 }
410
411 /* Newey and West's data-based bandwidth selection, based on
412 the exposition in their 1994 paper in Review of Economic
413 Studies, pp. 634-5. See also A. Hall, "Generalized Method
414 of Moments" (Oxford, 2005), p. 82.
415
416 Public because it's called from gmm.c (with @w NULL).
417 */
418
newey_west_bandwidth(const gretl_matrix * H,const gretl_matrix * w,int kern,int prewhitened,int * m,double * b)419 int newey_west_bandwidth (const gretl_matrix *H,
420 const gretl_matrix *w,
421 int kern, int prewhitened,
422 int *m, double *b)
423 {
424 /* kernel-specific parameter, c_\gamma */
425 const double cg[] = { 1.1447, 2.6614, 1.3221 };
426 double sigma;
427 double htw0, htwj, wi;
428 double s1, s0, g, p;
429 int T = H->rows;
430 int k = H->cols;
431 int n, i, j, t;
432 int err = 0;
433
434 n = lag_trunc_param(kern, prewhitened, T);
435
436 /* calculate sigma_0 */
437 sigma = 0;
438 for (t=1; t<T; t++) {
439 htw0 = 0.0;
440 for (i=0; i<k; i++) {
441 wi = w != NULL ? w->val[i] : 1.0;
442 htw0 += gretl_matrix_get(H, t, i) * wi;
443 }
444 sigma += htw0 * htw0;
445 }
446 sigma /= T-1;
447
448 /* initialize \hat{s}^{(1)} and \hat{s}^{(0)} */
449 s1 = 0;
450 s0 = sigma;
451
452 #if NW_DEBUG
453 fprintf(stderr, "T = %d, sigma0 = %g\n", T, sigma);
454 #endif
455
456 /* calculate \sigma_1 to \sigma_n and cumulate the \hat{s} values */
457 for (j=1; j<=n; j++) {
458 sigma = 0;
459 for (t=j+1; t<T; t++) {
460 htw0 = htwj = 0.0;
461 for (i=0; i<k; i++) {
462 wi = w != NULL ? w->val[i] : 1.0;
463 htw0 += gretl_matrix_get(H, t, i) * wi;
464 htwj += gretl_matrix_get(H, t-j, i) * wi;
465 }
466 sigma += htw0 * htwj;
467 }
468 sigma /= T-1;
469 if (kern == KERNEL_BARTLETT) {
470 s1 += 2*j*sigma;
471 } else {
472 s1 += 2*j*j*sigma;
473 }
474 s0 += 2 * sigma;
475 #if NW_DEBUG
476 fprintf(stderr, "sigma[%d] = %g\n", j, sigma);
477 #endif
478 }
479
480 p = (kern == KERNEL_BARTLETT)? 1.0/3.0 : 1.0/5.0;
481 g = cg[kern] * pow((s1 / s0) * (s1 / s0), p);
482 *b = g * pow((double) T, p);
483 *m = (int) floor(*b);
484
485 #if NW_DEBUG
486 fprintf(stderr, "s1=%g, s0=%g, gamma=%g, b=%g, m=%d\n",
487 s1, s0, g, *b, *m);
488 #endif
489
490 if (*b > T / 2.0) {
491 /* FIXME arbitrary truncation in case this method has gone
492 wonky! */
493 fprintf(stderr, "newey_west_bandwidth (PW=%d): invalid result %d (s^(0)=%g)\n",
494 prewhitened, *m, s0);
495 *b = T / 2.0;
496 *m = (int) floor(*b);
497 }
498
499 return err;
500 }
501
hac_recolor(gretl_matrix * XOX,gretl_matrix * A)502 static int hac_recolor (gretl_matrix *XOX, gretl_matrix *A)
503 {
504 gretl_matrix *S;
505 int k = XOX->rows;
506 int err = 0;
507
508 S = gretl_matrix_alloc(k, k);
509 if (S == NULL) {
510 err = E_ALLOC;
511 }
512
513 if (!err) {
514 double aij;
515 int i, j;
516
517 /* we won't need A hereafter, so make A into (I-A) */
518 for (i=0; i<k; i++) {
519 for (j=0; j<k; j++) {
520 aij = gretl_matrix_get(A, i, j);
521 aij = (i == j)? 1 - aij: -aij;
522 gretl_matrix_set(A, i, j, aij);
523 }
524 }
525 err = gretl_invert_matrix(A);
526 }
527
528 if (!err) {
529 /* NW: "S = (I-A)^{-1} * XOX * (I-A)^{-1}'", but note
530 that our A is the transpose of Newey and West's
531 */
532 err = gretl_matrix_qform(A, GRETL_MOD_TRANSPOSE,
533 XOX, S, GRETL_MOD_NONE);
534 }
535
536 if (!err) {
537 gretl_matrix_copy_values(XOX, S);
538 }
539
540 gretl_matrix_free(S);
541
542 return err;
543 }
544
revise_VAR_residuals(gretl_matrix * A,gretl_matrix * Y,gretl_matrix * X,gretl_matrix * E)545 static int revise_VAR_residuals (gretl_matrix *A,
546 gretl_matrix *Y,
547 gretl_matrix *X,
548 gretl_matrix *E)
549 {
550 gretl_matrix *Xt, *Pt;
551 double yti;
552 int k = X->cols;
553 int T = X->rows;
554 int i, t;
555
556 Xt = gretl_matrix_alloc(1, k);
557 Pt = gretl_matrix_alloc(1, k);
558
559 if (Xt == NULL || Pt == NULL) {
560 return E_ALLOC;
561 }
562
563 for (t=0; t<T; t++) {
564 for (i=0; i<k; i++) {
565 Xt->val[i] = gretl_matrix_get(X, t, i);
566 }
567 gretl_matrix_multiply(Xt, A, Pt);
568 for (i=0; i<k; i++) {
569 yti = gretl_matrix_get(Y, t, i);
570 gretl_matrix_set(E, t, i, yti - Pt->val[i]);
571 }
572 }
573
574 gretl_matrix_free(Xt);
575 gretl_matrix_free(Pt);
576
577 return 0;
578 }
579
maybe_limit_VAR_coeffs(gretl_matrix * A,gretl_matrix * Y,gretl_matrix * X,gretl_matrix * E)580 int maybe_limit_VAR_coeffs (gretl_matrix *A,
581 gretl_matrix *Y,
582 gretl_matrix *X,
583 gretl_matrix *E)
584 {
585 gretl_matrix *B = NULL;
586 gretl_matrix *D = NULL;
587 gretl_matrix *C = NULL;
588 gretl_matrix *T = NULL;
589 int i, k = A->rows;
590 int amod = 0;
591 int err;
592
593 err = gretl_matrix_SVD(A, &B, &D, &C, 0);
594
595 if (!err) {
596 for (i=0; i<k; i++) {
597 if (D->val[i] > 0.97) {
598 D->val[i] = 0.97;
599 amod = 1;
600 }
601 }
602 if (amod) {
603 /* "shrink" A */
604 T = gretl_matrix_dot_op(B, D, '*', &err);
605 gretl_matrix_multiply(T, C, A);
606 }
607 }
608
609 if (amod && X != NULL && Y != NULL && E != NULL) {
610 err = revise_VAR_residuals(A, Y, X, E);
611 }
612
613 gretl_matrix_free(B);
614 gretl_matrix_free(D);
615 gretl_matrix_free(C);
616 gretl_matrix_free(T);
617
618 return err;
619 }
620
nw_prewhiten(gretl_matrix * H,gretl_matrix ** pA)621 static int nw_prewhiten (gretl_matrix *H, gretl_matrix **pA)
622 {
623 gretl_matrix *H0, *H1, *A, *U;
624 gretl_matrix_block *B;
625 int T = H->rows;
626 int k = H->cols;
627 double h;
628 int t, j;
629 int err = 0;
630
631 B = gretl_matrix_block_new(&H0, T-1, k,
632 &H1, T-1, k,
633 &A, k, k,
634 &U, T-1, k,
635 NULL);
636 if (B == NULL) {
637 return E_ALLOC;
638 }
639
640 for (j=0; j<k; j++) {
641 for (t=1; t<T; t++) {
642 h = gretl_matrix_get(H, t, j);
643 gretl_matrix_set(H0, t-1, j, h);
644 h = gretl_matrix_get(H, t-1, j);
645 gretl_matrix_set(H1, t-1, j, h);
646 }
647 }
648
649 /* estimate VAR */
650 err = gretl_matrix_multi_ols(H0, H1, A, U, NULL);
651
652 /* apply the "0.97 limit" if needed */
653 if (!err) {
654 err = maybe_limit_VAR_coeffs(A, H0, H1, U);
655 }
656
657 if (!err) {
658 /* replace incoming H with VAR residuals */
659 for (j=0; j<k; j++) {
660 gretl_matrix_set(H, 0, j, 0.0);
661 for (t=1; t<T; t++) {
662 h = gretl_matrix_get(U, t-1, j);
663 gretl_matrix_set(H, t, j, h);
664 }
665 }
666 }
667
668 #if NW_DEBUG
669 gretl_matrix_print(A, "prewhite A-hat");
670 gretl_matrix_print(H, "prewhitened H");
671 #endif
672
673 if (!err) {
674 *pA = gretl_matrix_copy(A);
675 if (*pA == NULL) {
676 err = E_ALLOC;
677 }
678 }
679
680 gretl_matrix_block_destroy(B);
681
682 return err;
683 }
684
685 /* Form the matrix H, such that H_t = X_t * u_t. In
686 addition, if @pw is non-NULL and the first column
687 of @X is constant, then write into @pw a vector of
688 1s with the first element set to zero. Or if @u is
689 NULL, just copy X to H, but again fill @pw if
690 required.
691 */
692
newey_west_H(const gretl_matrix * X,const gretl_matrix * u,gretl_matrix ** pw,int * err)693 static gretl_matrix *newey_west_H (const gretl_matrix *X,
694 const gretl_matrix *u,
695 gretl_matrix **pw,
696 int *err)
697 {
698 gretl_matrix *H;
699 double x0 = X->val[0];
700 int T = X->rows;
701 int k = X->cols;
702 int make_w;
703 int j, t;
704
705 if (u == NULL) {
706 H = gretl_matrix_copy(X);
707 } else {
708 H = gretl_matrix_alloc(T, k);
709 }
710
711 if (H == NULL) {
712 *err = E_ALLOC;
713 return NULL;
714 }
715
716 /* tentative */
717 make_w = (pw != NULL && k > 1);
718
719 if (u == NULL) {
720 /* we just have to check for constancy */
721 if (make_w) {
722 for (t=1; t<T; t++) {
723 if (X->val[t] != x0) {
724 make_w = 0;
725 break;
726 }
727 }
728 }
729 } else {
730 /* @u is non-NULL */
731 double xtj;
732
733 for (j=0; j<k; j++) {
734 for (t=0; t<T; t++) {
735 xtj = gretl_matrix_get(X, t, j);
736 gretl_matrix_set(H, t, j, xtj * u->val[t]);
737 if (make_w && j == 0 && t > 0 && xtj != x0) {
738 /* first X column not constant */
739 make_w = 0;
740 }
741 }
742 }
743 }
744
745 if (make_w) {
746 *pw = gretl_unit_matrix_new(k, 1);
747 if (*pw != NULL) {
748 (*pw)->val[0] = 0;
749 }
750 }
751
752 return H;
753 }
754
755 /* HAC_XOX: compute the "sandwich filling" for the HAC estimator */
756
HAC_XOX(const gretl_matrix * X,const gretl_matrix * uhat,VCVInfo * vi,int use_prior,int * err)757 gretl_matrix *HAC_XOX (const gretl_matrix *X,
758 const gretl_matrix *uhat,
759 VCVInfo *vi, int use_prior,
760 int *err)
761 {
762 gretl_matrix *XOX = NULL;
763 gretl_matrix *Wtj = NULL;
764 gretl_matrix *Gj = NULL;
765 gretl_matrix *H = NULL;
766 gretl_matrix *A = NULL;
767 gretl_matrix *w = NULL;
768 int prewhiten;
769 int data_based;
770 int kern;
771 int T = X->rows;
772 int k = X->cols;
773 int p, j, t;
774 double bt = 0;
775
776 if (use_prior) {
777 kern = vi->vmin;
778 prewhiten = vi->flags & HAC_PREWHITEN;
779 if (kern == KERNEL_QS) {
780 bt = vi->bw;
781 p = T - 1;
782 } else {
783 p = vi->order;
784 }
785 data_based = 0;
786 } else {
787 kern = libset_get_int(HAC_KERNEL);
788 data_based = data_based_hac_bandwidth();
789 prewhiten = libset_get_bool(PREWHITEN);
790 }
791
792 #if NW_DEBUG
793 fprintf(stderr, "*** HAC: kern = %d, prewhiten = %d ***\n",
794 kern, prewhiten);
795 #endif
796
797 if (use_prior) {
798 H = newey_west_H(X, uhat, NULL, err);
799 } else if (prewhiten || data_based) {
800 H = newey_west_H(X, uhat, &w, err);
801 } else {
802 H = newey_west_H(X, uhat, NULL, err);
803 }
804
805 if (H == NULL) {
806 return NULL;
807 } else if (prewhiten) {
808 *err = nw_prewhiten(H, &A);
809 }
810
811 if (!*err) {
812 XOX = gretl_zero_matrix_new(k, k);
813 Wtj = gretl_matrix_alloc(k, k);
814 Gj = gretl_matrix_alloc(k, k);
815 if (XOX == NULL || Wtj == NULL || Gj == NULL) {
816 *err = E_ALLOC;
817 }
818 }
819
820 if (*err) {
821 goto bailout;
822 }
823
824 /* determine the bandwidth setting */
825
826 if (!use_prior) {
827 if (data_based) {
828 *err = newey_west_bandwidth(H, w, kern, prewhiten, &p, &bt);
829 if (*err) {
830 goto bailout;
831 }
832 } else if (kern == KERNEL_QS) {
833 bt = libset_get_double(QS_BANDWIDTH);
834 p = T - 1;
835 } else {
836 p = get_hac_lag(T);
837 }
838 }
839
840 if (!*err) {
841 double wj;
842
843 for (j=0; j<=p; j++) {
844 /* cumulate running sum of Gamma-hat terms */
845 gretl_matrix_zero(Gj);
846 for (t=j; t<T; t++) {
847 /* W(t)-transpose * W(t-j) */
848 wtw(Wtj, H, t, j);
849 gretl_matrix_add_to(Gj, Wtj);
850 }
851 if (j > 0) {
852 /* Gamma(j) = Gamma(j) + Gamma(j)-transpose */
853 gretl_matrix_add_self_transpose(Gj);
854 if (kern == KERNEL_QS) {
855 wj = qs_hac_weight(bt, j);
856 } else {
857 wj = hac_weight(kern, p, j);
858 }
859 gretl_matrix_multiply_by_scalar(Gj, wj);
860 }
861 gretl_matrix_add_to(XOX, Gj);
862 }
863 }
864
865 if (A != NULL) {
866 hac_recolor(XOX, A);
867 }
868
869 if (vi != NULL && !use_prior) {
870 vi->vmaj = VCV_HAC;
871 vi->vmin = kern;
872 vi->flags = prewhiten ? HAC_PREWHITEN : 0;
873
874 if (kern == KERNEL_QS) {
875 vi->order = 0;
876 vi->bw = bt;
877 } else {
878 vi->order = p;
879 vi->bw = NADBL;
880 }
881 }
882
883 bailout:
884
885 gretl_matrix_free(H);
886 gretl_matrix_free(Wtj);
887 gretl_matrix_free(Gj);
888 gretl_matrix_free(A);
889 gretl_matrix_free(w);
890
891 if (*err && XOX != NULL) {
892 gretl_matrix_free(XOX);
893 XOX = NULL;
894 }
895
896 return XOX;
897 }
898
newey_west_OPG(const gretl_matrix * G,int * err)899 gretl_matrix *newey_west_OPG (const gretl_matrix *G, int *err)
900 {
901 return HAC_XOX(G, NULL, NULL, 0, err);
902 }
903
904 /* To support the hansl function lrcovar(): note that
905 setting @demean to non-zero value means that @X should
906 be demeaned before proceeding.
907 */
908
long_run_covariance(const gretl_matrix * X,int demean,int * err)909 gretl_matrix *long_run_covariance (const gretl_matrix *X,
910 int demean, int *err)
911 {
912 gretl_matrix *V = NULL;
913
914 if (demean) {
915 gretl_matrix *Xd = gretl_matrix_copy(X);
916
917 if (Xd == NULL) {
918 *err = E_ALLOC;
919 } else {
920 gretl_matrix_center(Xd);
921 V = HAC_XOX(Xd, NULL, NULL, 0, err);
922 gretl_matrix_free(Xd);
923 }
924 } else {
925 V = HAC_XOX(X, NULL, NULL, 0, err);
926 }
927
928 if (V != NULL) {
929 gretl_matrix_divide_by_scalar(V, X->rows);
930 }
931
932 return V;
933 }
934
935 /* not a rigorous check, but catches the worst cases */
936
vcv_is_broken(const gretl_matrix * V)937 static int vcv_is_broken (const gretl_matrix *V)
938 {
939 if (V == NULL) {
940 return 1;
941 } else {
942 double v;
943 int i;
944
945 for (i=0; i<V->rows; i++) {
946 v = gretl_matrix_get(V, i, i);
947 if (v < 0 || na(v)) {
948 return 1;
949 }
950 }
951 }
952
953 return 0;
954 }
955
956 /* Calculate HAC covariance matrix. Algorithm and (basically)
957 notation taken from Davidson and MacKinnon (DM), Econometric
958 Theory and Methods, chapter 9.
959 */
960
qr_make_hac(MODEL * pmod,const DATASET * dset,gretl_matrix * XTXi)961 static int qr_make_hac (MODEL *pmod, const DATASET *dset,
962 gretl_matrix *XTXi)
963 {
964 gretl_matrix *X, *XOX, *V = NULL;
965 gretl_matrix umat;
966 VCVInfo vi;
967 int T = pmod->nobs;
968 int err = 0;
969
970 X = make_data_X(pmod, dset);
971 if (X == NULL) {
972 return E_ALLOC;
973 }
974
975 /* pmod->uhat is a full-length series: we must take an offset
976 into it, equal to the offset of the data on which the model
977 is actually estimated.
978 */
979 gretl_matrix_init(&umat);
980 umat.rows = T;
981 umat.cols = 1;
982 umat.val = pmod->uhat + pmod->t1;
983
984 XOX = HAC_XOX(X, &umat, &vi, 0, &err);
985
986 if (!err) {
987 V = gretl_matrix_alloc(XOX->rows, XOX->rows);
988 if (V == NULL) {
989 err = E_ALLOC;
990 }
991 }
992
993 if (!err) {
994 gretl_matrix_qform(XTXi, GRETL_MOD_TRANSPOSE, XOX,
995 V, GRETL_MOD_NONE);
996 if (vcv_is_broken(V)) {
997 err = E_JACOBIAN;
998 } else {
999 /* Transcribe vcv into triangular representation */
1000 err = qr_make_vcv(pmod, V, VCV_ROBUST);
1001 }
1002 }
1003
1004 if (!err) {
1005 gretl_model_set_full_vcv_info(pmod, vi.vmaj, vi.vmin,
1006 vi.order, vi.flags,
1007 vi.bw);
1008 if (!na(vi.bw)) {
1009 gretl_model_set_double(pmod, "hac_bw", vi.bw);
1010 } else {
1011 gretl_model_set_double(pmod, "hac_bw", (double) vi.order);
1012 }
1013 }
1014
1015 gretl_matrix_free(X);
1016 gretl_matrix_free(XOX);
1017 gretl_matrix_free(V);
1018
1019 return err;
1020 }
1021
1022 /* Multiply X transpose into D, treating D as if it were
1023 a diagonal matrix -- although in fact it is just a vector,
1024 for economy of storage. Result into R.
1025 */
1026
do_X_prime_diag(const gretl_matrix * X,const gretl_vector * D,gretl_matrix * R)1027 static void do_X_prime_diag (const gretl_matrix *X,
1028 const gretl_vector *D,
1029 gretl_matrix *R)
1030 {
1031 double x;
1032 int i, j;
1033
1034 for (i=0; i<R->rows; i++) {
1035 for (j=0; j<R->cols; j++) {
1036 x = gretl_matrix_get(X, j, i);
1037 gretl_matrix_set(R, i, j, x * D->val[j]);
1038 }
1039 }
1040 }
1041
1042 /* Heteroskedasticity-Consistent Covariance Matrices: See Davidson
1043 and MacKinnon, Econometric Theory and Methods, chapter 5, esp.
1044 page 200. Implements HC0, HC1, HC2 and HC3. Also added: the
1045 jackknife as described in MacKinnon and White (Journal of
1046 Econometrics, 1985), or "HC3a" in gretl parlance.
1047 */
1048
qr_make_hccme(MODEL * pmod,const DATASET * dset,gretl_matrix * Q,gretl_matrix * XTXi)1049 static int qr_make_hccme (MODEL *pmod, const DATASET *dset,
1050 gretl_matrix *Q, gretl_matrix *XTXi)
1051 {
1052 gretl_matrix *X;
1053 gretl_matrix *diag = NULL;
1054 gretl_matrix *tmp1 = NULL;
1055 gretl_matrix *tmp2 = NULL;
1056 gretl_matrix *V = NULL;
1057 gretl_matrix *u = NULL;
1058 int T = pmod->nobs;
1059 int k = pmod->list[0] - 1;
1060 int hc_version;
1061 int i, t;
1062 int err = 0;
1063
1064 X = make_data_X(pmod, dset);
1065 if (X == NULL) return 1;
1066
1067 hc_version = libset_get_int(HC_VERSION);
1068
1069 diag = gretl_vector_from_array(pmod->uhat + pmod->t1, T,
1070 GRETL_MOD_SQUARE);
1071 if (diag == NULL) {
1072 err = 1;
1073 goto bailout;
1074 }
1075 if (hc_version == 4) {
1076 u = gretl_vector_from_array(pmod->uhat + pmod->t1, T,
1077 GRETL_MOD_NONE);
1078 if (u == NULL) {
1079 err = 1;
1080 goto bailout;
1081 }
1082 }
1083
1084 tmp1 = gretl_matrix_alloc(k, T);
1085 tmp2 = gretl_matrix_alloc(k, k);
1086 V = gretl_matrix_alloc(k, k);
1087 if (tmp1 == NULL || tmp2 == NULL || V == NULL) {
1088 err = 1;
1089 goto bailout;
1090 }
1091
1092 gretl_model_set_vcv_info(pmod, VCV_HC, hc_version);
1093
1094 if (hc_version == 1) {
1095 for (t=0; t<T; t++) {
1096 diag->val[t] *= (double) T / (T - k);
1097 }
1098 } else if (hc_version > 1) {
1099 /* do the h_t calculations */
1100 for (t=0; t<T; t++) {
1101 double q, ht = 0.0;
1102
1103 for (i=0; i<k; i++) {
1104 q = gretl_matrix_get(Q, t, i);
1105 ht += q * q;
1106 }
1107 if (hc_version == 2) {
1108 diag->val[t] /= (1.0 - ht);
1109 } else {
1110 /* HC3+ */
1111 diag->val[t] /= (1.0 - ht) * (1.0 - ht);
1112 }
1113 if (hc_version == 4) {
1114 u->val[t] /= (1.0 - ht);
1115 }
1116 }
1117 }
1118
1119 do_X_prime_diag(X, diag, tmp1);
1120 gretl_matrix_multiply(tmp1, X, tmp2);
1121 gretl_matrix_qform(XTXi, GRETL_MOD_NONE, tmp2,
1122 V, GRETL_MOD_NONE);
1123
1124 if (hc_version == 4) {
1125 /* jackknife, per MacKinnon and White (1985) */
1126 double nfac = (T - 1.0) / T;
1127 gretl_matrix *g;
1128
1129 gretl_matrix_multiply_mod(XTXi, GRETL_MOD_NONE,
1130 X, GRETL_MOD_TRANSPOSE,
1131 tmp1, GRETL_MOD_NONE);
1132 g = gretl_matrix_reuse(diag, k, 1);
1133 gretl_matrix_multiply(tmp1, u, g);
1134 /* \gamma * \gamma' */
1135 gretl_matrix_multiply_mod(g, GRETL_MOD_NONE,
1136 g, GRETL_MOD_TRANSPOSE,
1137 tmp2, GRETL_MOD_NONE);
1138 gretl_matrix_divide_by_scalar(tmp2, T);
1139 gretl_matrix_subtract_from(V, tmp2);
1140 gretl_matrix_multiply_by_scalar(V, nfac);
1141 }
1142
1143 /* Transcribe vcv into triangular representation */
1144 err = qr_make_vcv(pmod, V, VCV_ROBUST);
1145
1146 bailout:
1147
1148 gretl_matrix_free(X);
1149 gretl_matrix_free(diag);
1150 gretl_matrix_free(tmp1);
1151 gretl_matrix_free(tmp2);
1152 gretl_matrix_free(V);
1153 gretl_matrix_free(u);
1154
1155 return err;
1156 }
1157
1158 /* Variant of qr_make_hccme for use when a model has
1159 been estimated via matrix methods. On input the
1160 vector @d should hold the squared residuals, and
1161 @h the diagonal elements of the "hat" matrix.
1162 On return @d will be overwritten (used as work-
1163 space), and if successful @V will contain the
1164 requested HCCME variant.
1165 */
1166
qr_matrix_hccme(const gretl_matrix * X,const gretl_matrix * h,const gretl_matrix * XTXi,gretl_matrix * d,gretl_matrix * VCV,int hc_version)1167 int qr_matrix_hccme (const gretl_matrix *X,
1168 const gretl_matrix *h,
1169 const gretl_matrix *XTXi,
1170 gretl_matrix *d,
1171 gretl_matrix *VCV,
1172 int hc_version)
1173 {
1174 gretl_matrix *tmp1 = NULL;
1175 gretl_matrix *tmp2 = NULL;
1176 int T = X->rows;
1177 int k = X->cols;
1178 int t, err = 0;
1179
1180 tmp1 = gretl_matrix_alloc(k, T);
1181 tmp2 = gretl_matrix_alloc(k, k);
1182
1183 if (tmp1 == NULL || tmp2 == NULL) {
1184 gretl_matrix_free(tmp1);
1185 gretl_matrix_free(tmp2);
1186 return E_ALLOC;
1187 }
1188
1189 if (hc_version == 1) {
1190 for (t=0; t<T; t++) {
1191 d->val[t] *= (double) T / (T - k);
1192 }
1193 } else if (hc_version > 1) {
1194 /* do the h_t calculations */
1195 for (t=0; t<T; t++) {
1196 double ht = h->val[t];
1197
1198 if (hc_version == 2) {
1199 d->val[t] /= (1.0 - ht);
1200 } else {
1201 /* HC3 */
1202 d->val[t] /= (1.0 - ht) * (1.0 - ht);
1203 }
1204 }
1205 }
1206
1207 do_X_prime_diag(X, d, tmp1);
1208 gretl_matrix_multiply(tmp1, X, tmp2);
1209 gretl_matrix_qform(XTXi, GRETL_MOD_NONE, tmp2,
1210 VCV, GRETL_MOD_NONE);
1211
1212 gretl_matrix_free(tmp1);
1213 gretl_matrix_free(tmp2);
1214
1215 return err;
1216 }
1217
qr_dw_stats(MODEL * pmod,const DATASET * dset,gretl_matrix * X,gretl_matrix * u)1218 static int qr_dw_stats (MODEL *pmod, const DATASET *dset,
1219 gretl_matrix *X, gretl_matrix *u)
1220 {
1221 double DW, pv;
1222 int t, s, err = 0;
1223
1224 get_data_X(X, pmod, dset);
1225
1226 for (s=0, t=pmod->t1; t<=pmod->t2; s++, t++) {
1227 gretl_vector_set(u, s, pmod->uhat[t]);
1228 }
1229
1230 pv = dw_pval(u, X, &DW, &err);
1231
1232 if (!err) {
1233 pmod->dw = DW;
1234 gretl_model_set_double(pmod, "dw_pval", pv);
1235 }
1236
1237 return err;
1238 }
1239
qr_make_regular_vcv(MODEL * pmod,gretl_matrix * v,gretlopt opt)1240 static int qr_make_regular_vcv (MODEL *pmod, gretl_matrix *v,
1241 gretlopt opt)
1242 {
1243 int flag = (opt & OPT_X)? VCV_XPX : VCV_SIMPLE;
1244
1245 return qr_make_vcv(pmod, v, flag);
1246 }
1247
get_model_data(MODEL * pmod,const DATASET * dset,gretl_matrix * Q,gretl_matrix * y)1248 static void get_model_data (MODEL *pmod, const DATASET *dset,
1249 gretl_matrix *Q, gretl_matrix *y)
1250 {
1251 double rho = gls_rho(pmod);
1252 int qdiff = (rho != 0.0);
1253 int pwe = (pmod->opt & OPT_P);
1254 double xt, pw1 = 0.0;
1255 int i, s, t;
1256
1257 if (pmod->missmask == NULL && !pwe && !qdiff && !pmod->nwt) {
1258 /* simple case: no missing values and no data transformation
1259 called for, so use faster procedure
1260 */
1261 int T = pmod->t2 - pmod->t1 + 1;
1262 size_t sz = T * sizeof(double);
1263
1264 s = 0;
1265 for (i=2; i<=pmod->list[0]; i++) {
1266 memcpy(Q->val + s, dset->Z[pmod->list[i]] + pmod->t1, sz);
1267 s += T;
1268 }
1269 if (y != NULL) {
1270 memcpy(y->val, dset->Z[pmod->list[1]] + pmod->t1, sz);
1271 }
1272 return;
1273 }
1274
1275 if (pwe) {
1276 pw1 = sqrt(1.0 - rho * rho);
1277 }
1278
1279 /* copy independent vars into matrix Q */
1280 s = 0;
1281 for (i=2; i<=pmod->list[0]; i++) {
1282 int vi = pmod->list[i];
1283
1284 for (t=pmod->t1; t<=pmod->t2; t++) {
1285 if (model_missing(pmod, t)) {
1286 continue;
1287 }
1288 xt = dset->Z[vi][t];
1289 if (pmod->nwt) {
1290 if (dset->Z[pmod->nwt][t] == 0.0) {
1291 continue;
1292 } else {
1293 xt *= sqrt(dset->Z[pmod->nwt][t]);
1294 }
1295 } else if (qdiff) {
1296 if (pwe && t == pmod->t1) {
1297 xt *= pw1;
1298 } else {
1299 xt -= rho * dset->Z[vi][t-1];
1300 }
1301 }
1302 Q->val[s++] = xt;
1303 }
1304 }
1305
1306 if (y != NULL) {
1307 /* copy dependent variable into y vector */
1308 int vy = pmod->list[1];
1309
1310 s = 0;
1311 for (t=pmod->t1; t<=pmod->t2; t++) {
1312 if (model_missing(pmod, t)) {
1313 continue;
1314 }
1315 xt = dset->Z[vy][t];
1316 if (pmod->nwt) {
1317 if (dset->Z[pmod->nwt][t] == 0.0) {
1318 continue;
1319 } else {
1320 xt *= sqrt(dset->Z[pmod->nwt][t]);
1321 }
1322 } else if (qdiff) {
1323 if (pwe && t == pmod->t1) {
1324 xt *= pw1;
1325 } else {
1326 xt -= rho * dset->Z[vy][t-1];
1327 }
1328 }
1329 y->val[s++] = xt;
1330 }
1331 }
1332 }
1333
1334 static int
allocate_model_arrays(MODEL * pmod,int k,int T)1335 allocate_model_arrays (MODEL *pmod, int k, int T)
1336 {
1337 if (pmod->sderr == NULL) {
1338 pmod->sderr = malloc(k * sizeof *pmod->sderr);
1339 }
1340
1341 if (pmod->yhat == NULL) {
1342 pmod->yhat = malloc(T * sizeof *pmod->yhat);
1343 }
1344
1345 if (pmod->uhat == NULL) {
1346 pmod->uhat = malloc(T * sizeof *pmod->uhat);
1347 }
1348
1349 if (pmod->sderr == NULL || pmod->yhat == NULL ||
1350 pmod->uhat == NULL) {
1351 return E_ALLOC;
1352 }
1353
1354 return 0;
1355 }
1356
1357 #define RCOND_WARN 1.0e-07
1358
1359 /* perform QR decomposition plus some additional tasks */
1360
QR_decomp_plus(gretl_matrix * Q,gretl_matrix * R,int * rank,int * warn)1361 static int QR_decomp_plus (gretl_matrix *Q, gretl_matrix *R,
1362 int *rank, int *warn)
1363 {
1364 integer k = gretl_matrix_rows(R);
1365 double rcond = 0;
1366 int r, err;
1367
1368 if (warn != NULL) {
1369 *warn = 0;
1370 }
1371
1372 /* basic decomposition */
1373 err = gretl_matrix_QR_decomp(Q, R);
1374 if (err) {
1375 return err;
1376 }
1377
1378 /* check rank of QR */
1379 r = gretl_check_QR_rank(R, &err, &rcond);
1380 if (err) {
1381 return err;
1382 }
1383
1384 if (r < k) {
1385 err = E_SINGULAR;
1386 } else {
1387 /* then invert the triangular R */
1388 char uplo = 'U';
1389 char diag = 'N';
1390 integer info = 0;
1391
1392 dtrtri_(&uplo, &diag, &k, R->val, &k, &info);
1393 if (info != 0) {
1394 fprintf(stderr, "dtrtri: info = %d\n", (int) info);
1395 err = 1;
1396 } else if (0 && rcond < RCOND_WARN && warn != NULL) {
1397 /* 2020-02-19: not sure we want this! */
1398 fprintf(stderr, "QR_decomp_plus: rcond = %g\n", rcond);
1399 *warn = 1;
1400 }
1401 }
1402
1403 if (rank != NULL) {
1404 *rank = r;
1405 }
1406
1407 return err;
1408 }
1409
1410 #define REDEBUG 0
1411
1412 static void
drop_redundant_vars(MODEL * pmod,DATASET * dset,gretl_matrix * R,int rank,gretlopt opt)1413 drop_redundant_vars (MODEL *pmod, DATASET *dset, gretl_matrix *R,
1414 int rank, gretlopt opt)
1415 {
1416 int *droplist = NULL;
1417 int i, vi, pos, nd;
1418 double d;
1419
1420 #if REDEBUG
1421 printlist(pmod->list, "pmod->list, into drop_redundant_vars");
1422 fprintf(stderr, "rank = %d\n", rank);
1423 #endif
1424
1425 pos = 2;
1426 nd = 0;
1427 for (i=0; i<R->rows; i++) {
1428 d = gretl_matrix_get(R, i, i);
1429 if (fabs(d) < R_DIAG_MIN) {
1430 vi = pmod->list[pos];
1431 gretl_list_append_term(&droplist, vi);
1432 fprintf(stderr, "dropping redundant variable %d (%s): d = %g\n",
1433 vi, dset->varname[vi], d);
1434 gretl_list_delete_at_pos(pmod->list, pos--);
1435 nd++;
1436 }
1437 pos++;
1438 }
1439
1440 pmod->ncoeff -= nd;
1441 pmod->dfd = pmod->nobs - pmod->ncoeff;
1442 pmod->dfn = pmod->ncoeff - pmod->ifc;
1443
1444 if (droplist != NULL) {
1445 gretl_model_set_list_as_data(pmod, "droplist", droplist);
1446 }
1447 }
1448
gretl_qr_regress(MODEL * pmod,DATASET * dset,gretlopt opt)1449 int gretl_qr_regress (MODEL *pmod, DATASET *dset, gretlopt opt)
1450 {
1451 integer T, k;
1452 gretl_matrix *y = NULL;
1453 gretl_matrix *Q = NULL;
1454 gretl_matrix *R = NULL;
1455 gretl_matrix *g = NULL;
1456 gretl_matrix *b = NULL;
1457 gretl_matrix *V = NULL;
1458 int rank, warn = 0, err = 0;
1459
1460 T = pmod->nobs; /* # of rows (observations) */
1461 k = pmod->list[0] - 1; /* # of cols (variables) */
1462
1463 y = gretl_matrix_alloc(T, 1);
1464 Q = gretl_matrix_alloc(T, k);
1465 R = gretl_matrix_alloc(k, k);
1466 V = gretl_matrix_alloc(k, k);
1467
1468 if (y == NULL || Q == NULL || R == NULL || V == NULL) {
1469 err = E_ALLOC;
1470 goto qr_cleanup;
1471 }
1472
1473 get_model_data(pmod, dset, Q, y);
1474 err = QR_decomp_plus(Q, R, &rank, &warn);
1475
1476 /* handling of near-perfect collinearity */
1477 if (err == E_SINGULAR && !(opt & OPT_Z)) {
1478 drop_redundant_vars(pmod, dset, R, rank, opt);
1479 k = pmod->list[0] - 1;
1480 gretl_matrix_reuse(Q, T, k);
1481 gretl_matrix_reuse(R, k, k);
1482 gretl_matrix_reuse(V, k, k);
1483 get_model_data(pmod, dset, Q, y);
1484 err = QR_decomp_plus(Q, R, NULL, &warn);
1485 if (!err) {
1486 maybe_shift_ldepvar(pmod, dset);
1487 }
1488 }
1489
1490 if (err) {
1491 goto qr_cleanup;
1492 }
1493
1494 /* allocate temporary arrays */
1495 g = gretl_matrix_alloc(k, 1);
1496 b = gretl_matrix_alloc(k, 1);
1497
1498 if (g == NULL || b == NULL) {
1499 err = E_ALLOC;
1500 } else {
1501 err = allocate_model_arrays(pmod, k, dset->n);
1502 }
1503
1504 if (err) {
1505 goto qr_cleanup;
1506 }
1507
1508 /* make "g" into gamma-hat */
1509 gretl_matrix_multiply_mod(Q, GRETL_MOD_TRANSPOSE,
1510 y, GRETL_MOD_NONE,
1511 g, GRETL_MOD_NONE);
1512
1513 /* OLS coefficients */
1514 gretl_matrix_multiply(R, g, b);
1515 pmod->coeff = gretl_matrix_steal_data(b);
1516
1517 /* write vector of fitted values into y */
1518 gretl_matrix_multiply(Q, g, y);
1519
1520 /* get vector of residuals and SSR */
1521 get_resids_and_SSR(pmod, dset, y, dset->n);
1522
1523 /* standard error of regression */
1524 if (T - k > 0) {
1525 if (pmod->opt & OPT_N) {
1526 /* no-df-corr */
1527 pmod->sigma = sqrt(pmod->ess / T);
1528 } else {
1529 pmod->sigma = sqrt(pmod->ess / (T - k));
1530 }
1531 } else {
1532 pmod->sigma = 0.0;
1533 }
1534
1535 /* create (X'X)^{-1} */
1536 gretl_matrix_multiply_mod(R, GRETL_MOD_NONE,
1537 R, GRETL_MOD_TRANSPOSE,
1538 V, GRETL_MOD_NONE);
1539
1540 /* VCV and standard errors */
1541 if (opt & OPT_R) {
1542 pmod->opt |= OPT_R;
1543 if (opt & OPT_C) {
1544 err = qr_make_cluster_vcv(pmod, pmod->ci, dset, V, opt);
1545 } else if ((opt & OPT_T) && !libset_get_bool(FORCE_HC)) {
1546 err = qr_make_hac(pmod, dset, V);
1547 } else {
1548 err = qr_make_hccme(pmod, dset, Q, V);
1549 }
1550 } else {
1551 err = qr_make_regular_vcv(pmod, V, opt);
1552 }
1553
1554 if ((opt & OPT_R) && err == E_JACOBIAN) {
1555 /* try fallback? */
1556 err = qr_make_regular_vcv(pmod, V, opt);
1557 if (!err) {
1558 gretl_model_set_int(pmod, "non-robust", 1);
1559 }
1560 }
1561
1562 if (!err) {
1563 /* get R^2, F-stat */
1564 qr_compute_stats(pmod, dset, T, opt);
1565
1566 /* D-W stat and p-value */
1567 if ((opt & OPT_I) && pmod->missmask == NULL) {
1568 qr_dw_stats(pmod, dset, Q, y);
1569 }
1570
1571 /* near-singularity? */
1572 if (warn) {
1573 gretl_model_set_int(pmod, "near-singular", 1);
1574 }
1575 }
1576
1577 qr_cleanup:
1578
1579 gretl_matrix_free(Q);
1580 gretl_matrix_free(R);
1581 gretl_matrix_free(y);
1582
1583 gretl_matrix_free(g);
1584 gretl_matrix_free(b);
1585 gretl_matrix_free(V);
1586
1587 pmod->errcode = err;
1588
1589 return err;
1590 }
1591
lapack_cholesky_regress(MODEL * pmod,const DATASET * dset,gretlopt opt)1592 int lapack_cholesky_regress (MODEL *pmod, const DATASET *dset,
1593 gretlopt opt)
1594 {
1595 integer T, k;
1596 gretl_matrix *y = NULL;
1597 gretl_matrix *X = NULL;
1598 gretl_matrix *b = NULL;
1599 gretl_matrix *XTX = NULL;
1600 int err = 0;
1601
1602 T = pmod->nobs; /* # of rows (observations) */
1603 k = pmod->list[0] - 1; /* # of cols (variables) */
1604
1605 y = gretl_matrix_alloc(T, 1);
1606 X = gretl_matrix_alloc(T, k);
1607 b = gretl_matrix_alloc(k, 1);
1608
1609 if (y == NULL || X == NULL || b == NULL) {
1610 err = E_ALLOC;
1611 goto ch_cleanup;
1612 }
1613
1614 get_model_data(pmod, dset, X, y);
1615
1616 XTX = gretl_matrix_XTX_new(X);
1617 if (XTX == NULL) {
1618 err = E_ALLOC;
1619 }
1620 if (!err) {
1621 err = gretl_matrix_multiply_mod(X, GRETL_MOD_TRANSPOSE,
1622 y, GRETL_MOD_NONE,
1623 b, GRETL_MOD_NONE);
1624 }
1625 if (!err) {
1626 err = gretl_cholesky_decomp_solve(XTX, b);
1627 }
1628
1629 if (!err) {
1630 err = allocate_model_arrays(pmod, k, dset->n);
1631 }
1632
1633 if (err) {
1634 goto ch_cleanup;
1635 }
1636
1637 /* write vector of fitted values into y */
1638 gretl_matrix_multiply(X, b, y);
1639
1640 /* OLS coefficients */
1641 pmod->coeff = gretl_matrix_steal_data(b);
1642
1643 /* get vector of residuals and SSR */
1644 get_resids_and_SSR(pmod, dset, y, dset->n);
1645
1646 /* standard error of regression */
1647 if (T - k > 0) {
1648 if (pmod->opt & OPT_N) {
1649 /* no-df-corr */
1650 pmod->sigma = sqrt(pmod->ess / T);
1651 } else {
1652 pmod->sigma = sqrt(pmod->ess / (T - k));
1653 }
1654 } else {
1655 pmod->sigma = 0.0;
1656 }
1657
1658 /* create (X'X)^{-1} */
1659 err = gretl_cholesky_invert(XTX);
1660 if (err) {
1661 goto ch_cleanup;
1662 }
1663
1664 /* VCV and standard errors */
1665 if (opt & OPT_R) {
1666 pmod->opt |= OPT_R;
1667 if (opt & OPT_C) {
1668 err = qr_make_cluster_vcv(pmod, OLS, dset, XTX, opt);
1669 } else if ((opt & OPT_T) && !libset_get_bool(FORCE_HC)) {
1670 err = qr_make_hac(pmod, dset, XTX);
1671 } else {
1672 err = qr_make_hccme(pmod, dset, X, XTX);
1673 }
1674 } else {
1675 err = qr_make_regular_vcv(pmod, XTX, opt);
1676 }
1677
1678 if ((opt & OPT_R) && err == E_JACOBIAN) {
1679 /* try fallback? */
1680 err = qr_make_regular_vcv(pmod, XTX, opt);
1681 if (!err) {
1682 gretl_model_set_int(pmod, "non-robust", 1);
1683 }
1684 }
1685
1686 if (!err) {
1687 /* get R^2, F-stat */
1688 qr_compute_stats(pmod, dset, T, opt);
1689
1690 /* D-W stat and p-value */
1691 if ((opt & OPT_I) && pmod->missmask == NULL) {
1692 qr_dw_stats(pmod, dset, X, y);
1693 }
1694 }
1695
1696 ch_cleanup:
1697
1698 gretl_matrix_free(X);
1699 gretl_matrix_free(y);
1700 gretl_matrix_free(b);
1701 gretl_matrix_free(XTX);
1702
1703 pmod->errcode = (err == E_NOTPD)? E_SINGULAR: err;
1704
1705 return err;
1706 }
1707
qr_tsls_vcv(MODEL * pmod,const DATASET * dset,gretlopt opt)1708 int qr_tsls_vcv (MODEL *pmod, const DATASET *dset, gretlopt opt)
1709 {
1710 gretl_matrix *Q = NULL;
1711 gretl_matrix *R = NULL;
1712 gretl_matrix *V = NULL;
1713 int k, err = 0;
1714
1715 k = pmod->list[0] - 1;
1716
1717 Q = make_data_X(pmod, dset);
1718 R = gretl_matrix_alloc(k, k);
1719 V = gretl_matrix_alloc(k, k);
1720
1721 if (Q == NULL || R == NULL || V == NULL) {
1722 err = E_ALLOC;
1723 goto qr_cleanup;
1724 }
1725
1726 err = QR_decomp_plus(Q, R, NULL, NULL);
1727 if (err) {
1728 goto qr_cleanup;
1729 }
1730
1731 /* create (X'X)^{-1} */
1732 gretl_matrix_multiply_mod(R, GRETL_MOD_NONE,
1733 R, GRETL_MOD_TRANSPOSE,
1734 V, GRETL_MOD_NONE);
1735
1736 /* VCV and standard errors */
1737 if (opt & OPT_R) {
1738 if (opt & OPT_C) {
1739 pmod->opt |= OPT_R;
1740 err = qr_make_cluster_vcv(pmod, IVREG, dset, V, opt);
1741 } else if (dataset_is_panel(dset)) {
1742 err = qr_make_regular_vcv(pmod, V, OPT_X);
1743 if (!err) {
1744 pmod->ci = IVREG; /* ? */
1745 err = panel_tsls_robust_vcv(pmod, dset);
1746 }
1747 } else if (dataset_is_time_series(dset) &&
1748 !libset_get_bool(FORCE_HC)) {
1749 pmod->opt |= OPT_R;
1750 err = qr_make_hac(pmod, dset, V);
1751 } else {
1752 pmod->opt |= OPT_R;
1753 err = qr_make_hccme(pmod, dset, Q, V);
1754 }
1755 } else {
1756 qr_make_regular_vcv(pmod, V, OPT_NONE);
1757 }
1758
1759 qr_cleanup:
1760
1761 gretl_matrix_free(Q);
1762 gretl_matrix_free(R);
1763 gretl_matrix_free(V);
1764
1765 pmod->errcode = err;
1766
1767 return err;
1768 }
1769
cval_count(MODEL * pmod,double cvi,const double * cZ)1770 static int cval_count (MODEL *pmod, double cvi, const double *cZ)
1771 {
1772 int t, cc = 0;
1773
1774 for (t=pmod->t1; t<=pmod->t2; t++) {
1775 if (!na(pmod->uhat[t]) && cZ[t] == cvi) {
1776 cc++;
1777 }
1778 }
1779
1780 return cc;
1781 }
1782
cval_count_max(MODEL * pmod,const gretl_matrix * cvals,const double * cZ)1783 static int cval_count_max (MODEL *pmod, const gretl_matrix *cvals,
1784 const double *cZ)
1785 {
1786 int n = gretl_vector_get_length(cvals);
1787 int i, cc, cmax = 0;
1788
1789 for (i=0; i<n; i++) {
1790 cc = cval_count(pmod, cvals->val[i], cZ);
1791 if (cc > cmax) {
1792 cmax = cc;
1793 }
1794 }
1795
1796 return cmax;
1797 }
1798
1799 #define CDEBUG 0
1800
cluster_vcv_calc(MODEL * pmod,int cvar,gretl_matrix * cvals,gretl_matrix * XX,const DATASET * dset,int * err)1801 static gretl_matrix *cluster_vcv_calc (MODEL *pmod,
1802 int cvar,
1803 gretl_matrix *cvals,
1804 gretl_matrix *XX,
1805 const DATASET *dset,
1806 int *err)
1807
1808 {
1809 gretl_matrix *V = NULL;
1810 gretl_matrix *W = NULL;
1811 gretl_matrix *XXW = NULL;
1812 gretl_vector *ei = NULL;
1813 gretl_matrix *Xi = NULL;
1814 gretl_vector *eXi = NULL;
1815 const double *cZ;
1816 const double *wgt = NULL;
1817 int n_c, M, N, k = pmod->ncoeff;
1818 int total_obs = 0;
1819 int i, j, v, t;
1820
1821 cZ = dset->Z[cvar];
1822 N = cval_count_max(pmod, cvals, cZ);
1823 #if CDEBUG
1824 fprintf(stderr, "max cval count = %d\n", N);
1825 #endif
1826
1827 V = gretl_matrix_alloc(k, k);
1828 W = gretl_zero_matrix_new(k, k);
1829 XXW = gretl_zero_matrix_new(k, k);
1830 ei = gretl_column_vector_alloc(N);
1831 Xi = gretl_matrix_alloc(N, k);
1832 eXi = gretl_vector_alloc(k);
1833
1834 if (V == NULL || W == NULL || XXW == NULL ||
1835 ei == NULL || Xi == NULL || eXi == NULL) {
1836 *err = E_ALLOC;
1837 goto bailout;
1838 }
1839
1840 if (pmod->nwt) {
1841 wgt = dset->Z[pmod->nwt];
1842 }
1843
1844 M = gretl_vector_get_length(cvals);
1845 n_c = 0;
1846
1847 for (i=0; i<M; i++) {
1848 double cvi = cvals->val[i];
1849 int Ni = cval_count(pmod, cvi, cZ);
1850 int s = 0;
1851
1852 if (Ni == 0) {
1853 continue;
1854 }
1855
1856 #if CDEBUG
1857 fprintf(stderr, "i=%d, cvi=%g, Ni=%d\n", i, cvi, Ni);
1858 #endif
1859 ei = gretl_matrix_reuse(ei, Ni, -1);
1860 Xi = gretl_matrix_reuse(Xi, Ni, -1);
1861
1862 for (t=pmod->t1; t<=pmod->t2; t++) {
1863 if (!na(pmod->uhat[t]) && cZ[t] == cvi) {
1864 gretl_vector_set(ei, s, pmod->uhat[t]);
1865 for (j=0; j<k; j++) {
1866 v = pmod->list[j+2];
1867 if (wgt != NULL) {
1868 gretl_matrix_set(Xi, s, j, dset->Z[v][t] * sqrt(wgt[t]));
1869 } else {
1870 gretl_matrix_set(Xi, s, j, dset->Z[v][t]);
1871 }
1872 }
1873 s++;
1874 }
1875 if (s == Ni) {
1876 /* we've filled this matrix */
1877 break;
1878 }
1879 }
1880
1881 gretl_matrix_multiply_mod(ei, GRETL_MOD_TRANSPOSE,
1882 Xi, GRETL_MOD_NONE,
1883 eXi, GRETL_MOD_NONE);
1884 gretl_matrix_multiply_mod(eXi, GRETL_MOD_TRANSPOSE,
1885 eXi, GRETL_MOD_NONE,
1886 W, GRETL_MOD_CUMULATE);
1887 #if CDEBUG > 1
1888 gretl_matrix_print(ei, "e(i)");
1889 gretl_matrix_print(Xi, "X(i)");
1890 gretl_matrix_print(W, "W");
1891 #endif
1892 n_c++;
1893 total_obs += s;
1894 }
1895
1896 if (n_c < 2) {
1897 gretl_errmsg_set("Invalid clustering variable");
1898 *err = E_DATA;
1899 goto bailout;
1900 } else if (total_obs < pmod->nobs) {
1901 *err = E_MISSDATA;
1902 goto bailout;
1903 }
1904
1905 /* form V(W) = (X'X)^{-1} W (X'X)^{-1} */
1906 gretl_matrix_multiply(XX, W, XXW);
1907 gretl_matrix_multiply(XXW, XX, V);
1908 gretl_matrix_xtr_symmetric(V);
1909
1910 #if CDEBUG
1911 gretl_matrix_print(XX, "X'X^{-1}");
1912 gretl_matrix_print(W, "W");
1913 gretl_matrix_print(V, "V");
1914 #endif
1915
1916 if (!(pmod->opt & OPT_N)) {
1917 /* apply df adjustment a la Stata */
1918 /* FIXME IVREG case? */
1919 double dfadj;
1920
1921 N = pmod->nobs;
1922 dfadj = (M/(M-1.0)) * (N-1.0)/(N-k);
1923 gretl_matrix_multiply_by_scalar(V, dfadj);
1924 #if CDEBUG > 1
1925 gretl_matrix_print(V, "V(adjusted)");
1926 #endif
1927
1928 }
1929
1930 bailout:
1931
1932 gretl_matrix_free(W);
1933 gretl_matrix_free(XXW);
1934 gretl_matrix_free(ei);
1935 gretl_matrix_free(Xi);
1936 gretl_matrix_free(eXi);
1937
1938 if (*err) {
1939 gretl_matrix_free(V);
1940 V = NULL;
1941 }
1942
1943 return V;
1944 }
1945
1946 /* Get the sorted values of the clustering series, @cvar, checking
1947 for missing values as we go. This is a little more complicated
1948 if there are interior missing values for the regressand or
1949 regressors: in that case we need to construct a temporary
1950 array to hold the relevant values of the clustering series.
1951 */
1952
cluster_var_values(const double * cvar,MODEL * pmod,int * err)1953 static gretl_matrix *cluster_var_values (const double *cvar,
1954 MODEL *pmod,
1955 int *err)
1956 {
1957 gretl_matrix *cvals = NULL;
1958 int t;
1959
1960 if (pmod->missmask != NULL) {
1961 double *ctmp = malloc(pmod->nobs * sizeof *ctmp);
1962
1963 if (ctmp == NULL) {
1964 *err = E_ALLOC;
1965 } else {
1966 int i = 0;
1967
1968 for (t=pmod->t1; t<=pmod->t2 && !*err; t++) {
1969 if (pmod->missmask[t] != '1') {
1970 if (na(cvar[t])) {
1971 *err = E_MISSDATA;
1972 } else {
1973 ctmp[i++] = cvar[t];
1974 }
1975 }
1976 }
1977
1978 if (!*err) {
1979 cvals = gretl_matrix_values(ctmp, pmod->nobs, OPT_S, err);
1980 }
1981 free(ctmp);
1982 }
1983 } else {
1984 /* no interior missing values for y or X */
1985 for (t=pmod->t1; t<=pmod->t2 && !*err; t++) {
1986 if (na(cvar[t])) {
1987 *err = E_MISSDATA;
1988 }
1989 }
1990 if (!*err) {
1991 cvals = gretl_matrix_values(cvar + pmod->t1, pmod->nobs,
1992 OPT_S, err);
1993 }
1994 }
1995
1996 return cvals;
1997 }
1998
1999 static int cluster_vcv_ci;
2000
2001 /**
2002 * qr_make_cluster_vcv:
2003 * @pmod: pointer to model.
2004 * @ci: command index (right now, OLS or IVREG).
2005 * @dset: pointer to dataset.
2006 * @XX: X'X matrix.
2007 *
2008 * Compute and set on @pmod a variance matrix that is "clustered"
2009 * by the value of a selected variable via the --cluster=foo
2010 * command-line option.
2011 *
2012 * Returns: 0 on success, non-zero code on error.
2013 */
2014
qr_make_cluster_vcv(MODEL * pmod,int ci,const DATASET * dset,gretl_matrix * XX,gretlopt opt)2015 static int qr_make_cluster_vcv (MODEL *pmod, int ci,
2016 const DATASET *dset,
2017 gretl_matrix *XX,
2018 gretlopt opt)
2019 {
2020 gretl_matrix *cvals = NULL;
2021 gretl_matrix *V = NULL;
2022 const char *cname;
2023 int cvar, n_c = 0;
2024 int err = 0;
2025
2026 if (pmod->ci != OLS && pmod->ci != IVREG && pmod->ci != WLS) {
2027 /* relax this? */
2028 return E_NOTIMP;
2029 }
2030
2031 if (cluster_vcv_ci != 0) {
2032 ci = cluster_vcv_ci;
2033 cluster_vcv_ci = 0;
2034 }
2035
2036 cname = get_optval_string(ci, OPT_C);
2037 if (cname == NULL) {
2038 gretl_errmsg_set("Got --cluster option but couldn't find varname");
2039 return E_PARSE;
2040 }
2041
2042 cvar = current_series_index(dset, cname);
2043 if (cvar < 1 || cvar >= dset->v) {
2044 err = E_UNKVAR;
2045 }
2046
2047 if (!err) {
2048 cvals = cluster_var_values(dset->Z[cvar], pmod, &err);
2049 if (!err) {
2050 n_c = gretl_vector_get_length(cvals);
2051 if (n_c < 2) {
2052 err = E_DATA;
2053 }
2054 }
2055 }
2056
2057 #if CDEBUG
2058 fprintf(stderr, "qr_make_cluster_vcv: err = %d\n", err);
2059 fprintf(stderr, "cluster var = %s (%d)\n", cname, cvar);
2060 gretl_matrix_print(cvals, "cvals");
2061 #endif
2062
2063 if (!err) {
2064 V = cluster_vcv_calc(pmod, cvar, cvals, XX, dset, &err);
2065 }
2066
2067 if (!err) {
2068 err = gretl_model_write_vcv(pmod, V);
2069 }
2070
2071 if (!err) {
2072 gretl_model_set_vcv_info(pmod, VCV_CLUSTER, cvar);
2073 gretl_model_set_int(pmod, "n_clusters", n_c);
2074 }
2075
2076 gretl_matrix_free(V);
2077 gretl_matrix_free(cvals);
2078
2079 return err;
2080 }
2081
set_cluster_vcv_ci(int ci)2082 int set_cluster_vcv_ci (int ci)
2083 {
2084 if (ci != 0 && !valid_short_opt(ci, 'c')) {
2085 return E_BADOPT;
2086 } else {
2087 cluster_vcv_ci = ci;
2088 return 0;
2089 }
2090 }
2091