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 "var.h"
22 #include "johansen.h"
23 #include "vartest.h"
24 #include "matrix_extra.h"
25 #include "libset.h"
26 #include "qr_estimate.h"
27 #include "varprint.h"
28
gretl_VAR_normality_test(const GRETL_VAR * var,gretlopt opt,PRN * prn)29 int gretl_VAR_normality_test (const GRETL_VAR *var,
30 gretlopt opt, PRN *prn)
31 {
32 int err = 0;
33
34 if (var->E == NULL || var->S == NULL) {
35 err = 1;
36 } else {
37 err = multivariate_normality_test(var->E, var->S,
38 opt, prn);
39 }
40
41 return err;
42 }
43
44 /* Multivariate version of Breusch-Godfrey test, as per
45 H. Lutkepohl, New Introduction to Multiple Time Series
46 Analysis (Springer, 2005) section 4.4, pp. 173-4.
47 */
48
49 static int
multivariate_autocorr_test(GRETL_VAR * var,int H,int autoH,gretlopt opt,PRN * prn)50 multivariate_autocorr_test (GRETL_VAR *var, int H,
51 int autoH, gretlopt opt,
52 PRN *prn)
53 {
54 const gretl_matrix *U;
55 char Fspec[32];
56 gretl_matrix *tests, *pvals;
57 gretl_matrix *B = NULL;
58 gretl_matrix *E = NULL;
59 gretl_matrix *X = NULL;
60 gretl_matrix *S = NULL;
61 double *targ;
62 double ulag, s, FRao;
63 double detUU, detEE;
64 double N, dfn, dfd;
65 int quiet = (opt & OPT_Q);
66 int h, h2, K = var->neqns;
67 int i, t, nx, K2 = K * K;
68 int g, lagcol, hw = 0;
69 int p = var->order;
70 int T = var->T;
71 int err = 0;
72
73 /* we need var->X and var->S below */
74 if (var == NULL || var->X == NULL || var->S == NULL) {
75 return E_DATA;
76 }
77
78 /* and also the matrix of residuals */
79 U = gretl_VAR_get_residual_matrix(var);
80 if (U == NULL) {
81 return E_DATA;
82 }
83
84 /* how many regressors are we going to have at max? */
85 g = var->ncoeff;
86 if (!var->ifc) {
87 /* we'll have to add an intercept */
88 g++;
89 }
90 nx = g + K * H;
91
92 if (nx >= T) {
93 /* not enough data to do this? */
94 if (autoH) {
95 H = floor((T - 1 - g) / (double) K);
96 if (H <= 0) {
97 return E_TOOFEW;
98 } else {
99 nx = g + K * H;
100 }
101 } else {
102 return E_TOOFEW;
103 }
104 }
105
106 tests = gretl_column_vector_alloc(H);
107 pvals = gretl_column_vector_alloc(H);
108
109 if (tests == NULL || pvals == NULL) {
110 return E_ALLOC;
111 }
112
113 /* calculate determinant of cross-equation Sigma */
114 S = gretl_matrix_copy(var->S);
115 if (S == NULL) {
116 err = E_ALLOC;
117 } else {
118 detUU = gretl_matrix_determinant(S, &err);
119 }
120
121 if (err) {
122 goto bailout;
123 }
124
125 /* allocate E, B, X to max size */
126 E = gretl_matrix_alloc(T, K);
127 B = gretl_matrix_alloc(nx, K);
128 X = gretl_matrix_alloc(T, nx);
129
130 if (E == NULL || B == NULL || X == NULL) {
131 err = E_ALLOC;
132 goto bailout;
133 }
134
135 if (var->ifc) {
136 targ = X->val;
137 } else {
138 /* insert constant */
139 for (t=0; t<T; t++) {
140 X->val[t] = 1.0;
141 }
142 targ = X->val + T;
143 }
144
145 /* copy the original regressors into X */
146 memcpy(targ, var->X->val, var->ncoeff * T * sizeof(double));
147
148 if (!quiet) {
149 pprintf(prn, "%s %d\n", _("Test for autocorrelation of order up to"), H);
150 hw = ceil(log10(H));
151 pputc(prn, '\n');
152 bufspace(9 + hw, prn);
153 pputs(prn, "Rao F Approx dist. p-value\n");
154 }
155
156 for (h=1; h<=H && !err; h++) {
157 nx = g + K * h;
158 gretl_matrix_reuse(B, nx, K);
159 gretl_matrix_reuse(X, T, nx);
160 h2 = h * h;
161 /* add next lag of U to the X matrix */
162 for (i=0; i<K; i++) {
163 lagcol = g + K*(h-1) + i;
164 for (t=0; t<T; t++) {
165 ulag = (t - h < 0)? 0.0 : gretl_matrix_get(U, t-h, i);
166 gretl_matrix_set(X, t, lagcol, ulag);
167 }
168 }
169 err = gretl_matrix_multi_SVD_ols(U, X, B, E, NULL);
170 if (!err) {
171 err = gretl_matrix_multiply_mod(E, GRETL_MOD_TRANSPOSE,
172 E, GRETL_MOD_NONE,
173 S, GRETL_MOD_NONE);
174 }
175 if (err) {
176 break;
177 }
178 gretl_matrix_divide_by_scalar(S, T);
179 /* calculate the test statistic */
180 s = sqrt((pow(K, 4) * h2 - 4.0) / (K2 + K2 * h2 - 5.0));
181 N = T - K*p - 1 - K*h - (K - K*h + 1) / 2.0;
182 dfn = K2 * h;
183 dfd = floor(N*s - 0.5*(K2 * h) + 1);
184 detEE = gretl_matrix_determinant(S, &err);
185 if (!err) {
186 FRao = pow(detUU / detEE, 1/s) - 1.0;
187 FRao *= dfd / dfn;
188 tests->val[h-1] = FRao;
189 pvals->val[h-1] = snedecor_cdf_comp(dfn, dfd, FRao);
190 if (!quiet) {
191 sprintf(Fspec, "F(%g, %g)", dfn, dfd);
192 pprintf(prn, "lag %*d %9.3f %-13s %5.4f\n", hw, h,
193 FRao, Fspec, pvals->val[h-1]);
194 }
195 }
196 }
197
198 if (!quiet) {
199 pputc(prn, '\n');
200 }
201
202 bailout:
203
204 gretl_matrix_free(E);
205 gretl_matrix_free(B);
206 gretl_matrix_free(X);
207 gretl_matrix_free(S);
208
209 if (!err) {
210 record_matrix_test_result(tests, pvals);
211 } else {
212 gretl_matrix_free(tests);
213 gretl_matrix_free(pvals);
214 }
215
216 return err;
217 }
218
219 static int
univariate_autocorr_test(GRETL_VAR * var,int order,DATASET * dset,gretlopt opt,PRN * prn)220 univariate_autocorr_test (GRETL_VAR *var, int order,
221 DATASET *dset, gretlopt opt,
222 PRN *prn)
223 {
224 MODEL *pmod;
225 gretl_matrix *tests, *pvals;
226 double lb;
227 int i, quiet = (opt & OPT_Q);
228 int err = 0;
229
230 tests = gretl_column_vector_alloc(var->neqns);
231 pvals = gretl_column_vector_alloc(var->neqns);
232
233 if (tests == NULL || pvals == NULL) {
234 err = E_ALLOC;
235 }
236
237 if (!quiet) {
238 pprintf(prn, "%s %d\n\n", _("Test for autocorrelation of order"),
239 order);
240 }
241
242 for (i=0; i<var->neqns && !err; i++) {
243 pmod = var->models[i];
244 if (!quiet) {
245 pprintf(prn, "%s %d:\n", _("Equation"), i + 1);
246 }
247 lb = ljung_box(order, pmod->t1, pmod->t2, pmod->uhat, &err);
248 if (!err) {
249 tests->val[i] = lb;
250 pvals->val[i] = chisq_cdf_comp(order, lb);
251 if (!(opt & OPT_Q)) {
252 pprintf(prn, "Ljung-Box Q' = %g %s = P(%s(%d) > %g) = %.3g\n",
253 lb, _("with p-value"), _("Chi-square"), order,
254 lb, pvals->val[i]);
255 pputc(prn, '\n');
256 }
257 }
258 }
259
260 if (!err) {
261 record_matrix_test_result(tests, pvals);
262 } else {
263 gretl_matrix_free(tests);
264 gretl_matrix_free(pvals);
265 }
266
267 return err;
268 }
269
gretl_VAR_autocorrelation_test(GRETL_VAR * var,int order,DATASET * dset,gretlopt opt,PRN * prn)270 int gretl_VAR_autocorrelation_test (GRETL_VAR *var, int order,
271 DATASET *dset, gretlopt opt,
272 PRN *prn)
273 {
274 int h = order;
275 int err;
276
277 if (order == 0) {
278 h = dset->pd;
279 }
280
281 if (opt & OPT_U) {
282 err = univariate_autocorr_test(var, h, dset, opt, prn);
283 } else {
284 int autoH = (order == 0);
285
286 err = multivariate_autocorr_test(var, h, autoH, opt, prn);
287 }
288
289 return err;
290 }
291
292 /* write the vech of @src into row @t of @targ */
293
vech_into_row(gretl_matrix * targ,int t,const gretl_matrix * src)294 static void vech_into_row (gretl_matrix *targ, int t,
295 const gretl_matrix *src)
296 {
297 double vij;
298 int j, i, k = 0;
299
300 for (j=0; j<src->cols; j++) {
301 for (i=0; i<=j; i++) {
302 vij = gretl_matrix_get(src, i, j);
303 gretl_matrix_set(targ, t, k++, vij);
304 }
305 }
306 }
307
308 /* write from @src into @targ, with a lag of @h,
309 starting at column @inicol and implicitly
310 discarding rows of @src that would go off the
311 end of @targ
312 */
313
inscribe_lag(gretl_matrix * targ,const gretl_matrix * src,int h,int inicol)314 static void inscribe_lag (gretl_matrix *targ,
315 const gretl_matrix *src,
316 int h, int inicol)
317 {
318 double xij;
319 int j, t, k = inicol;
320
321 for (j=0; j<src->cols; j++) {
322 for (t=h; t<targ->rows; t++) {
323 xij = gretl_matrix_get(src, t-h, j);
324 gretl_matrix_set(targ, t, k, xij);
325 }
326 k++;
327 }
328 }
329
multivariate_arch_test(GRETL_VAR * var,int H,int autoH,gretlopt opt,PRN * prn)330 static int multivariate_arch_test (GRETL_VAR *var, int H,
331 int autoH, gretlopt opt,
332 PRN *prn)
333 {
334 const gretl_matrix *U;
335 gretl_matrix_block *Bk;
336 gretl_matrix *tests, *pvals;
337 gretl_matrix *vU, *B, *X;
338 gretl_matrix *E, *S, *SS;
339 gretl_matrix *iC0, *ut, *uu;
340 double df, LM, K24;
341 int quiet = (opt & OPT_Q);
342 int h, K = var->neqns;
343 int i, t, nx, KK1;
344 int vdim, inicol;
345 int hw = 0;
346 int T = var->T;
347 int err = 0;
348
349 if (var == NULL) {
350 return E_DATA;
351 }
352
353 U = gretl_VAR_get_residual_matrix(var);
354 if (U == NULL) {
355 return E_DATA;
356 }
357
358 /* dimension of vech of variance */
359 vdim = (K * (K + 1)) / 2;
360
361 /* how many regressors are we going to have at max? */
362 nx = 1 + vdim * H;
363
364 if (nx >= T) {
365 /* not enough data to do this? */
366 if (autoH) {
367 H = floor((T - 1) / (double) vdim);
368 if (H <= 0) {
369 return E_TOOFEW;
370 } else {
371 nx = 1 + vdim * H;
372 }
373 } else {
374 return E_TOOFEW;
375 }
376 }
377
378 tests = gretl_column_vector_alloc(H);
379 pvals = gretl_column_vector_alloc(H);
380
381 if (tests == NULL || pvals == NULL) {
382 return E_ALLOC;
383 }
384
385 Bk = gretl_matrix_block_new(&vU, T, vdim,
386 &B, nx, vdim,
387 &E, T, vdim,
388 &X, T, nx,
389 &ut, 1, K,
390 &uu, K, K,
391 &S, vdim, vdim,
392 &SS, vdim, vdim,
393 &iC0, vdim, vdim,
394 NULL);
395 if (Bk == NULL) {
396 err = E_ALLOC;
397 goto bailout;
398 }
399
400 /* construct LHS (vU) */
401 for (t=0; t<T; t++) {
402 for (i=0; i<K; i++) {
403 ut->val[i] = gretl_matrix_get(U, t, i);
404 }
405 gretl_matrix_multiply_mod(ut, GRETL_MOD_TRANSPOSE,
406 ut, GRETL_MOD_NONE,
407 uu, GRETL_MOD_NONE);
408 vech_into_row(vU, t, uu);
409 }
410
411 /* initial setup of RHS (X) */
412 gretl_matrix_zero(X);
413 for (t=0; t<T; t++) {
414 X->val[t] = 1.0;
415 }
416
417 /* useful constants */
418 KK1 = K * (K+1);
419 K24 = KK1 * KK1 / 4.0;
420
421 if (!quiet) {
422 pprintf(prn, "%s %d\n", _("Test for ARCH of order up to"), H);
423 hw = ceil(log10(H));
424 pputc(prn, '\n');
425 bufspace(10 + hw, prn);
426 pputs(prn, "LM df p-value\n");
427 }
428
429 inicol = 1;
430
431 for (h=0; h<=H && !err; h++) {
432 nx = 1 + vdim * h;
433 gretl_matrix_reuse(B, nx, vdim);
434 gretl_matrix_reuse(X, T, nx);
435 if (h > 0) {
436 /* add next lag of vU to the X matrix */
437 inscribe_lag(X, vU, h, inicol);
438 inicol += vdim;
439 }
440 /* run aux regression */
441 err = gretl_matrix_multi_SVD_ols(vU, X, B, E, NULL);
442 if (!err) {
443 err = gretl_matrix_multiply_mod(E, GRETL_MOD_TRANSPOSE,
444 E, GRETL_MOD_NONE,
445 S, GRETL_MOD_NONE);
446 }
447 if (err) {
448 break;
449 }
450 gretl_matrix_divide_by_scalar(S, T);
451 if (h == 0) {
452 /* record baseline: inv(Cov0) */
453 gretl_matrix_copy_values(iC0, S);
454 err = gretl_invert_symmetric_matrix(iC0);
455 } else {
456 /* calculate the test statistic */
457 gretl_matrix_multiply(S, iC0, SS);
458 LM = T * (KK1 / 2.0 - gretl_matrix_trace(SS));
459 df = floor(h * K24);
460 tests->val[h-1] = LM;
461 pvals->val[h-1] = chisq_cdf_comp(df, LM);
462 if (!quiet) {
463 pprintf(prn, "lag %*d %9.3f %6d %11.4f\n", hw, h,
464 LM, (int) df, pvals->val[h-1]);
465 }
466 }
467 }
468
469 if (!quiet) {
470 pputc(prn, '\n');
471 }
472
473 bailout:
474
475 gretl_matrix_block_destroy(Bk);
476
477 if (!err) {
478 record_matrix_test_result(tests, pvals);
479 } else {
480 gretl_matrix_free(tests);
481 gretl_matrix_free(pvals);
482 }
483
484 return err;
485 }
486
univariate_arch_test(GRETL_VAR * var,int order,DATASET * dset,gretlopt opt,PRN * prn)487 static int univariate_arch_test (GRETL_VAR *var, int order,
488 DATASET *dset, gretlopt opt,
489 PRN *prn)
490 {
491 gretl_matrix *tests, *pvals;
492 int i, err = 0;
493
494 tests = gretl_column_vector_alloc(var->neqns);
495 pvals = gretl_column_vector_alloc(var->neqns);
496
497 if (tests == NULL || pvals == NULL) {
498 err = E_ALLOC;
499 } else if (!(opt & OPT_I)) {
500 pprintf(prn, "%s %d\n\n", _("Test for ARCH of order"), order);
501 }
502
503 for (i=0; i<var->neqns && !err; i++) {
504 if (!(opt & OPT_I)) {
505 pprintf(prn, "%s %d:\n", _("Equation"), i + 1);
506 }
507 /* add OPT_M for multi-equation output */
508 err = arch_test(var->models[i], order, dset, opt | OPT_M, prn);
509 if (!err) {
510 tests->val[i] = get_last_test_statistic();
511 pvals->val[i] = get_last_pvalue();
512 }
513 }
514
515 if (!err) {
516 record_matrix_test_result(tests, pvals);
517 } else {
518 gretl_matrix_free(tests);
519 gretl_matrix_free(pvals);
520 }
521
522 return err;
523 }
524
gretl_VAR_arch_test(GRETL_VAR * var,int order,DATASET * dset,gretlopt opt,PRN * prn)525 int gretl_VAR_arch_test (GRETL_VAR *var, int order,
526 DATASET *dset, gretlopt opt,
527 PRN *prn)
528 {
529 int h = order;
530 int err;
531
532 if (order == 0) {
533 h = dset->pd;
534 }
535
536 if (opt & OPT_U) {
537 err = univariate_arch_test(var, h, dset, opt, prn);
538 } else {
539 int autoH = (order == 0);
540
541 err = multivariate_arch_test(var, h, autoH, opt, prn);
542 }
543
544 return err;
545 }
546
547 static void
form_C0j(const GRETL_VAR * var,gretl_matrix * C0j,gretl_matrix * et,gretl_matrix * ej,int j)548 form_C0j (const GRETL_VAR *var, gretl_matrix *C0j,
549 gretl_matrix *et, gretl_matrix *ej,
550 int j)
551 {
552 int i, t;
553
554 gretl_matrix_zero(C0j);
555
556 for (t=j; t<var->T; t++) {
557 /* load e_t and e_{t-j} */
558 for (i=0; i<var->neqns; i++) {
559 et->val[i] = gretl_matrix_get(var->E, t, i);
560 ej->val[i] = gretl_matrix_get(var->E, t-j, i);
561 }
562 /* add e_t' * e_{t-j} */
563 gretl_matrix_multiply_mod(et, GRETL_MOD_TRANSPOSE,
564 ej, GRETL_MOD_NONE,
565 C0j, GRETL_MOD_CUMULATE);
566 }
567
568 gretl_matrix_divide_by_scalar(C0j, var->T);
569 }
570
571 /* See S. Johansen, Likelihood-Based Inference in Cointegrated
572 Vector Autoregressive Models, 1995, pp. 21-22
573 */
574
VAR_portmanteau_test(GRETL_VAR * var)575 int VAR_portmanteau_test (GRETL_VAR *var)
576 {
577 gretl_matrix_block *B;
578 gretl_matrix *C00, *C0j;
579 gretl_matrix *et, *ej;
580 gretl_matrix *L, *R, *Tmp;
581 int k, n = var->neqns;
582 double trj, LB = 0.0;
583 int s, j, err = 0;
584
585 var->LB = NADBL;
586 var->LBs = 0;
587
588 /* we'll do this only for unrestricted VARs */
589 if (var->ci == VECM && jrank(var) < var->neqns) {
590 return 0;
591 }
592
593 /* any guidance on the order for this test? */
594 s = var->T / 4;
595 if (s > 48) s = 48;
596
597 k = var->order + (var->ci == VECM);
598 if (s - k <= 0) {
599 /* no degrees of freedom */
600 return 0;
601 }
602
603 B = gretl_matrix_block_new(&C00, n, n,
604 &C0j, n, n,
605 &et, 1, n,
606 &ej, 1, n,
607 &L, n, n,
608 &R, n, n,
609 &Tmp, n, n,
610 NULL);
611
612 if (B == NULL) {
613 return E_ALLOC;
614 }
615
616 form_C0j(var, C00, et, ej, 0);
617 err = gretl_invert_symmetric_matrix(C00);
618
619 for (j=1; j<=s && !err; j++) {
620 form_C0j(var, C0j, et, ej, j);
621 gretl_matrix_multiply(C0j, C00, L);
622 gretl_matrix_multiply_mod(C0j, GRETL_MOD_TRANSPOSE,
623 C00, GRETL_MOD_NONE,
624 R, GRETL_MOD_NONE);
625 gretl_matrix_multiply(L, R, Tmp);
626 trj = gretl_matrix_trace(Tmp);
627 LB += (1.0 / (var->T - j)) * trj;
628 }
629
630 if (!err) {
631 LB *= var->T * (var->T + 2);
632 var->LB = LB;
633 var->LBs = s;
634 }
635
636 gretl_matrix_block_destroy(B);
637
638 return err;
639 }
640
VAR_LR_lag_test(GRETL_VAR * var,const gretl_matrix * E)641 int VAR_LR_lag_test (GRETL_VAR *var, const gretl_matrix *E)
642 {
643 double test_ldet;
644 int err = 0;
645
646 test_ldet = gretl_VAR_ldet(var, E, &err);
647
648 if (!err) {
649 double ll, AIC, BIC, HQC;
650 int T = var->T;
651 int g = var->neqns;
652 int m = var->ncoeff - g;
653 int k = g * m;
654
655 var->LR = T * (test_ldet - var->ldet);
656
657 ll = -(g * T / 2.0) * (LN_2_PI + 1) - (T / 2.0) * test_ldet;
658 AIC = (-2.0 * ll + 2.0 * k) / T;
659 BIC = (-2.0 * ll + k * log(T)) / T;
660 HQC = (-2.0 * ll + 2.0 * k * log(log(T))) / T;
661 var->Ivals[0] = AIC;
662 var->Ivals[1] = BIC;
663 var->Ivals[2] = HQC;
664 }
665
666 return err;
667 }
668
gretl_VAR_print_lagsel(int minlag,gretl_matrix * lltab,gretl_matrix * crittab,int * best_row,PRN * prn)669 static void gretl_VAR_print_lagsel (int minlag,
670 gretl_matrix *lltab,
671 gretl_matrix *crittab,
672 int *best_row,
673 PRN *prn)
674 {
675 int maxlag, nrows = gretl_matrix_rows(crittab);
676 double x;
677 int i, j;
678
679 maxlag = nrows + minlag - 1;
680
681 pprintf(prn, _("VAR system, maximum lag order %d"), maxlag);
682 pputs(prn, "\n\n");
683
684 pputs(prn, _("The asterisks below indicate the best (that is, minimized) values\n"
685 "of the respective information criteria, AIC = Akaike criterion,\n"
686 "BIC = Schwarz Bayesian criterion and HQC = Hannan-Quinn criterion."));
687 pputs(prn, "\n\n");
688
689 pputs(prn, _("lags loglik p(LR) AIC BIC HQC"));
690 pputs(prn, "\n\n");
691
692 for (i=0; i<nrows; i++) {
693 pprintf(prn, "%4d", minlag + i);
694 x = gretl_matrix_get(lltab, i, 0);
695 pprintf(prn, "%14.5f", x);
696 if (i > 0) {
697 x = gretl_matrix_get(lltab, i, 1);
698 pprintf(prn, "%9.5f", x);
699 } else {
700 pputs(prn, " ");
701 }
702 for (j=0; j<N_IVALS; j++) {
703 x = gretl_matrix_get(crittab, i, j);
704 pprintf(prn, "%12.6f", x);
705 if (i == best_row[j]) {
706 pputc(prn, '*');
707 } else {
708 pputc(prn, ' ');
709 }
710 }
711 pputc(prn, '\n');
712 }
713 pputc(prn, '\n');
714 }
715
lagsel_get_min_lag(int p,int * err)716 static int lagsel_get_min_lag (int p, int *err)
717 {
718 int m = get_optval_int(VAR, OPT_M, err);
719
720 if (!*err && (m < 0 || m > p)) {
721 *err = E_INVARG;
722 }
723
724 return m;
725 }
726
727 /* apparatus for selecting the optimal lag length for a VAR */
728
VAR_do_lagsel(GRETL_VAR * var,const DATASET * dset,gretlopt opt,PRN * prn)729 int VAR_do_lagsel (GRETL_VAR *var, const DATASET *dset,
730 gretlopt opt, PRN *prn)
731 {
732 gretl_matrix *crittab = NULL;
733 gretl_matrix *lltab = NULL;
734 gretl_matrix *E = NULL;
735 int p = var->order;
736 int r = p - 1;
737 int T = var->T;
738 int n = var->neqns;
739 /* initialize the "best" at the longest lag */
740 double best[N_IVALS] = { var->AIC, var->BIC, var->HQC };
741 int best_row[N_IVALS] = { r, r, r };
742 double crit[N_IVALS];
743 double LRtest;
744 double ldet = NADBL;
745 int cols0, minlag = 1;
746 int nrows, use_QR = 0;
747 int j, m = 0;
748 int err = 0;
749
750 /* number of cols in X that are not Y lags */
751 cols0 = var->ncoeff - p * n;
752
753 if (opt & OPT_M) {
754 minlag = lagsel_get_min_lag(p, &err);
755 }
756
757 if (p < minlag + 1) {
758 return 0;
759 }
760
761 E = gretl_matrix_alloc(T, n);
762 if (E == NULL) {
763 return E_ALLOC;
764 }
765
766 nrows = p - minlag + 1;
767 crittab = gretl_matrix_alloc(nrows, N_IVALS);
768 lltab = gretl_matrix_alloc(nrows, 2);
769
770 if (crittab == NULL || lltab == NULL) {
771 err = E_ALLOC;
772 goto bailout;
773 }
774
775 if (getenv("VAR_USE_QR") != NULL) {
776 use_QR = 1;
777 }
778
779 for (j=minlag; j<p && !err; j++) {
780 int jxcols = cols0 + j * n;
781
782 if (jxcols == 0) {
783 gretl_matrix_copy_values(E, var->Y);
784 } else {
785 VAR_fill_X(var, j, dset);
786
787 gretl_matrix_reuse(var->X, T, jxcols);
788 gretl_matrix_reuse(var->B, jxcols, n);
789
790 if (use_QR) {
791 err = gretl_matrix_QR_ols(var->Y, var->X, var->B,
792 E, NULL, NULL);
793 } else {
794 err = gretl_matrix_multi_ols(var->Y, var->X, var->B,
795 E, NULL);
796 }
797 }
798
799 if (!err) {
800 ldet = gretl_VAR_ldet(var, E, &err);
801 }
802
803 if (!err) {
804 double ll;
805 int q = var->ncoeff - (n * (p - j));
806 int c, k = n * q;
807
808 ll = -(n * T / 2.0) * (LN_2_PI + 1) - (T / 2.0) * ldet;
809 crit[0] = (-2.0 * ll + 2.0 * k) / T; /* AIC */
810 crit[1] = (-2.0 * ll + k * log(T)) / T; /* BIC */
811 crit[2] = (-2.0 * ll + 2.0 * k * log(log(T))) / T; /* HQC */
812
813 gretl_matrix_set(lltab, m, 0, ll);
814 if (j == minlag) {
815 gretl_matrix_set(lltab, m, 1, 0);
816 } else {
817 LRtest = 2.0 * (ll - gretl_matrix_get(lltab, m-1, 0));
818 gretl_matrix_set(lltab, m, 1, chisq_cdf_comp(n * n, LRtest));
819 }
820
821 for (c=0; c<N_IVALS; c++) {
822 gretl_matrix_set(crittab, m, c, crit[c]);
823 if (crit[c] < best[c]) {
824 best[c] = crit[c];
825 best_row[c] = m;
826 }
827 }
828
829 m++; /* increment table row */
830 }
831 }
832
833 if (!err) {
834 gretl_matrix_set(lltab, m, 0, var->ll);
835 LRtest = 2.0 * (var->ll - gretl_matrix_get(lltab, m - 1, 0));
836 gretl_matrix_set(lltab, m, 1, chisq_cdf_comp(n * n, LRtest));
837 gretl_matrix_set(crittab, m, 0, var->AIC);
838 gretl_matrix_set(crittab, m, 1, var->BIC);
839 gretl_matrix_set(crittab, m, 2, var->HQC);
840 if (!(opt & OPT_S)) {
841 gretl_VAR_print_lagsel(minlag, lltab, crittab, best_row, prn);
842 }
843 record_matrix_test_result(crittab, NULL);
844 crittab = NULL;
845 }
846
847 bailout:
848
849 gretl_matrix_free(crittab);
850 gretl_matrix_free(lltab);
851 gretl_matrix_free(E);
852
853 return err;
854 }
855
VAR_get_hvec(const gretl_matrix * X,const gretl_matrix * XTX,int * err)856 static gretl_matrix *VAR_get_hvec (const gretl_matrix *X,
857 const gretl_matrix *XTX,
858 int *err)
859 {
860 gretl_matrix *hvec;
861 gretl_matrix *Xt;
862 int t, T = X->rows;
863 int k = X->cols;
864 double ht;
865
866 Xt = gretl_matrix_alloc(1, k);
867 hvec = gretl_column_vector_alloc(T);
868
869 if (Xt == NULL || hvec == NULL) {
870 gretl_matrix_free(Xt);
871 gretl_matrix_free(hvec);
872 *err = E_ALLOC;
873 return NULL;
874 }
875
876 for (t=0; t<T && !*err; t++) {
877 gretl_matrix_get_row(X, t, Xt);
878 ht = gretl_scalar_qform(Xt, XTX, err);
879 gretl_vector_set(hvec, t, ht);
880 }
881
882 gretl_matrix_free(Xt);
883
884 return hvec;
885 }
886
var_hac_xox(GRETL_VAR * var,int k,VCVInfo * vi,int * err)887 static gretl_matrix *var_hac_xox (GRETL_VAR *var, int k,
888 VCVInfo *vi, int *err)
889 {
890 gretl_matrix *XOX = NULL;
891 gretl_matrix *uk;
892 int t;
893
894 uk = gretl_column_vector_alloc(var->T);
895
896 if (uk == NULL) {
897 *err = E_ALLOC;
898 } else {
899 for (t=0; t<var->T; t++) {
900 uk->val[t] = gretl_matrix_get(var->E, t, k);
901 }
902 XOX = HAC_XOX(var->X, uk, vi, 0, err);
903 gretl_matrix_free(uk);
904 }
905
906 return XOX;
907 }
908
var_hc_xox(GRETL_VAR * var,int k,int hcv,int * err)909 static gretl_matrix *var_hc_xox (GRETL_VAR *var, int k,
910 int hcv, int *err)
911 {
912 gretl_matrix *XOX;
913 gretl_matrix *hvec = NULL;
914 int T = var->T;
915 int g = var->ncoeff;
916
917 XOX = gretl_matrix_alloc(g, g);
918
919 if (XOX == NULL) {
920 *err = E_ALLOC;
921 } else if (hcv > 1) {
922 hvec = VAR_get_hvec(var->X, var->XTX, err);
923 }
924
925 if (!*err) {
926 /* form X' \Omega X */
927 double xti, xtj, xij;
928 double utk, u2, ht;
929 int i, j, t;
930
931 for (i=0; i<g; i++) {
932 for (j=i; j<g; j++) {
933 xij = 0.0;
934 for (t=0; t<T; t++) {
935 utk = gretl_matrix_get(var->E, t, k);
936 u2 = utk * utk;
937 if (hcv > 1) {
938 ht = gretl_vector_get(hvec, t);
939 u2 /= 1.0 - ht;
940 if (hcv > 2) {
941 u2 /= 1.0 - ht;
942 }
943 }
944 xti = gretl_matrix_get(var->X, t, i);
945 xtj = gretl_matrix_get(var->X, t, j);
946 xij += u2 * xti * xtj;
947 }
948 if (hcv == 1) {
949 xij *= (double) T / (T - g);
950 }
951 gretl_matrix_set(XOX, i, j, xij);
952 if (i != j) {
953 gretl_matrix_set(XOX, j, i, xij);
954 }
955 }
956 }
957 }
958
959 gretl_matrix_free(hvec);
960
961 return XOX;
962 }
963
964 /* (X'X)^{-1} * X'\Omega X * (X'X)^{-1} */
965
VAR_robust_vcv(GRETL_VAR * var,gretl_matrix * V,MODEL * pmod,int hcv,int k)966 static int VAR_robust_vcv (GRETL_VAR *var, gretl_matrix *V,
967 MODEL *pmod, int hcv, int k)
968 {
969 gretl_matrix *XOX = NULL;
970 VCVInfo vi = {0};
971 int err = 0;
972
973 if (var->robust == VAR_HAC) {
974 XOX = var_hac_xox(var, k, &vi, &err);
975 } else {
976 XOX = var_hc_xox(var, k, hcv, &err);
977 }
978
979 if (!err) {
980 gretl_matrix_qform(var->XTX, GRETL_MOD_TRANSPOSE, XOX,
981 V, GRETL_MOD_NONE);
982
983 if (var->robust == VAR_HAC) {
984 gretl_model_set_full_vcv_info(pmod, VCV_HAC, vi.vmin,
985 vi.order, vi.flags,
986 vi.bw);
987 } else {
988 gretl_model_set_vcv_info(pmod, VCV_HC, hcv);
989 }
990 }
991
992 gretl_matrix_free(XOX);
993
994 return err;
995 }
996
997 /* Run the various per-equation omit tests (all lags of each var in
998 turn, last lag of all vars) using the Wald method. We also
999 add the standard errors to the models here, since we have the
1000 per-equation covariance matrices to hand.
1001 */
1002
VAR_wald_omit_tests(GRETL_VAR * var)1003 int VAR_wald_omit_tests (GRETL_VAR *var)
1004 {
1005 gretl_matrix *V = NULL;
1006 gretl_matrix *C = NULL;
1007 gretl_vector *b = NULL;
1008 int hcv = libset_get_int(HC_VERSION);
1009 int p = (var->lags != NULL)? var->lags[0] : var->order;
1010 int n = var->neqns;
1011 int g = var->ncoeff;
1012 int dim = (p > n)? p : n;
1013 int i, j, k, m = 0;
1014 int any_F_err = 0;
1015 int err = 0;
1016
1017 if (var->ifc && var->robust && g - 1 > dim) {
1018 /* need bigger arrays for robust overall F-test */
1019 dim = g - 1;
1020 }
1021
1022 V = gretl_matrix_alloc(g, g);
1023 C = gretl_matrix_alloc(dim, dim);
1024 b = gretl_column_vector_alloc(dim);
1025
1026 if (V == NULL || C == NULL || b == NULL) {
1027 return E_ALLOC;
1028 }
1029
1030 for (i=0; i<var->neqns && !err; i++) {
1031 MODEL *pmod = var->models[i];
1032 int ii, jj, jpos, ipos = var->ifc;
1033 double w, vij;
1034 int F_err = 0;
1035
1036 if (var->robust) {
1037 err = VAR_robust_vcv(var, V, pmod, hcv, i);
1038 } else {
1039 gretl_matrix_copy_values(V, var->XTX);
1040 gretl_matrix_multiply_by_scalar(V, pmod->sigma * pmod->sigma);
1041 }
1042
1043 if (err) {
1044 break;
1045 }
1046
1047 /* set the per-model VCV */
1048 gretl_model_write_vcv(pmod, V);
1049
1050 gretl_matrix_reuse(C, p, p);
1051 gretl_matrix_reuse(b, p, 1);
1052
1053 /* exclusion of each var, all lags */
1054 if (p > 0) {
1055 for (j=0; j<var->neqns; j++) {
1056 double w = NADBL;
1057
1058 gretl_matrix_extract_matrix(C, V, ipos, ipos, GRETL_MOD_NONE);
1059 for (k=0; k<p; k++) {
1060 b->val[k] = pmod->coeff[k + ipos];
1061 }
1062 F_err = gretl_invert_symmetric_matrix(C);
1063 if (!F_err) {
1064 w = gretl_scalar_qform(b, C, &F_err);
1065 }
1066 if (F_err) {
1067 any_F_err = 1;
1068 var->Fvals[m++] = NADBL;
1069 } else {
1070 var->Fvals[m++] = w / p;
1071 }
1072 ipos += p;
1073 }
1074 }
1075
1076 /* exclusion of last lag, all vars? */
1077 if (p > 0) {
1078 gretl_matrix_reuse(C, n, n);
1079 gretl_matrix_reuse(b, n, 1);
1080
1081 ipos = var->ifc + p - 1;
1082 for (ii=0; ii<n; ii++) {
1083 jpos = var->ifc + p - 1;
1084 for (jj=0; jj<n; jj++) {
1085 vij = gretl_matrix_get(V, ipos, jpos);
1086 gretl_matrix_set(C, ii, jj, vij);
1087 jpos += p;
1088 }
1089 b->val[ii] = pmod->coeff[ipos];
1090 ipos += p;
1091 }
1092
1093 F_err = gretl_invert_symmetric_matrix(C);
1094 if (!F_err) {
1095 w = gretl_scalar_qform(b, C, &F_err);
1096 }
1097 if (F_err) {
1098 any_F_err = 1;
1099 var->Fvals[m++] = NADBL;
1100 } else {
1101 var->Fvals[m++] = w / n;
1102 }
1103 }
1104
1105 /* exclusion of all but const? */
1106 if (var->ifc && var->robust) {
1107 gretl_matrix_reuse(C, g-1, g-1);
1108 gretl_matrix_reuse(b, g-1, 1);
1109
1110 gretl_matrix_extract_matrix(C, V, 1, 1, GRETL_MOD_NONE);
1111 for (k=0; k<g-1; k++) {
1112 b->val[k] = pmod->coeff[k+1];
1113 }
1114 F_err = gretl_invert_symmetric_matrix(C);
1115 if (!F_err) {
1116 w = gretl_scalar_qform(b, C, &F_err);
1117 }
1118 if (F_err) {
1119 any_F_err = 1;
1120 pmod->fstt = NADBL;
1121 } else {
1122 pmod->fstt = w / (g-1);
1123 }
1124 }
1125 }
1126
1127 gretl_matrix_free(V);
1128 gretl_matrix_free(C);
1129 gretl_matrix_free(b);
1130
1131 if (!err && any_F_err) {
1132 fprintf(stderr, "*** Warning: some F-tests could not be computed\n");
1133 }
1134
1135 return err;
1136 }
1137
1138 #define VO_DEBUG 0
1139
gretl_VAR_get_exo_list(const GRETL_VAR * var)1140 const int *gretl_VAR_get_exo_list (const GRETL_VAR *var)
1141 {
1142 return (var == NULL)? NULL : var->xlist;
1143 }
1144
gretl_VAR_get_endo_list(const GRETL_VAR * var)1145 const int *gretl_VAR_get_endo_list (const GRETL_VAR *var)
1146 {
1147 return (var == NULL)? NULL : var->ylist;
1148 }
1149
1150 /* Based on the specification stored in the original VAR struct,
1151 plus a new exogenous list, constitute a full VAR list.
1152 */
1153
build_VAR_list(const GRETL_VAR * var,int * exolist,int * err)1154 static int *build_VAR_list (const GRETL_VAR *var, int *exolist, int *err)
1155 {
1156 return VAR_list_composite(var->ylist, exolist, var->rlist);
1157 }
1158
gretl_VAR_real_omit_test(const GRETL_VAR * orig,const GRETL_VAR * new,const DATASET * dset,gretlopt opt,PRN * prn)1159 static int gretl_VAR_real_omit_test (const GRETL_VAR *orig,
1160 const GRETL_VAR *new,
1161 const DATASET *dset,
1162 gretlopt opt,
1163 PRN *prn)
1164 {
1165 int *omitlist = NULL;
1166 double LR, pval;
1167 int df, nr = 0;
1168 int err = 0;
1169
1170 #if VO_DEBUG
1171 fprintf(stderr, "gretl_VAR_real_omit_test: about to diff lists\n");
1172 printlist(orig->xlist, "orig xlist");
1173 printlist(new->xlist, "new xlist");
1174 #endif
1175
1176 if (orig->xlist != NULL) {
1177 if (new->xlist == NULL) {
1178 omitlist = gretl_list_copy(orig->xlist);
1179 } else {
1180 omitlist = gretl_list_diff_new(orig->xlist, new->xlist, 1);
1181 }
1182 if (omitlist == NULL) {
1183 return E_ALLOC;
1184 }
1185 nr = omitlist[0];
1186 }
1187 if (opt & OPT_E) {
1188 /* omitting seasonals */
1189 nr += dset->pd + 1;
1190 }
1191 if (opt & OPT_T) {
1192 /* omitting trend */
1193 nr++;
1194 }
1195
1196 LR = orig->T * (new->ldet - orig->ldet);
1197 df = orig->neqns * nr;
1198 pval = chisq_cdf_comp(df, LR);
1199
1200 record_test_result(LR, pval);
1201
1202 pprintf(prn, "%s:\n", _("Test on the original VAR"));
1203
1204 print_add_omit_null(omitlist, dset, opt | OPT_S, prn);
1205
1206 pprintf(prn, " %s: %s(%d) = %g, ", _("LR test"),
1207 _("Chi-square"), df, LR);
1208 pprintf(prn, _("with p-value = %g\n"), pval);
1209
1210 free(omitlist);
1211
1212 return err;
1213 }
1214
VAR_omit_check(GRETL_VAR * var,const int * omitlist,gretlopt opt)1215 static int VAR_omit_check (GRETL_VAR *var, const int *omitlist,
1216 gretlopt opt)
1217 {
1218 int nx = omitlist == NULL ? 0 : omitlist[0];
1219 int err = 0;
1220
1221 if (nx == 0 && !(opt & (OPT_T | OPT_E))) {
1222 /* nothing to be omitted */
1223 err = E_NOOMIT;
1224 } else if ((opt & OPT_T) && !(var->detflags & DET_TREND)) {
1225 /* can't have the --trend option with no auto-trend */
1226 err = E_BADOPT;
1227 } else if ((opt & OPT_E) && !(var->detflags & DET_SEAS)) {
1228 /* can't have the --seasonals option with no auto-seasonals */
1229 err = E_BADOPT;
1230 }
1231
1232 return err;
1233 }
1234
1235 /**
1236 * gretl_VAR_omit_test:
1237 * @var: pointer to original VAR.
1238 * @omitlist: list of variables to omit from original model.
1239 * @dset: dataset struct.
1240 * @opt: option flags.
1241 * @prn: gretl printing struct.
1242 * @err: location to receive error code.
1243 *
1244 * Re-estimates a given VAR after removing the variables
1245 * specified in @omitlist, and reports per-equation F-tests
1246 * and system-wide LR tests for the null hypothesis that
1247 * the omitted variables have zero parameters.
1248 *
1249 * Returns: restricted VAR on success, %NULL on error.
1250 */
1251
gretl_VAR_omit_test(GRETL_VAR * var,const int * omitlist,DATASET * dset,gretlopt opt,PRN * prn,int * err)1252 GRETL_VAR *gretl_VAR_omit_test (GRETL_VAR *var, const int *omitlist,
1253 DATASET *dset, gretlopt opt,
1254 PRN *prn, int *err)
1255 {
1256 GRETL_VAR *vnew = NULL;
1257 gretlopt varopt = OPT_NONE;
1258 int smpl_t1 = dset->t1;
1259 int smpl_t2 = dset->t2;
1260 int *tmplist = NULL;
1261 int *varlist = NULL;
1262 int c1 = 0;
1263
1264 if (var == NULL) {
1265 *err = E_DATA;
1266 return NULL;
1267 }
1268
1269 *err = VAR_omit_check(var, omitlist, opt);
1270 if (*err) {
1271 return NULL;
1272 }
1273
1274 if (var->ifc) {
1275 c1 = !gretl_list_const_pos(omitlist, 1, dset);
1276 }
1277
1278 if (omitlist != NULL && omitlist[0] > 0) {
1279 /* create reduced exogenous vars list for test VAR */
1280 tmplist = gretl_list_omit(var->xlist, omitlist, 1, err);
1281 if (tmplist == NULL) {
1282 goto bailout;
1283 }
1284 } else if (var->xlist != NULL) {
1285 tmplist = gretl_list_copy(var->xlist);
1286 if (tmplist == NULL) {
1287 goto bailout;
1288 }
1289 }
1290
1291 /* create full input VAR list for test VAR */
1292 varlist = build_VAR_list(var, tmplist, err);
1293 if (varlist == NULL) {
1294 goto bailout;
1295 }
1296
1297 if ((var->detflags & DET_SEAS) && !(opt & OPT_E)) {
1298 varopt |= OPT_D;
1299 }
1300
1301 if ((var->detflags & DET_TREND) && !(opt & OPT_T)) {
1302 varopt |= OPT_T;
1303 }
1304
1305 /* If the original VAR did not include a constant, we need to
1306 pass OPT_N to the test VAR to suppress the constant.
1307 We also need to pass OPT_N in case the constant was
1308 present originally but is now to be omitted.
1309 */
1310 if (var->ifc == 0 || c1 == 0) {
1311 varopt |= OPT_N;
1312 }
1313
1314 /* impose as sample range the estimation range of the
1315 original VAR */
1316 dset->t1 = var->t1;
1317 dset->t2 = var->t2;
1318
1319 vnew = gretl_VAR(var->order, var->lags, varlist, dset,
1320 varopt, NULL, err);
1321
1322 if (vnew != NULL) {
1323 /* do the actual test(s) */
1324 *err = gretl_VAR_real_omit_test(var, vnew, dset, opt, prn);
1325 if (!*err && prn != NULL) {
1326 gretl_VAR_print(vnew, dset, OPT_NONE, prn);
1327 }
1328 }
1329
1330 /* put back into dset what was there on input */
1331 dset->t1 = smpl_t1;
1332 dset->t2 = smpl_t2;
1333
1334 bailout:
1335
1336 free(tmplist);
1337 free(varlist);
1338
1339 return vnew;
1340 }
1341
set_wald_R_ones(gretl_matrix * R,int k,int r0,int neq,int stride)1342 static int set_wald_R_ones (gretl_matrix *R, int k, int r0, int neq,
1343 int stride)
1344 {
1345 if (k < 0) {
1346 return E_DATA;
1347 } else {
1348 int i, rmax = r0 + neq;
1349
1350 for (i=r0; i<rmax; i++) {
1351 gretl_matrix_set(R, i, k, 1.0);
1352 k += stride;
1353 }
1354 return 0;
1355 }
1356 }
1357
1358 /**
1359 * gretl_VAR_wald_omit_test:
1360 * @var: pointer to original.
1361 * @omitlist: list of variables to omit from original model.
1362 * @dset: dataset struct.
1363 * @opt: option flags.
1364 * @prn: gretl printing struct.
1365 *
1366 * Runs a Wald test for the joint omission of the exogenous
1367 * variables specified in @omitlist.
1368 *
1369 * Returns: 0 on successful completion, non-zero error code
1370 * otherwise.
1371 */
1372
gretl_VAR_wald_omit_test(GRETL_VAR * var,const int * omitlist,DATASET * dset,gretlopt opt,PRN * prn)1373 int gretl_VAR_wald_omit_test (GRETL_VAR *var, const int *omitlist,
1374 DATASET *dset, gretlopt opt,
1375 PRN *prn)
1376 {
1377 gretl_matrix *B = NULL;
1378 gretl_matrix *S = NULL;
1379 gretl_matrix *V = NULL;
1380 gretl_matrix *R = NULL;
1381 gretl_matrix *RB = NULL;
1382 gretl_matrix *RVR = NULL;
1383 double x, test = NADBL;
1384 int i, j, neq, eqk, row0;
1385 int pos, nr = 0, nx_omit = 0;
1386 int err = 0;
1387
1388 if (var == NULL || var->B == NULL ||
1389 var->S == NULL || var->XTX == NULL) {
1390 return E_DATA;
1391 }
1392
1393 err = VAR_omit_check(var, omitlist, opt);
1394 if (err) {
1395 return err;
1396 }
1397
1398 if (omitlist != NULL) {
1399 nx_omit = omitlist[0];
1400 }
1401
1402 B = gretl_matrix_vectorize_new(var->B);
1403 if (B == NULL) {
1404 return E_ALLOC;
1405 }
1406
1407 neq = var->neqns;
1408
1409 if (nx_omit > 0) {
1410 /* omitting user-specified exogenous vars */
1411 nr += nx_omit * neq;
1412 }
1413 if (opt & OPT_T) {
1414 /* omit trend */
1415 nr += neq;
1416 }
1417 if (opt & OPT_E) {
1418 /* omit seasonals */
1419 nr += (dset->pd - 1) * neq;
1420 }
1421
1422 S = gretl_zero_matrix_new(neq, neq);
1423 R = gretl_zero_matrix_new(nr, B->rows);
1424 RB = gretl_matrix_alloc(nr, 1);
1425 RVR = gretl_matrix_alloc(nr, nr);
1426
1427 if (S == NULL || R == NULL || RB == NULL || RVR == NULL) {
1428 err = E_ALLOC;
1429 goto bailout;
1430 }
1431
1432 for (i=0; i<neq; i++) {
1433 x = gretl_matrix_get(var->S, i, i);
1434 gretl_matrix_set(S, i, i, x);
1435 }
1436
1437 V = gretl_matrix_kronecker_product_new(S, var->XTX, &err);
1438 if (err) {
1439 goto bailout;
1440 }
1441
1442 eqk = var->B->rows;
1443 pos = var->ifc + neq * var->order;
1444 row0 = 0;
1445
1446 if (nx_omit > 0) {
1447 /* @omitlist holds ID numbers of exog variables */
1448 int vi;
1449
1450 for (i=1; i<=omitlist[0]; i++) {
1451 vi = omitlist[i];
1452 if (vi == 0) {
1453 j = 0;
1454 } else {
1455 j = in_gretl_list(var->xlist, vi);
1456 j += pos - 1;
1457 }
1458 set_wald_R_ones(R, j, row0, neq, eqk);
1459 row0 += neq;
1460 }
1461 }
1462 if (opt & OPT_E) {
1463 /* omit auto-seasonals */
1464 for (i=1; i<dset->pd; i++) {
1465 j = pos + i - 1;
1466 set_wald_R_ones(R, j, row0, neq, eqk);
1467 row0 += neq;
1468 }
1469 }
1470 if (opt & OPT_T) {
1471 /* omit auto-trend (comes as last param) */
1472 j = eqk - 1;
1473 set_wald_R_ones(R, j, row0, neq, eqk);
1474 row0 += neq;
1475 }
1476
1477 gretl_matrix_multiply(R, B, RB);
1478 gretl_matrix_qform(R, GRETL_MOD_NONE, V,
1479 RVR, GRETL_MOD_NONE);
1480 err = gretl_invert_symmetric_matrix(RVR);
1481
1482 if (!err) {
1483 test = gretl_scalar_qform(RB, RVR, &err);
1484 if (!err) {
1485 double pval = chisq_cdf_comp(nr, test);
1486
1487 record_test_result(test, pval);
1488
1489 if (!(opt & OPT_I)) {
1490 /* not silent */
1491 pprintf(prn, "%s:\n", _("Test on VAR"));
1492 print_add_omit_null(omitlist, dset, opt | OPT_S, prn);
1493 pprintf(prn, " %s: %s(%d) = %g, %s %g\n\n", _("Wald test"),
1494 _("Chi-square"), nr, test,
1495 _("p-value"), pval);
1496 }
1497 }
1498 }
1499
1500 bailout:
1501
1502 gretl_matrix_free(B);
1503 gretl_matrix_free(S);
1504 gretl_matrix_free(V);
1505 gretl_matrix_free(R);
1506 gretl_matrix_free(RB);
1507 gretl_matrix_free(RVR);
1508
1509 return err;
1510 }
1511