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 #define FULL_XML_HEADERS
21
22 #include "libgretl.h"
23 #include "var.h"
24 #include "johansen.h"
25 #include "varprint.h"
26 #include "vartest.h"
27 #include "libset.h"
28 #include "transforms.h"
29 #include "gretl_xml.h"
30 #include "matrix_extra.h"
31 #include "system.h"
32
33 #define VDEBUG 0
34
35 #define VAR_SE_DFCORR 1
36
37 enum {
38 VAR_ESTIMATE = 1,
39 VAR_LAGSEL,
40 VECM_ESTIMATE,
41 VECM_CTEST
42 };
43
44 static JohansenInfo *
45 johansen_info_new (GRETL_VAR *var, int rank, gretlopt opt);
46
47 static gretl_matrix *
48 gretl_VAR_get_fcast_se (GRETL_VAR *var, int periods);
49
VAR_add_models(GRETL_VAR * var,const DATASET * dset)50 static int VAR_add_models (GRETL_VAR *var, const DATASET *dset)
51 {
52 int n = var->neqns;
53 int i, err = 0;
54
55 if (var->models != NULL) {
56 for (i=0; i<n; i++) {
57 clear_model(var->models[i]);
58 }
59 } else {
60 var->models = gretl_model_array_new(n);
61 if (var->models == NULL) {
62 err = E_ALLOC;
63 }
64 }
65
66 if (var->models != NULL) {
67 for (i=0; i<n; i++) {
68 gretl_model_smpl_init(var->models[i], dset);
69 }
70 }
71
72 return err;
73 }
74
75 /* Here we allocate storage for a full-size companion
76 matrix and fill in the "boilerplate" rows at the
77 bottom. The initial (VAR coefficient) rows are left
78 as zeros.
79 */
80
VAR_allocate_companion_matrix(GRETL_VAR * var)81 static int VAR_allocate_companion_matrix (GRETL_VAR *var)
82 {
83 int n, dim, err = 0;
84
85 if (var->A != NULL) {
86 /* the job is already done */
87 return 0;
88 }
89
90 n = var->neqns;
91 dim = n * effective_order(var);
92 var->A = gretl_zero_matrix_new(dim, dim);
93
94 if (var->A == NULL) {
95 err = E_ALLOC;
96 } else {
97 int j, jmax = dim - n;
98
99 for (j=0; j<jmax; j++) {
100 gretl_matrix_set(var->A, j+n, j, 1);
101 }
102 }
103
104 return err;
105 }
106
107 /* Given an m x n matrix @A, with m < n, construct the full
108 n x n companion matrix.
109 */
110
companionize(const gretl_matrix * A,int * err)111 static gretl_matrix *companionize (const gretl_matrix *A,
112 int *err)
113 {
114 gretl_matrix *ret = NULL;
115 int m = A->rows, n = A->cols;
116 int d = n - m;
117
118 if (d < 0) {
119 *err = E_INVARG;
120 return NULL;
121 }
122
123 if (d == 0) {
124 /* shouldn't get here, but just in case */
125 ret = gretl_matrix_copy(A);
126 } else {
127 double aij;
128 int i, j;
129
130 ret = gretl_zero_matrix_new(n, n);
131 if (ret != NULL) {
132 for (j=0; j<n; j++) {
133 for (i=0; i<m; i++) {
134 aij = gretl_matrix_get(A, i, j);
135 gretl_matrix_set(ret, i, j, aij);
136 }
137 if (j < d) {
138 gretl_matrix_set(ret, j+m, j, 1);
139 }
140 }
141 }
142 }
143
144 if (ret == NULL) {
145 *err = E_ALLOC;
146 }
147
148 return ret;
149 }
150
151 /* Given the full (square) companion matrix, return a version
152 without the trailing rows holding (I ~ 0). The inverse
153 operation to companionize() above.
154 */
155
decompanionize(const gretl_matrix * A,int neqns,GretlMatrixMod mod)156 gretl_matrix *decompanionize (const gretl_matrix *A, int neqns,
157 GretlMatrixMod mod)
158 {
159 gretl_matrix *ret;
160 int np = A->cols;
161
162 if (mod == GRETL_MOD_TRANSPOSE) {
163 ret = gretl_matrix_alloc(np, neqns);
164 } else {
165 ret = gretl_matrix_alloc(neqns, np);
166 }
167
168 if (ret != NULL) {
169 double aij;
170 int i, j;
171
172 for (j=0; j<np; j++) {
173 for (i=0; i<neqns; i++) {
174 aij = gretl_matrix_get(A, i, j);
175 if (mod == GRETL_MOD_TRANSPOSE) {
176 gretl_matrix_set(ret, j, i, aij);
177 } else {
178 gretl_matrix_set(ret, i, j, aij);
179 }
180 }
181 }
182 }
183
184 return ret;
185 }
186
VAR_allocate_cholesky_matrix(GRETL_VAR * var)187 static int VAR_allocate_cholesky_matrix (GRETL_VAR *var)
188 {
189 int n = var->neqns;
190 int err = 0;
191
192 if (var->C != NULL) {
193 /* the job is already done */
194 return 0;
195 }
196
197 var->C = gretl_zero_matrix_new(n, n);
198 if (var->C == NULL) {
199 err = E_ALLOC;
200 }
201
202 return err;
203 }
204
VAR_allocate_residuals_matrix(GRETL_VAR * var)205 static int VAR_allocate_residuals_matrix (GRETL_VAR *var)
206 {
207 int err = 0;
208
209 if (var->E != NULL) {
210 /* the job is already done */
211 return 0;
212 }
213
214 var->E = gretl_zero_matrix_new(var->T, var->neqns);
215 if (var->E == NULL) {
216 err = E_ALLOC;
217 }
218
219 return err;
220 }
221
VAR_add_XTX_matrix(GRETL_VAR * var)222 static int VAR_add_XTX_matrix (GRETL_VAR *var)
223 {
224 int k = var->X->cols;
225 int err = 0;
226
227 var->XTX = gretl_zero_matrix_new(k, k);
228
229 if (var->XTX == NULL) {
230 err = E_ALLOC;
231 } else {
232 gretl_matrix_multiply_mod(var->X, GRETL_MOD_TRANSPOSE,
233 var->X, GRETL_MOD_NONE,
234 var->XTX, GRETL_MOD_NONE);
235 err = gretl_invert_matrix(var->XTX);
236 }
237
238 return err;
239 }
240
n_restricted_terms(const GRETL_VAR * v)241 int n_restricted_terms (const GRETL_VAR *v)
242 {
243 int n = 0;
244
245 if (v->jinfo != NULL &&
246 (v->jinfo->code == J_REST_CONST ||
247 v->jinfo->code == J_REST_TREND)) {
248 n++;
249 }
250
251 if (v->rlist != NULL) {
252 n += v->rlist[0];
253 }
254
255 return n;
256 }
257
gretl_VAR_clear(GRETL_VAR * var)258 void gretl_VAR_clear (GRETL_VAR *var)
259 {
260 var->ci = 0;
261 var->refcount = 0;
262 var->err = 0;
263 var->neqns = var->order = 0;
264 var->t1 = var->t2 = var->T = var->df = 0;
265 var->ifc = var->ncoeff = 0;
266 var->detflags = 0;
267 var->robust = 0;
268 var->LBs = 0;
269 var->xcols = 0;
270
271 var->lags = NULL;
272 var->ylist = NULL;
273 var->xlist = NULL;
274 var->rlist = NULL;
275
276 var->Y = NULL;
277 var->X = NULL;
278 var->B = NULL;
279 var->XTX = NULL;
280 var->A = NULL;
281 var->L = NULL;
282 var->E = NULL;
283 var->C = NULL;
284 var->S = NULL;
285 var->F = NULL;
286 var->V = NULL;
287 var->ord = NULL;
288
289 var->models = NULL;
290 var->Fvals = NULL;
291 var->Ivals = NULL;
292
293 var->ll = var->ldet = var->LR = NADBL;
294 var->AIC = var->BIC = var->HQC = NADBL;
295 var->LB = NADBL;
296
297 var->jinfo = NULL;
298 var->name = NULL;
299 }
300
301 #define lag_wanted(v, i) (v->lags == NULL || in_gretl_list(v->lags, i))
302
303 /* Construct the common X matrix (composed of lags of the core
304 variables plus other terms if applicable). This is in
305 common between VARs and VECMs
306 */
307
VAR_fill_X(GRETL_VAR * v,int p,const DATASET * dset)308 void VAR_fill_X (GRETL_VAR *v, int p, const DATASET *dset)
309 {
310 const double *x;
311 int diff = (v->ci == VECM);
312 int i, j, s, t, vi;
313 int k = 0; /* X column index */
314
315 /* const first */
316 if (v->detflags & DET_CONST) {
317 s = 0;
318 for (t=v->t1; t<=v->t2; t++) {
319 gretl_matrix_set(v->X, s++, k, 1.0);
320 }
321 k++;
322 }
323
324 /* add lagged Ys */
325 for (i=0; i<v->neqns; i++) {
326 vi = v->ylist[i+1];
327 for (j=1; j<=p; j++) {
328 if (!lag_wanted(v, j)) {
329 continue;
330 }
331 x = dset->Z[vi];
332 s = 0;
333 for (t=v->t1; t<=v->t2; t++) {
334 if (diff) {
335 gretl_matrix_set(v->X, s++, k, x[t-j] - x[t-j-1]);
336 } else {
337 gretl_matrix_set(v->X, s++, k, x[t-j]);
338 }
339 }
340 k++;
341 }
342 }
343
344 /* add any exogenous vars */
345 if (v->xlist != NULL) {
346 for (i=1; i<=v->xlist[0]; i++) {
347 vi = v->xlist[i];
348 s = 0;
349 for (t=v->t1; t<=v->t2; t++) {
350 gretl_matrix_set(v->X, s++, k, dset->Z[vi][t]);
351 }
352 k++;
353 }
354 }
355
356 /* add other deterministic terms */
357 if (v->detflags & DET_SEAS) {
358 int per = get_subperiod(v->t1, dset, NULL);
359 int pd1 = dset->pd - 1;
360 double s0, s1;
361
362 if (v->ci == VECM) {
363 s1 = 1 - 1.0 / dset->pd;
364 s0 = s1 - 1;
365 } else {
366 s1 = 1;
367 s0 = 0;
368 }
369
370 for (t=0; t<v->T; t++) {
371 for (i=0; i<pd1; i++) {
372 gretl_matrix_set(v->X, t, k+i, (per == i)? s1 : s0);
373 }
374 per = (per < pd1)? per + 1 : 0;
375 }
376 k += pd1;
377 }
378
379 if (v->detflags & DET_TREND) {
380 s = 0;
381 for (t=v->t1; t<=v->t2; t++) {
382 gretl_matrix_set(v->X, s++, k, (double) (t + 1));
383 }
384 k++;
385 }
386
387 if (v->X != NULL) {
388 gretl_matrix_set_t1(v->X, v->t1);
389 gretl_matrix_set_t2(v->X, v->t2);
390 }
391
392 #if VDEBUG
393 gretl_matrix_print(v->X, "X");
394 #endif
395 }
396
397 /* construct the matrix of dependent variables for a plain VAR */
398
VAR_fill_Y(GRETL_VAR * v,const DATASET * dset)399 static void VAR_fill_Y (GRETL_VAR *v, const DATASET *dset)
400 {
401 int i, vi, s, t;
402
403 for (i=0; i<v->neqns; i++) {
404 vi = v->ylist[i+1];
405 s = 0;
406 for (t=v->t1; t<=v->t2; t++) {
407 gretl_matrix_set(v->Y, s++, i, dset->Z[vi][t]);
408 }
409 }
410
411 if (v->Y != NULL) {
412 gretl_matrix_set_t1(v->Y, v->t1);
413 gretl_matrix_set_t2(v->Y, v->t2);
414 }
415
416 #if VDEBUG
417 gretl_matrix_print(v->Y, "Y");
418 #endif
419 }
420
421 /* construct the combined matrix of dependent variables for a VECM */
422
VECM_fill_Y(GRETL_VAR * v,const DATASET * dset,gretl_matrix * Y)423 static void VECM_fill_Y (GRETL_VAR *v, const DATASET *dset,
424 gretl_matrix *Y)
425 {
426 const double *yi;
427 int i, vi, s, t;
428 int k = 0;
429
430 for (i=0; i<v->neqns; i++) {
431 vi = v->ylist[i+1];
432 yi = dset->Z[vi];
433 k = i + v->neqns;
434 s = 0;
435 for (t=v->t1; t<=v->t2; t++) {
436 gretl_matrix_set(Y, s, i, yi[t] - yi[t-1]);
437 gretl_matrix_set(Y, s, k, yi[t-1]);
438 s++;
439 }
440 }
441
442 if (auto_restr(v)) {
443 int trend = (v->jinfo->code == J_REST_TREND);
444
445 k++;
446 for (t=0; t<v->T; t++) {
447 gretl_matrix_set(Y, t, k, trend ? (v->t1 + t) : 1);
448 }
449 }
450
451 if (v->rlist != NULL) {
452 /* There's room for debate here (?) over whether we should
453 use the current value or the first lag of restricted
454 exogenous variables. But using the current value agrees
455 with Ox. See also VAR_set_sample().
456 */
457 for (i=0; i<v->rlist[0]; i++) {
458 vi = v->rlist[i+1];
459 k++;
460 s = 0;
461 for (t=v->t1; t<=v->t2; t++) {
462 gretl_matrix_set(Y, s++, k, dset->Z[vi][t]); /* was t-1 */
463 }
464 }
465 }
466
467 #if VDEBUG
468 gretl_matrix_print(Y, "VECM Y");
469 #endif
470 }
471
472 /* Split the user-supplied list, if need be, and construct the lists
473 of endogenous and (possibly) exogenous vars. Note that
474 deterministic terms are handled separately, via option flags.
475 */
476
VAR_make_lists(GRETL_VAR * v,const int * list,const DATASET * dset)477 static int VAR_make_lists (GRETL_VAR *v, const int *list,
478 const DATASET *dset)
479 {
480 int err = 0;
481
482 #if VDEBUG
483 printlist(list, "VAR_make_lists: incoming list");
484 #endif
485
486 if (gretl_list_has_separator(list)) {
487 err = gretl_list_split_on_separator(list, &v->ylist, &v->xlist);
488 } else {
489 v->ylist = gretl_list_copy(list);
490 if (v->ylist == NULL) {
491 err = E_ALLOC;
492 }
493 }
494
495 if (!err && (v->ylist == NULL || v->ylist[0] < 1)) {
496 /* first test for at least 1 endog var */
497 err = E_ARGS;
498 }
499
500 if (!err && v->ci == VECM) {
501 if (v->xlist != NULL && gretl_list_has_separator(v->xlist)) {
502 int *tmp = NULL;
503
504 err = gretl_list_split_on_separator(v->xlist, &tmp, &v->rlist);
505 if (!err) {
506 free(v->xlist);
507 v->xlist = tmp;
508 }
509 }
510 }
511
512 if (!err) {
513 gretl_list_purge_const(v->ylist, dset);
514 if (v->xlist != NULL) {
515 gretl_list_purge_const(v->xlist, dset);
516 }
517 if (v->rlist != NULL) {
518 gretl_list_purge_const(v->rlist, dset);
519 }
520 }
521
522 if (!err && (v->ylist == NULL || v->ylist[0] < 1)) {
523 /* re-test after (possibly) losing const */
524 err = E_ARGS;
525 }
526
527 #if VDEBUG
528 printlist(v->ylist, "v->ylist");
529 printlist(v->xlist, "v->xlist");
530 printlist(v->rlist, "v->rlist");
531 #endif
532
533 return err;
534 }
535
536 /* Starting from the given sample range, construct the feasible
537 estimation range for a VAR or VECM. Flag an error if there
538 are missing values within the sample period.
539 Also see system_adjust_t1t2 in lib/src/system.c
540 (perhaps the two could be consolidated at some point)
541 */
542
VAR_set_sample(GRETL_VAR * v,const DATASET * dset)543 static int VAR_set_sample (GRETL_VAR *v, const DATASET *dset)
544 {
545 int diff = (v->ci == VECM)? 1 : 0;
546 int i, vi, t, err = 0;
547
548 /* advance t1 if needed */
549
550 for (t=dset->t1; t<=dset->t2; t++) {
551 int miss = 0, p, s = t - (v->order + diff);
552
553 for (i=1; i<=v->ylist[0] && !miss; i++) {
554 vi = v->ylist[i];
555 if (na(dset->Z[vi][t]) || s < 0) {
556 v->t1 += 1;
557 miss = 1;
558 }
559 for (p=s; p<t && !miss; p++) {
560 if (na(dset->Z[vi][p])) {
561 v->t1 += 1;
562 miss = 1;
563 }
564 }
565 }
566
567 if (v->xlist != NULL && !miss) {
568 for (i=1; i<=v->xlist[0] && !miss; i++) {
569 vi = v->xlist[i];
570 if (na(dset->Z[vi][t])) {
571 v->t1 += 1;
572 miss = 1;
573 }
574 }
575 }
576
577 /* In estimating a VECM we need the level (or first lag?)
578 of each restricted exogenous var to serve as dependent
579 variable in one of the initial OLS equations. See also
580 VECM_fill_Y().
581 */
582
583 if (v->rlist != NULL && !miss) {
584 for (i=1; i<=v->rlist[0] && !miss; i++) {
585 vi = v->rlist[i];
586 if (na(dset->Z[vi][t])) {
587 v->t1 += 1;
588 miss = 1;
589 }
590 }
591 }
592
593 if (!miss) {
594 break;
595 }
596 }
597
598 /* retard t2 if needed */
599
600 for (t=dset->t2; t>=v->t1; t--) {
601 int miss = 0;
602
603 for (i=1; i<=v->ylist[0] && !miss; i++) {
604 vi = v->ylist[i];
605 if (na(dset->Z[vi][t])) {
606 v->t2 -= 1;
607 miss = 1;
608 }
609 }
610
611 if (v->xlist != NULL && !miss) {
612 for (i=1; i<=v->xlist[0] && !miss; i++) {
613 vi = v->xlist[i];
614 if (na(dset->Z[vi][t])) {
615 v->t2 -= 1;
616 miss = 1;
617 }
618 }
619 }
620
621 if (v->rlist != NULL && !miss) {
622 for (i=1; i<=v->rlist[0] && !miss; i++) {
623 vi = v->rlist[i];
624 if (na(dset->Z[vi][t])) {
625 v->t2 -= 1;
626 miss = 1;
627 }
628 }
629 }
630
631 if (!miss) {
632 break;
633 }
634 }
635
636 /* reject sample in case of internal missing values */
637
638 for (t=v->t1; t<=v->t2 && !err; t++) {
639
640 for (i=1; i<=v->ylist[0] && !err; i++) {
641 vi = v->ylist[i];
642 if (na(dset->Z[vi][t])) {
643 err = E_MISSDATA;
644 }
645 }
646
647 if (v->xlist != NULL && !err) {
648 for (i=1; i<=v->xlist[0] && !err; i++) {
649 vi = v->xlist[i];
650 if (na(dset->Z[vi][t])) {
651 err = E_MISSDATA;
652 }
653 }
654 }
655
656 if (v->rlist != NULL && !err) {
657 for (i=1; i<=v->rlist[0] && !err; i++) {
658 vi = v->rlist[i];
659 if (na(dset->Z[vi][t])) {
660 err = E_MISSDATA;
661 }
662 }
663 }
664 }
665
666 #if 0
667 fprintf(stderr, "VAR_set_sample: t1=%d, t2=%d\n", v->t1, v->t2);
668 #endif
669
670 return err;
671 }
672
673 /* Account for deterministic terms and check for
674 non-negative degrees of freedom.
675 */
676
VAR_check_df_etc(GRETL_VAR * v,const DATASET * dset,gretlopt opt)677 static int VAR_check_df_etc (GRETL_VAR *v, const DATASET *dset,
678 gretlopt opt)
679 {
680 int nl, err = 0;
681
682 v->T = v->t2 - v->t1 + 1;
683
684 nl = var_n_lags(v);
685 v->ncoeff = nl * v->neqns;
686
687 if (v->xlist != NULL) {
688 /* user-specified exogenous vars */
689 v->ncoeff += v->xlist[0];
690 }
691
692 if (v->ci == VECM) {
693 if (!(opt & OPT_N) && !(opt & OPT_R)) {
694 v->detflags |= DET_CONST;
695 v->ncoeff += 1;
696 v->ifc = 1;
697 }
698 } else if (!(opt & OPT_N)) {
699 v->detflags |= DET_CONST;
700 v->ncoeff += 1;
701 v->ifc = 1;
702 }
703
704 if ((opt & OPT_D) && dset->pd != 1) {
705 v->detflags |= DET_SEAS;
706 v->ncoeff += dset->pd - 1;
707 }
708
709 if (opt & OPT_T) {
710 v->detflags |= DET_TREND;
711 v->ncoeff += 1;
712 }
713
714 v->df = v->T - v->ncoeff;
715
716 if (v->df < 0) {
717 err = E_DF;
718 }
719
720 return err;
721 }
722
VAR_add_basic_matrices(GRETL_VAR * v)723 static int VAR_add_basic_matrices (GRETL_VAR *v)
724 {
725 int err = 0;
726
727 v->Y = gretl_matrix_alloc(v->T, v->neqns);
728 v->E = gretl_matrix_alloc(v->T, v->neqns);
729
730 if (v->Y == NULL || v->E == NULL) {
731 err = E_ALLOC;
732 }
733
734 if (!err && v->ncoeff > 0) {
735 v->X = gretl_matrix_alloc(v->T, v->ncoeff);
736 v->B = gretl_matrix_alloc(v->ncoeff, v->neqns);
737 if (v->X == NULL || v->B == NULL) {
738 err = E_ALLOC;
739 } else {
740 /* record the number of columns of X */
741 v->xcols = v->X->cols;
742 }
743 }
744
745 return err;
746 }
747
set_to_NA(double * x,int n)748 static void set_to_NA (double *x, int n)
749 {
750 int i;
751
752 for (i=0; i<n; i++) {
753 x[i] = NADBL;
754 }
755 }
756
757 /* main function for constructing a new VAR struct, which
758 may be used for estimating a VAR, a VECM, or various
759 auxiliary tasks */
760
gretl_VAR_new(int code,int order,int rank,const int * lags,const int * list,const DATASET * dset,gretlopt opt,int * errp)761 static GRETL_VAR *gretl_VAR_new (int code, int order, int rank,
762 const int *lags,
763 const int *list,
764 const DATASET *dset,
765 gretlopt opt, int *errp)
766 {
767 GRETL_VAR *var;
768 int ci, err = 0;
769
770 ci = (code >= VECM_ESTIMATE)? VECM : VAR;
771
772 if (order < 0) {
773 gretl_errmsg_sprintf(_("Invalid lag order %d"), order);
774 *errp = E_INVARG;
775 return NULL;
776 }
777
778 var = malloc(sizeof *var);
779 if (var == NULL) {
780 *errp = E_ALLOC;
781 return NULL;
782 }
783
784 gretl_VAR_clear(var);
785
786 var->ci = ci;
787 var->order = order;
788 var->lags = gretl_list_copy(lags);
789
790 if (ci == VAR) {
791 if (opt & OPT_H) {
792 var->robust = VAR_HAC;
793 } else if (opt & OPT_R) {
794 var->robust = VAR_HC;
795 }
796 }
797
798 err = VAR_make_lists(var, list, dset);
799
800 if (!err && rank > var->ylist[0]) {
801 gretl_errmsg_sprintf(_("vecm: rank %d is out of bounds"), rank);
802 err = E_DATA;
803 }
804
805 if (!err) {
806 var->neqns = var->ylist[0];
807 var->t1 = dset->t1;
808 var->t2 = dset->t2;
809 }
810
811 /* FIXME below: some of these allocations are redundant
812 for some uses of the VAR struct */
813
814 if (!err) {
815 err = VAR_set_sample(var, dset);
816 }
817
818 if (!err) {
819 err = VAR_check_df_etc(var, dset, opt);
820 }
821
822 if (!err) {
823 err = VAR_add_basic_matrices(var);
824 }
825
826 if (!err && var->ci == VAR) {
827 err = VAR_allocate_companion_matrix(var);
828 }
829
830 if (!err && var->ci == VAR) {
831 err = VAR_allocate_cholesky_matrix(var);
832 }
833
834 if (!err && code != VAR_LAGSEL) {
835 err = VAR_add_models(var, dset);
836 }
837
838 if (!err && code == VAR_ESTIMATE) {
839 int m = var->neqns * var->neqns + var->neqns;
840
841 var->Fvals = malloc(m * sizeof *var->Fvals);
842 if (var->Fvals == NULL) {
843 err = 1;
844 } else {
845 set_to_NA(var->Fvals, m);
846 }
847 }
848
849 if (!err && code == VAR_ESTIMATE) {
850 var->Ivals = malloc(N_IVALS * sizeof *var->Ivals);
851 if (var->Ivals == NULL) {
852 err = 1;
853 } else {
854 set_to_NA(var->Ivals, N_IVALS);
855 }
856 }
857
858 if (!err && var->ci == VECM) {
859 var->jinfo = johansen_info_new(var, rank, opt);
860 if (var->jinfo == NULL) {
861 err = E_ALLOC;
862 } else if (var->detflags & DET_SEAS) {
863 var->jinfo->seasonals = dset->pd - 1;
864 }
865 }
866
867 if (!err) {
868 if (var->ci == VAR) {
869 /* this is deferred for a VECM */
870 VAR_fill_Y(var, dset);
871 }
872 VAR_fill_X(var, order, dset);
873 }
874
875 if (err) {
876 gretl_VAR_free(var);
877 var = NULL;
878 }
879
880 *errp = err;
881
882 return var;
883 }
884
johansen_info_free(JohansenInfo * jv)885 static void johansen_info_free (JohansenInfo *jv)
886 {
887 gretl_matrix_free(jv->R0);
888 gretl_matrix_free(jv->R1);
889
890 gretl_matrix_free(jv->S00);
891 gretl_matrix_free(jv->S11);
892 gretl_matrix_free(jv->S01);
893
894 gretl_matrix_free(jv->evals);
895 gretl_matrix_free(jv->Beta);
896 gretl_matrix_free(jv->Alpha);
897 gretl_matrix_free(jv->Bvar);
898 gretl_matrix_free(jv->Bse);
899 gretl_matrix_free(jv->Ase);
900 gretl_matrix_free(jv->R);
901 gretl_matrix_free(jv->q);
902 gretl_matrix_free(jv->Ra);
903 gretl_matrix_free(jv->qa);
904
905 gretl_matrix_free(jv->YY);
906 gretl_matrix_free(jv->RR);
907 gretl_matrix_free(jv->BB);
908
909 free(jv);
910 }
911
gretl_VAR_free(GRETL_VAR * var)912 void gretl_VAR_free (GRETL_VAR *var)
913 {
914 if (var == NULL) return;
915
916 #if VDEBUG
917 fprintf(stderr, "gretl_VAR_free: var = %p, refcount = %d\n",
918 (void *) var, var->refcount);
919 #endif
920
921 var->refcount -= 1;
922 if (var->refcount > 0) {
923 return;
924 }
925
926 free(var->lags);
927 free(var->ylist);
928 free(var->xlist);
929 free(var->rlist);
930
931 gretl_matrix_free(var->Y);
932 gretl_matrix_free(var->X);
933 gretl_matrix_free(var->B);
934 gretl_matrix_free(var->XTX);
935
936 gretl_matrix_free(var->A);
937 gretl_matrix_free(var->L);
938 gretl_matrix_free(var->E);
939 gretl_matrix_free(var->C);
940 gretl_matrix_free(var->S);
941 gretl_matrix_free(var->F);
942 gretl_matrix_free(var->V);
943 gretl_matrix_free(var->ord);
944
945 free(var->Fvals);
946 free(var->Ivals);
947 free(var->name);
948
949 if (var->models != NULL) {
950 gretl_model_array_destroy(var->models, var->neqns);
951 }
952
953 if (var->jinfo != NULL) {
954 johansen_info_free(var->jinfo);
955 }
956
957 free(var);
958
959 #if VDEBUG
960 fprintf(stderr, "gretl_VAR_free: done\n");
961 fflush(stderr);
962 #endif
963 }
964
VAR_add_fcast_variance(GRETL_VAR * var,gretl_matrix * F,int n_static)965 static int VAR_add_fcast_variance (GRETL_VAR *var, gretl_matrix *F,
966 int n_static)
967 {
968 gretl_matrix *se = NULL;
969 double ftj, vti;
970 int k, n = var->neqns;
971 int i, j, s;
972 int err = 0;
973
974 if (n_static < 0) {
975 fprintf(stderr, "n_static = %d\n", n_static);
976 n_static = 0;
977 }
978
979 k = F->rows - n_static;
980
981 if (k > 0) {
982 se = gretl_VAR_get_fcast_se(var, k);
983 if (se == NULL) {
984 err = E_ALLOC;
985 goto bailout;
986 }
987 }
988
989 for (j=0; j<n; j++) {
990 for (s=0; s<F->rows; s++) {
991 ftj = gretl_matrix_get(F, s, j);
992 if (na(ftj)) {
993 gretl_matrix_set(F, s, n + j, NADBL);
994 } else {
995 i = s - n_static;
996 if (i < 0) {
997 vti = sqrt(gretl_matrix_get(var->S, j, j));
998 } else {
999 vti = gretl_matrix_get(se, i, j);
1000 }
1001 gretl_matrix_set(F, s, n + j, vti);
1002 }
1003 }
1004 }
1005
1006 if (se != NULL) {
1007 gretl_matrix_free(se);
1008 }
1009
1010 bailout:
1011
1012 if (err) {
1013 for (i=0; i<n; i++) {
1014 for (s=0; s<F->rows; s++) {
1015 gretl_matrix_set(F, s, n + i, NADBL);
1016 }
1017 }
1018 }
1019
1020 return err;
1021 }
1022
1023 /* determine start of dynamic portion of forecast */
1024
VAR_get_tdyn(GRETL_VAR * var,int t1,int t2,gretlopt opt)1025 static int VAR_get_tdyn (GRETL_VAR *var, int t1, int t2, gretlopt opt)
1026 {
1027 int td;
1028
1029 if (opt & OPT_D) {
1030 /* force dynamic */
1031 td = t1;
1032 } else if (opt & OPT_S) {
1033 /* force static */
1034 td = t2 + 1;
1035 } else {
1036 /* out of sample */
1037 td = var->t2 + 1;
1038 }
1039
1040 return td;
1041 }
1042
1043 static int
VAR_add_forecast(GRETL_VAR * var,int t1,int t2,const DATASET * dset,gretlopt opt)1044 VAR_add_forecast (GRETL_VAR *var, int t1, int t2,
1045 const DATASET *dset,
1046 gretlopt opt)
1047 {
1048 const MODEL *pmod;
1049 double fti, xti, xtid;
1050 int tdyn, nf = t2 - t1 + 1;
1051 int i, j, k, s, t, p;
1052 int lag, vj, m = 0;
1053 int fcols;
1054
1055 fcols = 2 * var->neqns;
1056
1057 /* rows = number of forecast periods; cols = 1 to hold forecast
1058 for each variable, plus 1 to hold variance for each variable
1059 if forecast is dynamic.
1060 */
1061
1062 var->F = gretl_zero_matrix_new(nf, fcols);
1063 if (var->F == NULL) {
1064 return E_ALLOC;
1065 }
1066
1067 if (var->detflags & DET_SEAS) {
1068 m = get_subperiod(t1, dset, NULL);
1069 }
1070
1071 /* start of dynamic portion of forecast? */
1072 tdyn = VAR_get_tdyn(var, t1, t2, opt);
1073
1074 #if VDEBUG
1075 fprintf(stderr, "var fcast: t1=%d, tdyn=%d, t2=%d\n", t1, tdyn, t2);
1076 #endif
1077
1078 for (t=t1, s=0; t<=t2; t++, s++) {
1079 int miss = 0;
1080
1081 for (i=0; i<var->neqns; i++) {
1082 pmod = var->models[i];
1083 fti = 0.0;
1084 k = 0;
1085
1086 /* constant, if present */
1087 if (pmod->ifc) {
1088 fti += pmod->coeff[k++];
1089 }
1090
1091 /* lags of stochastic vars */
1092 for (j=0; j<var->neqns; j++) {
1093 vj = var->ylist[j+1];
1094 for (lag=1; lag<=var->order; lag++) {
1095 if (!lag_wanted(var, lag)) {
1096 continue;
1097 }
1098 p = t - lag;
1099 xtid = NADBL;
1100 if (p < tdyn || s - lag < 0) {
1101 /* use actual data if possible */
1102 if (p < 0) {
1103 xti = NADBL;
1104 } else {
1105 xti = dset->Z[vj][p];
1106 }
1107 } else {
1108 /* prior forecast value preferred */
1109 if (p >= 0) {
1110 xtid = dset->Z[vj][p];
1111 }
1112 xti = gretl_matrix_get(var->F, s-lag, j);
1113 }
1114 if (!na(xti)) {
1115 fti += pmod->coeff[k] * xti;
1116 } else if (!na(xtid)) {
1117 fti += pmod->coeff[k] * xtid;
1118 } else {
1119 miss = 1;
1120 }
1121 k++;
1122 }
1123 }
1124
1125 /* exogenous vars, if any */
1126 if (!miss && var->xlist != NULL) {
1127 for (j=1; j<=var->xlist[0]; j++) {
1128 vj = var->xlist[j];
1129 xti = dset->Z[vj][t];
1130 if (na(xti)) {
1131 miss = 1;
1132 } else {
1133 fti += pmod->coeff[k] * xti;
1134 }
1135 k++;
1136 }
1137 }
1138
1139 /* other deterministics, if present */
1140 if (!miss) {
1141 if (var->detflags & DET_SEAS) {
1142 if (m < dset->pd - 1) {
1143 fti += pmod->coeff[k+m];
1144 }
1145 }
1146 if (var->detflags & DET_TREND) {
1147 k = pmod->ncoeff - 1;
1148 fti += pmod->coeff[k] * (t + 1);
1149 }
1150 }
1151
1152 if (miss) {
1153 fti = NADBL;
1154 }
1155
1156 gretl_matrix_set(var->F, s, i, fti);
1157 }
1158
1159 if (var->detflags & DET_SEAS) {
1160 m = (m < dset->pd - 1)? m + 1 : 0;
1161 }
1162 }
1163
1164 VAR_add_fcast_variance(var, var->F, tdyn - t1);
1165
1166 if (var->F != NULL) {
1167 gretl_matrix_set_t1(var->F, t1);
1168 gretl_matrix_set_t2(var->F, t2);
1169 }
1170
1171 #if VDEBUG
1172 gretl_matrix_print(var->F, "var->F");
1173 #endif
1174
1175 return 0;
1176 }
1177
1178 static int
VECM_add_forecast(GRETL_VAR * var,int t1,int t2,const DATASET * dset,gretlopt opt)1179 VECM_add_forecast (GRETL_VAR *var, int t1, int t2,
1180 const DATASET *dset, gretlopt opt)
1181 {
1182 gretl_matrix *B = NULL;
1183 double s0 = 0, s1 = 1;
1184 int order = effective_order(var);
1185 int nexo = (var->xlist != NULL)? var->xlist[0] : 0;
1186 int nseas = var->jinfo->seasonals;
1187 int tdyn, nf = t2 - t1 + 1;
1188 int i, j, k, vj, s, t;
1189 int fcols, m = 0;
1190
1191 fcols = 2 * var->neqns;
1192
1193 var->F = gretl_zero_matrix_new(nf, fcols);
1194 if (var->F == NULL) {
1195 return E_ALLOC;
1196 }
1197
1198 B = VAR_coeff_matrix_from_VECM(var, 0);
1199 if (B == NULL) {
1200 gretl_matrix_free(var->F);
1201 var->F = NULL;
1202 return E_ALLOC;
1203 }
1204
1205 /* start of dynamic portion of forecast? */
1206 tdyn = VAR_get_tdyn(var, t1, t2, opt);
1207
1208 if (nseas > 0) {
1209 m = get_subperiod(t1, dset, NULL);
1210 s1 -= 1.0 / dset->pd;
1211 s0 = s1 - 1;
1212 }
1213
1214 for (t=t1, s=0; t<=t2; t++, s++) {
1215
1216 for (i=0; i<var->neqns; i++) {
1217 double bij, xtj, xtjd, fti = 0.0;
1218 int col = 0;
1219
1220 /* unrestricted constant, if present */
1221 if (var->ifc) {
1222 bij = gretl_matrix_get(B, i, col++);
1223 fti += bij;
1224 }
1225
1226 /* lags of endogenous vars */
1227 for (j=0; j<var->neqns; j++) {
1228 vj = var->ylist[j+1];
1229 for (k=1; k<=order; k++) {
1230 /* FIXME gappy lags */
1231 if (t - k < 0) {
1232 fti = NADBL;
1233 break;
1234 }
1235 bij = gretl_matrix_get(B, i, col++);
1236 xtjd = NADBL;
1237 if (t < tdyn || s - k < 0) {
1238 /* pre-forecast value */
1239 xtj = dset->Z[vj][t-k];
1240 } else {
1241 /* prior forecast value preferred */
1242 xtjd = dset->Z[vj][t-k];
1243 xtj = gretl_matrix_get(var->F, s-k, j);
1244 }
1245 if (!na(xtj)) {
1246 fti += bij * xtj;
1247 } else if (!na(xtjd)) {
1248 fti += bij * xtjd;
1249 } else {
1250 fti = NADBL;
1251 break;
1252 }
1253 }
1254 if (na(fti)) {
1255 break;
1256 }
1257 }
1258
1259 if (na(fti)) {
1260 goto set_fcast;
1261 }
1262
1263 /* exogenous vars, if present */
1264 for (j=0; j<nexo; j++) {
1265 vj = var->xlist[j+1];
1266 xtj = dset->Z[vj][t];
1267 if (na(xtj)) {
1268 fti = NADBL;
1269 } else {
1270 bij = gretl_matrix_get(B, i, col++);
1271 fti += bij * xtj;
1272 }
1273 }
1274
1275 if (na(fti)) {
1276 goto set_fcast;
1277 }
1278
1279 /* seasonals, if present */
1280 for (j=0; j<nseas; j++) {
1281 xtj = (m == j)? s1 : s0;
1282 bij = gretl_matrix_get(B, i, col++);
1283 fti += bij * xtj;
1284 }
1285
1286 if (jcode(var) == J_UNREST_TREND) {
1287 /* unrestricted trend */
1288 bij = gretl_matrix_get(B, i, col++);
1289 fti += bij * (t + 1);
1290 } else if (jcode(var) == J_REST_CONST) {
1291 /* restricted constant */
1292 fti += gretl_matrix_get(B, i, col++);
1293 } else if (jcode(var) == J_REST_TREND) {
1294 /* restricted trend */
1295 bij = gretl_matrix_get(B, i, col++);
1296 fti += bij * t;
1297 }
1298
1299 /* restricted exog vars */
1300 if (var->rlist != NULL) {
1301 for (j=0; j<var->rlist[0]; j++) {
1302 vj = var->rlist[j+1];
1303 xtj = dset->Z[vj][t-1]; /* ?? */
1304 if (na(xtj)) {
1305 fti = NADBL;
1306 } else {
1307 bij = gretl_matrix_get(B, i, col++);
1308 fti += bij * xtj;
1309 }
1310 }
1311 }
1312
1313 set_fcast:
1314 gretl_matrix_set(var->F, s, i, fti);
1315 }
1316
1317 if (nseas > 0) {
1318 m = (m < dset->pd - 1)? m + 1 : 0;
1319 }
1320 }
1321
1322 gretl_matrix_free(B);
1323
1324 VAR_add_fcast_variance(var, var->F, tdyn - t1);
1325
1326 gretl_matrix_set_t1(var->F, t1);
1327 gretl_matrix_set_t2(var->F, t2);
1328
1329 #if VDEBUG
1330 gretl_matrix_print(var->F, "var->F");
1331 #endif
1332
1333 return 0;
1334 }
1335
1336 const gretl_matrix *
gretl_VAR_get_forecast_matrix(GRETL_VAR * var,int t1,int t2,DATASET * dset,gretlopt opt,int * err)1337 gretl_VAR_get_forecast_matrix (GRETL_VAR *var, int t1, int t2,
1338 DATASET *dset, gretlopt opt,
1339 int *err)
1340 {
1341 if (var->F != NULL) {
1342 /* there's a forecast attached, but it may not be what we want */
1343 gretl_matrix_free(var->F);
1344 var->F = NULL;
1345 }
1346
1347 if (var->ci == VECM) {
1348 *err = VECM_add_forecast(var, t1, t2, dset, opt);
1349 } else {
1350 *err = VAR_add_forecast(var, t1, t2, dset, opt);
1351 }
1352
1353 return var->F;
1354 }
1355
VAR_dw_rho(MODEL * pmod)1356 static void VAR_dw_rho (MODEL *pmod)
1357 {
1358 double ut, u1;
1359 double xx = 0;
1360 double ut1 = 0, u11 = 0;
1361 int t, s;
1362
1363 for (t=pmod->t1; t<=pmod->t2; t++) {
1364 s = t - 1;
1365 if (s >= 0 && !na(pmod->uhat[s])) {
1366 ut = pmod->uhat[t];
1367 u1 = pmod->uhat[s];
1368 xx += (ut - u1) * (ut - u1);
1369 ut1 += ut * u1;
1370 u11 += u1 * u1;
1371 }
1372 }
1373
1374 pmod->dw = xx / pmod->ess;
1375 pmod->rho = ut1 / u11;
1376 }
1377
1378 /* set basic per-model statistics after estimation has been
1379 done by matrix methods */
1380
set_VAR_model_stats(GRETL_VAR * var,int i)1381 int set_VAR_model_stats (GRETL_VAR *var, int i)
1382 {
1383 MODEL *pmod = var->models[i];
1384 const double *y;
1385 double u, x, SSR = 0, TSS = 0;
1386 int t;
1387
1388 /* get the appropriate column of var->Y as a plain array */
1389 y = var->Y->val + i * var->T;
1390
1391 pmod->ybar = gretl_mean(0, var->T - 1, y);
1392 pmod->sdy = gretl_stddev(0, var->T - 1, y);
1393
1394 for (t=0; t<var->T; t++) {
1395 u = gretl_matrix_get(var->E, t, i);
1396 SSR += u * u;
1397 x = (var->ifc)? (y[t] - pmod->ybar) : y[t];
1398 TSS += x * x;
1399 pmod->uhat[t + pmod->t1] = u;
1400 pmod->yhat[t + pmod->t1] = y[t] - u;
1401 }
1402
1403 pmod->ess = SSR;
1404 #if VAR_SE_DFCORR
1405 pmod->sigma = sqrt(SSR / pmod->dfd);
1406 #else
1407 pmod->sigma = sqrt(SSR / var->T);
1408 #endif
1409 pmod->tss = TSS;
1410 pmod->rsq = 1.0 - SSR / TSS;
1411 pmod->adjrsq = 1.0 - (SSR / (pmod->dfd)) / (TSS / (pmod->nobs - 1));
1412 pmod->fstt = ((TSS - SSR) / pmod->dfn) / (SSR / pmod->dfd);
1413
1414 pmod->lnL = NADBL;
1415
1416 VAR_dw_rho(pmod);
1417
1418 return 0;
1419 }
1420
gretl_VAR_do_error_decomp(const gretl_matrix * S,gretl_matrix * C,const gretl_matrix * ord)1421 int gretl_VAR_do_error_decomp (const gretl_matrix *S,
1422 gretl_matrix *C,
1423 const gretl_matrix *ord)
1424 {
1425 int g = gretl_matrix_rows(S);
1426 gretl_matrix *tmp = NULL;
1427 double x;
1428 int i, j, r, c;
1429 int err = 0;
1430
1431 /* copy cross-equation covariance matrix */
1432 tmp = gretl_matrix_copy(S);
1433 if (tmp == NULL) {
1434 err = E_ALLOC;
1435 }
1436
1437 if (ord != NULL) {
1438 for (i=0; i<g; i++) {
1439 r = ord->val[i];
1440 for (j=0; j<g; j++) {
1441 c = ord->val[j];
1442 x = gretl_matrix_get(S, r, c);
1443 gretl_matrix_set(tmp, i, j, x);
1444 }
1445 }
1446 }
1447
1448 /* lower-triangularize and decompose */
1449 if (!err) {
1450 for (i=0; i<g-1; i++) {
1451 for (j=i+1; j<g; j++) {
1452 gretl_matrix_set(tmp, i, j, 0.0);
1453 }
1454 }
1455 err = gretl_matrix_cholesky_decomp(tmp);
1456 }
1457
1458 /* write the decomposition (lower triangle) into the C matrix */
1459 if (!err) {
1460 for (i=0; i<g; i++) {
1461 r = (ord != NULL)? ord->val[i] : i;
1462 for (j=0; j<=i; j++) {
1463 c = (ord != NULL)? ord->val[j] : j;
1464 x = gretl_matrix_get(tmp, i, j);
1465 gretl_matrix_set(C, r, c, x);
1466 }
1467 }
1468 }
1469
1470 if (tmp != NULL) {
1471 gretl_matrix_free(tmp);
1472 }
1473
1474 return err;
1475 }
1476
1477 const gretl_matrix *
gretl_VAR_get_residual_matrix(const GRETL_VAR * var)1478 gretl_VAR_get_residual_matrix (const GRETL_VAR *var)
1479 {
1480 if (var->E != NULL) {
1481 gretl_matrix_set_t1(var->E, var->t1);
1482 gretl_matrix_set_t2(var->E, var->t2);
1483 }
1484
1485 return var->E;
1486 }
1487
1488 #define samesize(p,q) (p->rows == q->rows && p->cols == q->cols)
1489
1490 gretl_matrix *
gretl_VAR_get_fitted_matrix(const GRETL_VAR * var)1491 gretl_VAR_get_fitted_matrix (const GRETL_VAR *var)
1492 {
1493 gretl_matrix *Yh = NULL;
1494
1495 if (var->Y != NULL && var->E != NULL && samesize(var->Y, var->E)) {
1496 Yh = gretl_matrix_copy(var->Y);
1497 if (Yh != NULL) {
1498 gretl_matrix_subtract_from(Yh, var->E);
1499 gretl_matrix_set_t1(Yh, var->t1);
1500 gretl_matrix_set_t2(Yh, var->t2);
1501 }
1502 }
1503
1504 return Yh;
1505 }
1506
1507 gretl_matrix *
gretl_VAR_get_vma_matrix(const GRETL_VAR * var,const DATASET * dset,int * err)1508 gretl_VAR_get_vma_matrix (const GRETL_VAR *var, const DATASET *dset,
1509 int *err)
1510 {
1511 int h = default_VAR_horizon(dset);
1512 int ar, n = var->neqns;
1513 int n2 = n * n;
1514 gretl_matrix *VMA = NULL;
1515 gretl_matrix *Tmp1, *Tmp2;
1516 double x;
1517 int i, j, ii, jj;
1518
1519 if (var->A == NULL) {
1520 /* companion matrix is absent */
1521 *err = E_BADSTAT;
1522 return NULL;
1523 }
1524
1525 ar = var->A->rows;
1526
1527 Tmp1 = gretl_identity_matrix_new(ar);
1528 Tmp2 = gretl_matrix_alloc(ar, ar);
1529 if (Tmp1 == NULL || Tmp2 == NULL) {
1530 *err = E_ALLOC;
1531 goto bailout;
1532 }
1533
1534 VMA = gretl_zero_matrix_new(h, n2);
1535 if (VMA == NULL) {
1536 *err = E_ALLOC;
1537 goto bailout;
1538 }
1539
1540 /* compose first row of VMA = vec(I(n))' */
1541 for (i=0; i<n2; i+=n+1) {
1542 gretl_matrix_set(VMA, 0, i, 1.0);
1543 }
1544
1545 for (i=1; i<h; i++) {
1546 /* Tmp <- A * Tmp */
1547 gretl_matrix_multiply(var->A, Tmp1, Tmp2);
1548 gretl_matrix_copy_values(Tmp1, Tmp2);
1549 /* VMA |= vec(Tmp[1:n,1:n])' */
1550 ii = jj = 0;
1551 for (j=0; j<n2; j++) {
1552 x = gretl_matrix_get(Tmp1, ii++, jj);
1553 gretl_matrix_set(VMA, i, j, x);
1554 if (ii == n) {
1555 jj++;
1556 ii = 0;
1557 }
1558 }
1559 }
1560
1561 bailout:
1562
1563 gretl_matrix_free(Tmp1);
1564 gretl_matrix_free(Tmp2);
1565
1566 return VMA;
1567 }
1568
1569 static gretl_matrix *
gretl_VAR_get_portmanteau_test(const GRETL_VAR * var,int * err)1570 gretl_VAR_get_portmanteau_test (const GRETL_VAR *var, int *err)
1571 {
1572 gretl_matrix *m = NULL;
1573
1574 if (na(var->LB)) {
1575 *err = E_BADSTAT;
1576 } else {
1577 m = gretl_matrix_alloc(1, 4);
1578 if (m == NULL) {
1579 *err = E_ALLOC;
1580 } else {
1581 int k = var->order + (var->ci == VECM);
1582 int df = var->neqns * var->neqns * (var->LBs - k);
1583 double pv = chisq_cdf_comp(df, var->LB);
1584
1585 m->val[0] = var->LB;
1586 m->val[1] = var->LBs;
1587 m->val[2] = df;
1588 m->val[3] = pv;
1589 *err = 0;
1590 }
1591 }
1592
1593 return m;
1594 }
1595
gretl_VAR_get_variable_number(const GRETL_VAR * var,int k)1596 int gretl_VAR_get_variable_number (const GRETL_VAR *var, int k)
1597 {
1598 int vnum = 0;
1599
1600 if (var->models != NULL && k >= 0 && k < var->neqns) {
1601 if (var->models[k] != NULL && var->models[k]->list != NULL) {
1602 vnum = var->models[k]->list[1];
1603 }
1604 }
1605
1606 return vnum;
1607 }
1608
gretl_VAR_get_n_equations(const GRETL_VAR * var)1609 int gretl_VAR_get_n_equations (const GRETL_VAR *var)
1610 {
1611 return (var == NULL)? 0 : var->neqns;
1612 }
1613
gretl_VAR_get_t1(const GRETL_VAR * var)1614 int gretl_VAR_get_t1 (const GRETL_VAR *var)
1615 {
1616 return (var == NULL)? 0 : var->t1;
1617 }
1618
gretl_VAR_get_t2(const GRETL_VAR * var)1619 int gretl_VAR_get_t2 (const GRETL_VAR *var)
1620 {
1621 return (var == NULL)? 0 : var->t2;
1622 }
1623
gretl_var_get_sample(const GRETL_VAR * var,int * t1,int * t2)1624 int gretl_var_get_sample (const GRETL_VAR *var, int *t1, int *t2)
1625 {
1626 if (var != NULL && var->models != NULL && var->models[0] != NULL) {
1627 MODEL *pmod = var->models[0];
1628
1629 if (pmod->smpl.t1 >= 0 && pmod->smpl.t2 > pmod->smpl.t1) {
1630 *t1 = pmod->smpl.t1;
1631 *t2 = pmod->smpl.t2;
1632 return 0;
1633 }
1634 }
1635
1636 return E_DATA;
1637 }
1638
gretl_VAR_get_model(const GRETL_VAR * var,int i)1639 const MODEL *gretl_VAR_get_model (const GRETL_VAR *var, int i)
1640 {
1641 if (var != NULL && var->models != NULL && i < var->neqns) {
1642 return var->models[i];
1643 } else {
1644 return NULL;
1645 }
1646 }
1647
periods_from_pd(int pd)1648 static int periods_from_pd (int pd)
1649 {
1650 int periods = 10;
1651
1652 if (pd == 4) {
1653 /* quarterly: try 5 years */
1654 periods = 20;
1655 } else if (pd == 12) {
1656 /* monthly: two years */
1657 periods = 24;
1658 } else if (pd == 7 || pd == 6 || pd == 5) {
1659 /* daily: three weeks */
1660 periods = 3 * pd;
1661 }
1662
1663 return periods;
1664 }
1665
default_VAR_horizon(const DATASET * dset)1666 int default_VAR_horizon (const DATASET *dset)
1667 {
1668 int h = libset_get_int(HORIZON);
1669
1670 if (h <= 0) {
1671 if (dset != NULL) {
1672 h = periods_from_pd(dset->pd);
1673 } else {
1674 h = 20; /* default */
1675 }
1676 }
1677
1678 return h;
1679 }
1680
reorder_responses(const GRETL_VAR * var,int * err)1681 gretl_matrix *reorder_responses (const GRETL_VAR *var, int *err)
1682 {
1683 gretl_matrix *S, *C;
1684 int i, j, r, c;
1685 double x;
1686
1687 S = gretl_matrix_copy(var->S);
1688 C = gretl_matrix_alloc(var->neqns, var->neqns);
1689
1690 if (S == NULL || C == NULL) {
1691 gretl_matrix_free(S);
1692 gretl_matrix_free(C);
1693 *err = E_ALLOC;
1694 return NULL;
1695 }
1696
1697 /* rearrange @S */
1698 for (i=0; i<var->neqns; i++) {
1699 r = var->ord->val[i];
1700 for (j=0; j<var->neqns; j++) {
1701 c = var->ord->val[j];
1702 x = gretl_matrix_get(var->S, r, c);
1703 gretl_matrix_set(S, i, j, x);
1704 }
1705 }
1706
1707 /* decompose rearranged @S */
1708 gretl_matrix_cholesky_decomp(S);
1709
1710 /* C = rearranged decomposed S */
1711 for (i=0; i<var->neqns; i++) {
1712 r = var->ord->val[i];
1713 for (j=0; j<var->neqns; j++) {
1714 c = var->ord->val[j];
1715 x = gretl_matrix_get(S, i, j);
1716 gretl_matrix_set(C, r, c, x);
1717 }
1718 }
1719
1720 gretl_matrix_free(S);
1721
1722 return C;
1723 }
1724
copy_north_west(gretl_matrix * targ,const gretl_matrix * src,int add)1725 void copy_north_west (gretl_matrix *targ,
1726 const gretl_matrix *src,
1727 int add)
1728 {
1729 int i, j, k = src->cols;
1730 double xij;
1731
1732 if (!add) {
1733 gretl_matrix_zero(targ);
1734 }
1735
1736 for (j=0; j<k; j++) {
1737 for (i=0; i<k; i++) {
1738 xij = gretl_matrix_get(src, i, j);
1739 if (add) {
1740 xij += gretl_matrix_get(targ, i, j);
1741 }
1742 gretl_matrix_set(targ, i, j, xij);
1743 }
1744 }
1745 }
1746
real_point_responses(const gretl_matrix * A,const gretl_matrix * C,gretl_matrix * resp,int targ,int shock)1747 static int real_point_responses (const gretl_matrix *A,
1748 const gretl_matrix *C,
1749 gretl_matrix *resp,
1750 int targ, int shock)
1751 {
1752 gretl_matrix *rtmp;
1753 gretl_matrix *ctmp;
1754 double rij;
1755 int n = C->rows;
1756 int np = A->cols;
1757 int i, j, k, t;
1758
1759 /* workspace */
1760 rtmp = gretl_matrix_alloc(np, n);
1761 ctmp = gretl_zero_matrix_new(np, n);
1762 if (rtmp == NULL || ctmp == NULL) {
1763 gretl_matrix_free(rtmp);
1764 gretl_matrix_free(ctmp);
1765 return E_ALLOC;
1766 }
1767
1768 for (t=0; t<resp->rows; t++) {
1769 if (t == 0) {
1770 /* initial estimated responses */
1771 copy_north_west(rtmp, C, 0);
1772 } else {
1773 /* calculate further estimated responses */
1774 gretl_matrix_multiply(A, rtmp, ctmp);
1775 gretl_matrix_copy_values(rtmp, ctmp);
1776 }
1777 if (resp->cols == 1) {
1778 resp->val[t] = gretl_matrix_get(rtmp, targ, shock);
1779 } else if (shock >= 0) {
1780 /* all targets for one shock */
1781 for (i=0; i<n; i++) {
1782 rij = gretl_matrix_get(rtmp, i, shock);
1783 gretl_matrix_set(resp, t, i, rij);
1784 }
1785 } else if (targ >= 0) {
1786 /* all shocks for one target */
1787 for (j=0; j<n; j++) {
1788 rij = gretl_matrix_get(rtmp, targ, j);
1789 gretl_matrix_set(resp, t, j, rij);
1790 }
1791 } else {
1792 /* all shocks, all targets */
1793 k = 0;
1794 for (i=0; i<n; i++) {
1795 for (j=0; j<n; j++) {
1796 rij = gretl_matrix_get(rtmp, i, j);
1797 gretl_matrix_set(resp, t, k++, rij);
1798 }
1799 }
1800 }
1801 }
1802
1803 gretl_matrix_free(rtmp);
1804 gretl_matrix_free(ctmp);
1805
1806 return 0;
1807 }
1808
1809 /* Calculate the VMA representation based on the A and C matrices,
1810 for a given horizon.
1811 */
1812
vma_rep(gretl_matrix * A,gretl_matrix * C,int horizon,int * err)1813 gretl_matrix *vma_rep (gretl_matrix *A, gretl_matrix *C,
1814 int horizon, int *err)
1815 {
1816 gretl_matrix *resp = NULL;
1817 gretl_matrix *realC = C;
1818 gretl_matrix *realA = A;
1819 int neqns, dim;
1820
1821 if (horizon <= 0) {
1822 *err = E_INVARG;
1823 return NULL;
1824 }
1825
1826 neqns = A->rows;
1827 dim = A->cols;
1828
1829 if (realC == NULL) {
1830 /* manufacture "plain" @C for convenience */
1831 realC = gretl_identity_matrix_new(neqns);
1832 if (realC == NULL) {
1833 *err = E_ALLOC;
1834 }
1835 }
1836
1837 if (!*err && dim > neqns) {
1838 /* the incoming @A is just the top @neqns rows
1839 of the companion matrix
1840 */
1841 realA = companionize(A, err);
1842 }
1843
1844 if (!*err) {
1845 int nresp = neqns * neqns;
1846
1847 resp = gretl_matrix_alloc(horizon, nresp);
1848
1849 if (resp == NULL) {
1850 *err = E_ALLOC;
1851 } else {
1852 if (neqns == 1) {
1853 /* "all" means "first" in this case */
1854 *err = real_point_responses(realA, realC, resp, 0, 0);
1855 } else {
1856 *err = real_point_responses(realA, realC, resp, -1, -1);
1857 }
1858 }
1859 }
1860
1861 if (realC != C) {
1862 gretl_matrix_free(realC);
1863 }
1864 if (realA != A) {
1865 gretl_matrix_free(realA);
1866 }
1867 if (*err && resp != NULL) {
1868 gretl_matrix_free(resp);
1869 resp = NULL;
1870 }
1871
1872 return resp;
1873 }
1874
1875 static gretl_matrix *
gretl_VAR_get_point_responses(GRETL_VAR * var,int targ,int shock,int periods,int * err)1876 gretl_VAR_get_point_responses (GRETL_VAR *var, int targ, int shock,
1877 int periods, int *err)
1878 {
1879 gretl_matrix *resp = NULL;
1880 gretl_matrix *C = var->C;
1881 int nresp = 1;
1882
1883 if (shock >= var->neqns) {
1884 fprintf(stderr, "Shock variable out of bounds\n");
1885 *err = E_INVARG;
1886 }
1887 if (!*err && targ >= var->neqns) {
1888 fprintf(stderr, "Target variable out of bounds\n");
1889 *err = E_INVARG;
1890 }
1891 if (!*err && periods <= 0) {
1892 fprintf(stderr, "Invalid number of periods\n");
1893 *err = E_INVARG;
1894 }
1895
1896 if (!*err && var->ord != NULL) {
1897 C = reorder_responses(var, err);
1898 }
1899
1900 if (*err) {
1901 return NULL;
1902 }
1903
1904 if (targ < 0 && shock < 0) {
1905 /* compute all */
1906 nresp = var->neqns * var->neqns;
1907 } else if (targ < 0 || shock < 0) {
1908 nresp = var->neqns;
1909 }
1910
1911 resp = gretl_matrix_alloc(periods, nresp);
1912
1913 if (resp == NULL) {
1914 *err = E_ALLOC;
1915 } else {
1916 *err = real_point_responses(var->A, C, resp, targ, shock);
1917 }
1918
1919 if (C != var->C) {
1920 gretl_matrix_free(C);
1921 }
1922
1923 if (*err && resp != NULL) {
1924 gretl_matrix_free(resp);
1925 resp = NULL;
1926 }
1927
1928 return resp;
1929 }
1930
1931 /**
1932 * gretl_VAR_get_impulse_response:
1933 * @var: pointer to VAR struct.
1934 * @targ: index of the target or response variable.
1935 * @shock: index of the source or shock variable.
1936 * @periods: number of periods over which to compute the response.
1937 * @alpha: determines confidence level for bootstrap interval.
1938 * @dset: dataset struct.
1939 * @err: location to receive error code.
1940 *
1941 * Computes the response of @targ to a perturbation of @shock
1942 * in the context of @var: @targ and @shock are zero-based indices
1943 * relative to the structure of @var. For example if @targ = 0 and
1944 * @shock = 1, we compute the response of the dependent variable in
1945 * the first VAR equation to a perturbation of the variable that
1946 * appears as dependent in the second VAR equation.
1947 *
1948 * If @alpha is zero, the response matrix returned is a column vector
1949 * of length @periods, giving the point estimate of the response
1950 * function. Otherwise, the response matrix returned has
1951 * three columns, containing the point estimate, the @alpha / 2
1952 * quantile and the 1 - @alpha / 2 quantile, where the quantiles
1953 * are based on bootstrap replication, with resampling of the
1954 * original residuals with replacement.
1955 *
1956 * Returns: matrix containing the estimated impulse responses,
1957 * with or without a confidence interval.
1958 */
1959
1960 gretl_matrix *
gretl_VAR_get_impulse_response(GRETL_VAR * var,int targ,int shock,int periods,double alpha,const DATASET * dset,int * err)1961 gretl_VAR_get_impulse_response (GRETL_VAR *var,
1962 int targ, int shock,
1963 int periods, double alpha,
1964 const DATASET *dset,
1965 int *err)
1966 {
1967 gretl_matrix *point = NULL;
1968 gretl_matrix *full = NULL;
1969 gretl_matrix *ret = NULL;
1970
1971 if (periods == 0) {
1972 if (dset != NULL) {
1973 periods = default_VAR_horizon(dset);
1974 } else {
1975 *err = E_DATA;
1976 return NULL;
1977 }
1978 }
1979
1980 if (alpha != 0 && (alpha < 0.01 || alpha > 0.6)) {
1981 *err = E_DATA;
1982 }
1983
1984 point = gretl_VAR_get_point_responses(var, targ, shock, periods, err);
1985
1986 if (dset == NULL || dset->Z == NULL || alpha == 0.0) {
1987 /* no data matrix provided, or no alpha given:
1988 just return point estimate */
1989 ret = point;
1990 } else if (point != NULL) {
1991 full = irf_bootstrap(var, targ, shock, periods,
1992 alpha, point, dset, err);
1993 gretl_matrix_free(point);
1994 ret = full;
1995 }
1996
1997 return ret;
1998 }
1999
2000 static gretl_matrix *
gretl_VAR_get_fcast_se(GRETL_VAR * var,int periods)2001 gretl_VAR_get_fcast_se (GRETL_VAR *var, int periods)
2002 {
2003 int k = var->neqns * effective_order(var);
2004 gretl_matrix *vtmp = NULL;
2005 gretl_matrix *vt = NULL;
2006 gretl_matrix *se = NULL;
2007 double vti;
2008 int i, t;
2009
2010 if (periods <= 0) {
2011 fprintf(stderr, "Invalid number of periods\n");
2012 return NULL;
2013 }
2014
2015 se = gretl_zero_matrix_new(periods, var->neqns);
2016 vt = gretl_matrix_alloc(k, k);
2017 vtmp = gretl_matrix_alloc(k, k);
2018
2019 if (se == NULL || vt == NULL || vtmp == NULL) {
2020 gretl_matrix_free(se);
2021 gretl_matrix_free(vt);
2022 gretl_matrix_free(vtmp);
2023 return NULL;
2024 }
2025
2026 for (t=0; t<periods; t++) {
2027 if (t == 0) {
2028 /* initial variance */
2029 copy_north_west(vt, var->S, 0);
2030 } else {
2031 /* calculate further variances */
2032 gretl_matrix_copy_values(vtmp, vt);
2033 gretl_matrix_qform(var->A, GRETL_MOD_NONE,
2034 vtmp, vt, GRETL_MOD_NONE);
2035 copy_north_west(vt, var->S, 1);
2036 }
2037 for (i=0; i<var->neqns; i++) {
2038 vti = gretl_matrix_get(vt, i, i);
2039 gretl_matrix_set(se, t, i, sqrt(vti));
2040 }
2041 }
2042
2043 gretl_matrix_free(vt);
2044 gretl_matrix_free(vtmp);
2045
2046 return se;
2047 }
2048
2049 gretl_matrix *
gretl_VAR_get_fcast_decomp(const GRETL_VAR * var,int targ,int periods,int * err)2050 gretl_VAR_get_fcast_decomp (const GRETL_VAR *var,
2051 int targ, int periods,
2052 int *err)
2053 {
2054 int n = var->neqns;
2055 int k = n * effective_order(var);
2056 gretl_matrix_block *B;
2057 gretl_matrix *idx, *vtmp;
2058 gretl_matrix *cic, *vt;
2059 gretl_matrix *vd = NULL;
2060 gretl_matrix *C = var->C;
2061 int i, t;
2062
2063 *err = 0;
2064
2065 if (targ >= n) {
2066 fprintf(stderr, "Target variable out of bounds\n");
2067 *err = E_DATA;
2068 }
2069
2070 if (!*err && periods <= 0) {
2071 fprintf(stderr, "Invalid number of periods\n");
2072 *err = E_DATA;
2073 }
2074
2075 if (!*err && var->ord != NULL) {
2076 C = reorder_responses(var, err);
2077 }
2078
2079 if (*err) {
2080 return NULL;
2081 }
2082
2083 B = gretl_matrix_block_new(&idx, n, n,
2084 &cic, n, n,
2085 &vt, k, k,
2086 &vtmp, k, k,
2087 NULL);
2088 if (B == NULL) {
2089 *err = E_ALLOC;
2090 goto bailout;
2091 }
2092
2093 vd = gretl_zero_matrix_new(periods, n + 1);
2094 if (vd == NULL) {
2095 *err = E_ALLOC;
2096 goto bailout;
2097 }
2098
2099 gretl_matrix_zero(idx);
2100
2101 for (i=0; i<n && !*err; i++) {
2102 double vti;
2103
2104 /* adjust index matrix */
2105 if (i > 0) {
2106 gretl_matrix_set(idx, i-1, i-1, 0.0);
2107 }
2108 gretl_matrix_set(idx, i, i, 1.0);
2109
2110 for (t=0; t<periods && !*err; t++) {
2111 if (t == 0) {
2112 /* calculate initial variances */
2113 *err = gretl_matrix_qform(C, GRETL_MOD_NONE,
2114 idx, cic, GRETL_MOD_NONE);
2115 copy_north_west(vt, cic, 0);
2116 } else {
2117 /* calculate further variances */
2118 gretl_matrix_copy_values(vtmp, vt);
2119 *err = gretl_matrix_qform(var->A, GRETL_MOD_NONE,
2120 vtmp, vt, GRETL_MOD_NONE);
2121 copy_north_west(vt, cic, 1);
2122 }
2123 if (!*err) {
2124 vti = gretl_matrix_get(vt, targ, targ);
2125 gretl_matrix_set(vd, t, i, vti);
2126 }
2127 }
2128 }
2129
2130 for (t=0; t<periods && !*err; t++) {
2131 double vi, vtot = 0.0;
2132
2133 for (i=0; i<n; i++) {
2134 vtot += gretl_matrix_get(vd, t, i);
2135 }
2136 /* normalize variance contributions as % shares */
2137 for (i=0; i<n; i++) {
2138 vi = gretl_matrix_get(vd, t, i);
2139 gretl_matrix_set(vd, t, i, 100.0 * vi / vtot);
2140 }
2141
2142 gretl_matrix_set(vd, t, var->neqns, sqrt(vtot));
2143 }
2144
2145 bailout:
2146
2147 gretl_matrix_block_destroy(B);
2148
2149 if (C != var->C) {
2150 gretl_matrix_free(C);
2151 }
2152
2153 if (*err && vd != NULL) {
2154 gretl_matrix_free(vd);
2155 vd = NULL;
2156 }
2157
2158 return vd;
2159 }
2160
2161 gretl_matrix *
gretl_VAR_get_FEVD_matrix(const GRETL_VAR * var,int targ,int shock,int horizon,const DATASET * dset,int * err)2162 gretl_VAR_get_FEVD_matrix (const GRETL_VAR *var,
2163 int targ, int shock,
2164 int horizon,
2165 const DATASET *dset,
2166 int *err)
2167 {
2168 gretl_matrix *vd, *V;
2169 double vjk;
2170 int h = horizon;
2171 int n = var->neqns;
2172 int imin, imax;
2173 int i, j, k, kk;
2174
2175 if (h <= 0) {
2176 h = default_VAR_horizon(dset);
2177 }
2178
2179 if (targ < 0) {
2180 /* do the whole thing */
2181 k = n * n;
2182 imin = 0;
2183 imax = n;
2184 } else {
2185 /* got a specific target */
2186 k = n;
2187 imin = targ;
2188 imax = targ + 1;
2189 }
2190
2191 V = gretl_matrix_alloc(h, k);
2192 if (V == NULL) {
2193 *err = E_ALLOC;
2194 return NULL;
2195 }
2196
2197 kk = 0;
2198 for (i=imin; i<imax && !*err; i++) {
2199 vd = gretl_VAR_get_fcast_decomp(var, i, h, err);
2200 if (!*err) {
2201 for (k=0; k<n; k++) {
2202 for (j=0; j<h; j++) {
2203 vjk = gretl_matrix_get(vd, j, k);
2204 gretl_matrix_set(V, j, kk, vjk / 100.0);
2205 }
2206 kk++;
2207 }
2208 gretl_matrix_free(vd);
2209 }
2210 }
2211
2212 /* If @shock is specific, not all, we now proceed to carve
2213 out the columns of @V that are actually wanted. This is
2214 not as efficient as it might be -- we're computing more
2215 than we need above -- but it's simple and will do for
2216 the present.
2217 */
2218
2219 if (!*err && shock >= 0) {
2220 gretl_matrix *V2 = NULL;
2221
2222 if (targ >= 0) {
2223 /* just one column wanted */
2224 V2 = gretl_column_vector_alloc(h);
2225 if (V2 != NULL) {
2226 for (j=0; j<h; j++) {
2227 V2->val[j] = gretl_matrix_get(V, j, shock);
2228 }
2229 }
2230 } else {
2231 /* n columns wanted */
2232 V2 = gretl_matrix_alloc(h, n);
2233 if (V2 != NULL) {
2234 k = shock;
2235 for (j=0; j<n; j++) {
2236 for (i=0; i<h; i++) {
2237 vjk = gretl_matrix_get(V, i, k);
2238 gretl_matrix_set(V2, i, j, vjk);
2239 }
2240 k += n;
2241 }
2242 }
2243 }
2244 if (V2 == NULL) {
2245 *err = E_ALLOC;
2246 } else {
2247 gretl_matrix_free(V);
2248 V = V2;
2249 }
2250 }
2251
2252 if (*err && V != NULL) {
2253 gretl_matrix_free(V);
2254 V = NULL;
2255 }
2256
2257 return V;
2258 }
2259
2260 gretl_matrix *
gretl_VAR_get_full_FEVD_matrix(const GRETL_VAR * var,const DATASET * dset,int * err)2261 gretl_VAR_get_full_FEVD_matrix (const GRETL_VAR *var, const DATASET *dset,
2262 int *err)
2263 {
2264 return gretl_VAR_get_FEVD_matrix(var, -1, -1, -1, dset, err);
2265 }
2266
var_max_order(const int * list,const DATASET * dset)2267 int var_max_order (const int *list, const DATASET *dset)
2268 {
2269 int T = dset->t2 - dset->t1 + 1;
2270 int nstoch = 0, ndet = 0;
2271 int gotsep = 0;
2272 int order = 1;
2273 int i;
2274
2275 for (i=1; i<=list[0]; i++) {
2276 if (list[i] == LISTSEP) {
2277 gotsep = 1;
2278 continue;
2279 }
2280 if (!gotsep) {
2281 nstoch++;
2282 } else {
2283 ndet++;
2284 }
2285 }
2286
2287 order = (T - ndet) / nstoch;
2288
2289 while (order > 0) {
2290 int t1 = (order > dset->t1)? order : dset->t1;
2291
2292 T = dset->t2 - t1 + 1;
2293 if (nstoch * order + ndet > T) {
2294 order--;
2295 } else {
2296 break;
2297 }
2298 }
2299
2300 return order - 1;
2301 }
2302
VAR_add_roots(GRETL_VAR * var)2303 static int VAR_add_roots (GRETL_VAR *var)
2304 {
2305 int err = 0;
2306
2307 if (var->A == NULL) {
2308 fprintf(stderr, "VAR_add_roots: var->A is missing\n");
2309 return E_DATA;
2310 }
2311
2312 /* save eigenvalues of companion form matrix */
2313 var->L = gretl_general_matrix_eigenvals(var->A, &err);
2314 #if 0
2315 gretl_matrix_print(var->A, "Companion form matrix");
2316 gretl_matrix_print(var->L, "Eigenvalues");
2317 #endif
2318
2319 if (err) {
2320 gretl_matrix_free(var->L);
2321 var->L = NULL;
2322 }
2323
2324 return err;
2325 }
2326
gretl_VAR_get_roots(GRETL_VAR * var,int * err)2327 const gretl_matrix *gretl_VAR_get_roots (GRETL_VAR *var, int *err)
2328 {
2329 if (var == NULL) {
2330 fprintf(stderr, "gretl_VAR_get_roots: VAR is NULL\n");
2331 *err = E_DATA;
2332 return NULL;
2333 }
2334
2335 if (var->L == NULL) {
2336 /* roots not computed yet */
2337 *err = VAR_add_roots(var);
2338 }
2339
2340 return var->L;
2341 }
2342
2343 /* used in context of LR tests: see vartest.c */
2344
gretl_VAR_ldet(GRETL_VAR * var,const gretl_matrix * E,int * err)2345 double gretl_VAR_ldet (GRETL_VAR *var, const gretl_matrix *E,
2346 int *err)
2347 {
2348 gretl_matrix *S = NULL;
2349 double ldet = NADBL;
2350
2351 S = gretl_matrix_alloc(var->neqns, var->neqns);
2352
2353 if (S == NULL) {
2354 *err = E_ALLOC;
2355 } else {
2356 gretl_matrix_multiply_mod(E, GRETL_MOD_TRANSPOSE,
2357 E, GRETL_MOD_NONE,
2358 S, GRETL_MOD_NONE);
2359 gretl_matrix_divide_by_scalar(S, var->T);
2360 ldet = gretl_vcv_log_determinant(S, err);
2361 gretl_matrix_free(S);
2362 }
2363
2364 return ldet;
2365 }
2366
2367 /* identify columns of the VAR X matrix that contain the final
2368 lag of an endogenous variable */
2369
omit_column(GRETL_VAR * var,int nl,int j)2370 static int omit_column (GRETL_VAR *var, int nl, int j)
2371 {
2372 if (var->ifc) {
2373 return j % nl == 0 && j < 1 + var->neqns * nl;
2374 } else {
2375 return (j+1) % nl == 0 && j < var->neqns * nl;
2376 }
2377 }
2378
2379 /* make and record residuals for LR test on last lag */
2380
VAR_short_residuals(GRETL_VAR * var,int * err)2381 static gretl_matrix *VAR_short_residuals (GRETL_VAR *var, int *err)
2382 {
2383 gretl_matrix *X = NULL;
2384 gretl_matrix *B = NULL;
2385 gretl_matrix *E = NULL;
2386 double x;
2387 /* note: removing one lag from each equation */
2388 int g = var->ncoeff - var->neqns;
2389 int j, t, k, nl;
2390
2391 E = gretl_matrix_alloc(var->T, var->neqns);
2392 X = gretl_matrix_alloc(var->T, g);
2393 B = gretl_matrix_alloc(g, var->neqns);
2394
2395 if (E == NULL || X == NULL || B == NULL) {
2396 *err = E_ALLOC;
2397 goto bailout;
2398 }
2399
2400 nl = var_n_lags(var);
2401
2402 k = 0;
2403 for (j=0; j<var->ncoeff; j++) {
2404 /* loop across the cols of var->X */
2405 if (j > 0 && omit_column(var, nl, j)) {
2406 continue;
2407 }
2408 for (t=0; t<var->T; t++) {
2409 x = gretl_matrix_get(var->X, t, j);
2410 gretl_matrix_set(X, t, k, x);
2411 }
2412 k++;
2413 }
2414
2415 /* put residuals from "short" estimation into E */
2416 *err = gretl_matrix_multi_ols(var->Y, X, B, E, NULL);
2417
2418 bailout:
2419
2420 gretl_matrix_free(X);
2421 gretl_matrix_free(B);
2422
2423 if (*err) {
2424 gretl_matrix_free(E);
2425 E = NULL;
2426 }
2427
2428 return E;
2429 }
2430
VAR_add_stats(GRETL_VAR * var,int code)2431 static int VAR_add_stats (GRETL_VAR *var, int code)
2432 {
2433 gretl_matrix *E1 = NULL;
2434 int err = 0;
2435
2436 if (var->order > 1 && code == VAR_ESTIMATE) {
2437 E1 = VAR_short_residuals(var, &err);
2438 }
2439
2440 var->S = gretl_matrix_alloc(var->neqns, var->neqns);
2441 if (var->S == NULL) {
2442 err = E_ALLOC;
2443 }
2444
2445 if (!err) {
2446 gretl_matrix_multiply_mod(var->E, GRETL_MOD_TRANSPOSE,
2447 var->E, GRETL_MOD_NONE,
2448 var->S, GRETL_MOD_NONE);
2449 /* for computing log-determinant, don't apply df
2450 correction (?) */
2451 gretl_matrix_divide_by_scalar(var->S, var->T);
2452 }
2453
2454 if (!err) {
2455 var->ldet = gretl_vcv_log_determinant(var->S, &err);
2456
2457 #if VAR_S_DFCORR
2458 /* Hmm, should we df-adjust var->S here? Note that this
2459 will affect the impulse response output; see also
2460 irfboot.c.
2461 */
2462 if (!err) {
2463 double cfac = var->T / (double) var->df;
2464
2465 gretl_matrix_multiply_by_scalar(var->S, cfac);
2466 }
2467 #endif
2468 }
2469
2470 if (!err) {
2471 int T = var->T;
2472 int g = var->neqns;
2473 int p = var->ncoeff;
2474 int k = g * p;
2475
2476 var->ll = -(g * T / 2.0) * (LN_2_PI + 1) - (T / 2.0) * var->ldet;
2477 var->AIC = (-2.0 * var->ll + 2.0 * k) / T;
2478 var->BIC = (-2.0 * var->ll + k * log(T)) / T;
2479 var->HQC = (-2.0 * var->ll + 2.0 * k * log(log(T))) / T;
2480 }
2481
2482 if (E1 != NULL) {
2483 if (!err) {
2484 VAR_LR_lag_test(var, E1);
2485 }
2486 gretl_matrix_free(E1);
2487 }
2488
2489 if (!err) {
2490 VAR_portmanteau_test(var);
2491 }
2492
2493 return err;
2494 }
2495
gretl_VAR_param_names(GRETL_VAR * v,char ** params,const DATASET * dset)2496 void gretl_VAR_param_names (GRETL_VAR *v, char **params,
2497 const DATASET *dset)
2498 {
2499 char lagstr[12];
2500 int i, j, n, k = 0;
2501
2502 if (v->detflags & DET_CONST) {
2503 strcpy(params[k++], dset->varname[0]);
2504 }
2505
2506 for (i=1; i<=v->ylist[0]; i++) {
2507 for (j=1; j<=v->order; j++) {
2508 if (!lag_wanted(v, j)) {
2509 continue;
2510 }
2511 sprintf(lagstr, "_%d", j);
2512 n = strlen(lagstr);
2513 if (v->ci == VECM) {
2514 strcpy(params[k], "d_");
2515 n += 2;
2516 }
2517 strncat(params[k], dset->varname[v->ylist[i]],
2518 VNAMELEN - n - 1);
2519 strncat(params[k], lagstr, n);
2520 k++;
2521 }
2522 }
2523
2524 if (v->xlist != NULL) {
2525 for (i=1; i<=v->xlist[0]; i++) {
2526 strcpy(params[k++], dset->varname[v->xlist[i]]);
2527 }
2528 }
2529
2530 if (v->detflags & DET_SEAS) {
2531 for (i=1; i<dset->pd; i++) {
2532 sprintf(params[k++], "S%d", i);
2533 }
2534 }
2535
2536 if (v->detflags & DET_TREND) {
2537 strcpy(params[k++], "time");
2538 }
2539
2540 if (v->ci == VECM) {
2541 int rank = jrank(v);
2542
2543 for (i=0; i<rank; i++) {
2544 sprintf(params[k++], "EC%d", i+1);
2545 }
2546 }
2547 }
2548
2549 /* public because called from irfboot.c */
2550
VAR_write_A_matrix(GRETL_VAR * v,GretlMatrixMod mod)2551 void VAR_write_A_matrix (GRETL_VAR *v, GretlMatrixMod mod)
2552 {
2553 int i, ii, j, k, lag;
2554 int n = v->neqns;
2555 int dim = n * v->order;
2556 double bij;
2557
2558 for (j=0; j<n; j++) {
2559 k = lag = ii = 0;
2560 for (i=0; i<dim; i++) {
2561 if (lag_wanted(v, lag+1)) {
2562 bij = gretl_matrix_get(v->B, ii + v->ifc, j);
2563 ii++;
2564 } else {
2565 bij = 0;
2566 }
2567 if (mod == GRETL_MOD_TRANSPOSE) {
2568 gretl_matrix_set(v->A, n * lag + k, j, bij);
2569 } else {
2570 gretl_matrix_set(v->A, j, n * lag + k, bij);
2571 }
2572 if (lag < v->order - 1) {
2573 lag++;
2574 } else {
2575 lag = 0;
2576 k++;
2577 }
2578 }
2579 }
2580 }
2581
VAR_depvar_name(GRETL_VAR * var,int i,const char * yname)2582 static int VAR_depvar_name (GRETL_VAR *var, int i, const char *yname)
2583 {
2584 MODEL *pmod = var->models[i];
2585
2586 if (var->ci == VAR) {
2587 pmod->depvar = gretl_strdup(yname);
2588 } else {
2589 pmod->depvar = malloc(VNAMELEN);
2590 if (pmod->depvar != NULL) {
2591 strcpy(pmod->depvar, "d_");
2592 strncat(pmod->depvar, yname, VNAMELEN - 3);
2593 }
2594 }
2595
2596 return (pmod->depvar == NULL)? E_ALLOC: 0;
2597 }
2598
2599 /* transcribe the per-equation output from a VAR or VECM into
2600 MODEL structs, if wanted */
2601
transcribe_VAR_models(GRETL_VAR * var,const DATASET * dset,const gretl_matrix * XTX)2602 int transcribe_VAR_models (GRETL_VAR *var,
2603 const DATASET *dset,
2604 const gretl_matrix *XTX)
2605 {
2606 MODEL *pmod;
2607 char **params = NULL;
2608 int yno, N = dset->n;
2609 int ecm = (var->ci == VECM);
2610 double x;
2611 int i, j, jmax;
2612 int err = 0;
2613
2614 params = strings_array_new_with_length(var->ncoeff, VNAMELEN);
2615 if (params == NULL) {
2616 return E_ALLOC;
2617 }
2618
2619 gretl_VAR_param_names(var, params, dset);
2620
2621 jmax = (var->B != NULL)? var->B->rows : 0;
2622
2623 for (i=0; i<var->neqns && !err; i++) {
2624 yno = var->ylist[i+1];
2625
2626 pmod = var->models[i];
2627 pmod->ID = i + 1;
2628 pmod->ci = (ecm)? OLS : VAR;
2629 pmod->aux = (ecm)? AUX_VECM: AUX_VAR;
2630
2631 pmod->full_n = N;
2632 pmod->nobs = var->T;
2633 pmod->t1 = var->t1;
2634 pmod->t2 = var->t2;
2635 pmod->ncoeff = var->ncoeff;
2636 pmod->ifc = var->ifc;
2637 pmod->dfn = var->ncoeff - pmod->ifc;
2638 pmod->dfd = (ecm)? var->df : pmod->nobs - pmod->ncoeff;
2639
2640 err = gretl_model_allocate_storage(pmod);
2641
2642 if (!err) {
2643 VAR_depvar_name(var, i, dset->varname[yno]);
2644 if (i == 0) {
2645 pmod->params = params;
2646 } else {
2647 pmod->params = strings_array_dup(params, var->ncoeff);
2648 if (pmod->params == NULL) {
2649 err = E_ALLOC;
2650 }
2651 }
2652 }
2653
2654 if (!err) {
2655 pmod->nparams = var->ncoeff;
2656 pmod->list = gretl_list_new(1);
2657 if (pmod->list == NULL) {
2658 err = E_ALLOC;
2659 }
2660 }
2661
2662 if (!err) {
2663 pmod->list[1] = yno;
2664 set_VAR_model_stats(var, i);
2665 for (j=0; j<jmax; j++) {
2666 pmod->coeff[j] = gretl_matrix_get(var->B, j, i);
2667 if (XTX != NULL) {
2668 if (XTX->rows <= var->ncoeff) {
2669 x = gretl_matrix_get(XTX, j, j);
2670 pmod->sderr[j] = pmod->sigma * sqrt(x);
2671 } else {
2672 int jj = i * var->ncoeff + j;
2673
2674 x = gretl_matrix_get(XTX, jj, jj);
2675 pmod->sderr[j] = pmod->sigma * sqrt(x);
2676 }
2677 }
2678 }
2679 }
2680 }
2681
2682 return err;
2683 }
2684
lags_from_laglist(const int * llist,int * err)2685 static int *lags_from_laglist (const int *llist, int *err)
2686 {
2687 int *lags = NULL;
2688
2689 if (llist[0] == 0) {
2690 *err = E_DATA;
2691 } else {
2692 lags = gretl_list_copy(llist);
2693 if (lags == NULL) {
2694 *err = E_ALLOC;
2695 }
2696 }
2697
2698 if (lags != NULL) {
2699 gretl_list_sort(lags);
2700 if (lags[1] < 1) {
2701 *err = E_DATA;
2702 free(lags);
2703 lags = NULL;
2704 }
2705 }
2706
2707 return lags;
2708 }
2709
2710 /**
2711 * gretl_VAR:
2712 * @order: lag order for the VAR.
2713 * @laglist: specific list of lags, or NULL.
2714 * @list: specification for the first model in the set.
2715 * @dset: dataset struct.
2716 * @opt: if includes %OPT_R, use robust VCV;
2717 * if includes %OPT_H, use HAC VCV;
2718 * if includes %OPT_I, print impulse responses;
2719 * if includes %OPT_F, print forecast variance decompositions;
2720 * if includes %OPT_D, add seasonal dummies;
2721 * if includes %OPT_M, set min lag for lag-selection.
2722 * if includes %OPT_N, do not include a constant.
2723 * if includes %OPT_Q, do not show individual regressions.
2724 * if includes %OPT_T, include a linear trend.
2725 * if includes %OPT_L, test for optimal lag length (only).
2726 * if includes %OPT_S, silent (no printing).
2727 * @prn: gretl printing struct.
2728 * @err: location to receive error code.
2729 *
2730 * Estimate a vector auto-regression (VAR), print and save
2731 * the results.
2732 *
2733 * Returns: pointer to VAR struct, which may be %NULL on error.
2734 */
2735
gretl_VAR(int order,int * laglist,int * list,const DATASET * dset,gretlopt opt,PRN * prn,int * err)2736 GRETL_VAR *gretl_VAR (int order, int *laglist, int *list,
2737 const DATASET *dset, gretlopt opt,
2738 PRN *prn, int *err)
2739 {
2740 GRETL_VAR *var = NULL;
2741 int code = (opt & OPT_L)? VAR_LAGSEL : VAR_ESTIMATE;
2742 int *lags = NULL;
2743
2744 if ((opt & OPT_M) && !(opt & OPT_L)) {
2745 gretl_errmsg_sprintf("%s: inapplicable option", "--minlag");
2746 *err = E_BADOPT;
2747 return NULL;
2748 }
2749
2750 if (laglist != NULL) {
2751 lags = lags_from_laglist(laglist, err);
2752 if (*err) {
2753 return NULL;
2754 }
2755 }
2756
2757 /* allocation and initial set-up */
2758 var = gretl_VAR_new(code, order, 0, lags, list, dset,
2759 opt, err);
2760 if (var == NULL) {
2761 return NULL;
2762 }
2763
2764 /* run the regressions */
2765 if (getenv("VAR_USE_QR") != NULL) {
2766 *err = gretl_matrix_QR_ols(var->Y, var->X,
2767 var->B, var->E,
2768 &var->XTX, NULL);
2769 } else {
2770 /* use Cholesky or QR as needed */
2771 *err = gretl_matrix_multi_ols(var->Y, var->X,
2772 var->B, var->E,
2773 &var->XTX);
2774 }
2775
2776 if (!*err) {
2777 if (code == VAR_LAGSEL) {
2778 /* doing lag-length selection */
2779 *err = VAR_add_stats(var, code);
2780 if (!*err) {
2781 *err = VAR_do_lagsel(var, dset, opt, prn);
2782 }
2783 } else {
2784 /* regular VAR estimation */
2785 *err = transcribe_VAR_models(var, dset, NULL);
2786
2787 if (!*err) {
2788 VAR_write_A_matrix(var, GRETL_MOD_NONE);
2789 /* note: the following also sets std. errors */
2790 *err = VAR_wald_omit_tests(var);
2791 }
2792
2793 if (!*err) {
2794 *err = VAR_add_stats(var, code);
2795 }
2796
2797 if (!*err) {
2798 *err = gretl_VAR_do_error_decomp(var->S, var->C, NULL);
2799 }
2800
2801 if (!*err && prn != NULL) {
2802 gretl_VAR_print(var, dset, opt, prn);
2803 }
2804 }
2805 }
2806
2807 if (code == VAR_LAGSEL || (*err && var != NULL)) {
2808 gretl_VAR_free(var);
2809 var = NULL;
2810 }
2811
2812 return var;
2813 }
2814
2815 static void
print_johansen_sigmas(const JohansenInfo * jv,PRN * prn)2816 print_johansen_sigmas (const JohansenInfo *jv, PRN *prn)
2817 {
2818 int i, j;
2819
2820 pprintf(prn, "%s\n\n", _("Sample variance-covariance matrices for residuals"));
2821
2822 pprintf(prn, " %s (S00)\n\n", _("VAR system in first differences"));
2823 for (i=0; i<jv->S00->rows; i++) {
2824 for (j=0; j<jv->S00->rows; j++) {
2825 pprintf(prn, "%#12.5g", gretl_matrix_get(jv->S00, i, j));
2826 }
2827 pputc(prn, '\n');
2828 }
2829
2830 pprintf(prn, "\n %s (S11)\n\n", _("System with levels as dependent variable"));
2831 for (i=0; i<jv->S11->rows; i++) {
2832 for (j=0; j<jv->S11->rows; j++) {
2833 pprintf(prn, "%#12.5g", gretl_matrix_get(jv->S11, i, j));
2834 }
2835 pputc(prn, '\n');
2836 }
2837
2838 pprintf(prn, "\n %s (S01)\n\n", _("Cross-products"));
2839 for (i=0; i<jv->S01->rows; i++) {
2840 for (j=0; j<jv->S01->cols; j++) {
2841 pprintf(prn, "%#12.5g", gretl_matrix_get(jv->S01, i, j));
2842 }
2843 pputc(prn, '\n');
2844 }
2845
2846 pputc(prn, '\n');
2847 }
2848
2849 /* called from gretl_restriction_finalize() if the target
2850 is a VECM */
2851
gretl_VECM_test(GRETL_VAR * vecm,gretl_restriction * rset,const DATASET * dset,gretlopt opt,PRN * prn)2852 int gretl_VECM_test (GRETL_VAR *vecm,
2853 gretl_restriction *rset,
2854 const DATASET *dset,
2855 gretlopt opt,
2856 PRN *prn)
2857 {
2858 int (*jfun) (GRETL_VAR *, gretl_restriction *,
2859 const DATASET *, gretlopt opt, PRN *);
2860 int err = 0;
2861
2862 if (vecm->jinfo == NULL || rset == NULL) {
2863 return E_DATA;
2864 }
2865
2866 gretl_error_clear();
2867
2868 jfun = get_plugin_function("vecm_test_restriction");
2869
2870 if (jfun == NULL) {
2871 err = 1;
2872 } else {
2873 err = (*jfun) (vecm, rset, dset, opt, prn);
2874 }
2875
2876 return err;
2877 }
2878
var_add_full_vcv(GRETL_VAR * var)2879 static int var_add_full_vcv (GRETL_VAR *var)
2880 {
2881 gretl_matrix *V;
2882 MODEL *pmod;
2883 double vij;
2884 int i, k, nc;
2885 int bi, bj;
2886 int vi, vj;
2887 int err = 0;
2888
2889 nc = var->neqns * var->ncoeff;
2890 V = gretl_zero_matrix_new(nc, nc);
2891 if (V == NULL) {
2892 return E_ALLOC;
2893 }
2894
2895 k = var->ncoeff;
2896
2897 vi = vj = 0;
2898 for (i=0; i<var->neqns; i++) {
2899 pmod = var->models[i];
2900 if (pmod->vcv == NULL) {
2901 err = makevcv(pmod, pmod->sigma);
2902 if (err) {
2903 break;
2904 }
2905 }
2906 for (bi=0; bi<k; bi++) {
2907 for (bj=0; bj<k; bj++) {
2908 vij = gretl_model_get_vcv_element(pmod, bi, bj, k);
2909 gretl_matrix_set(V, vi+bi, vj+bj, vij);
2910 }
2911 }
2912 vi += k;
2913 vj += k;
2914 }
2915
2916 if (err) {
2917 gretl_matrix_free(V);
2918 } else {
2919 var->V = V;
2920 }
2921
2922 return err;
2923 }
2924
2925 /* called from gretl_restriction_finalize() if the target
2926 is a plain VAR */
2927
gretl_VAR_test(GRETL_VAR * var,gretl_restriction * rset,const DATASET * dset,gretlopt opt,PRN * prn)2928 int gretl_VAR_test (GRETL_VAR *var,
2929 gretl_restriction *rset,
2930 const DATASET *dset,
2931 gretlopt opt,
2932 PRN *prn)
2933 {
2934 int err = 0;
2935
2936 if (rset == NULL) {
2937 return E_DATA;
2938 }
2939
2940 gretl_error_clear();
2941
2942 if (var->V == NULL) {
2943 err = var_add_full_vcv(var);
2944 }
2945
2946 if (!err) {
2947 const gretl_matrix *R = rset_get_R_matrix(rset);
2948 const gretl_matrix *q = rset_get_q_matrix(rset);
2949 int dfu = var->T - var->ncoeff;
2950 int r = var->B->rows;
2951 int c = var->B->cols;
2952
2953 /* pass var->B as a column vector */
2954 gretl_matrix_reuse(var->B, r*c, 1);
2955 err = multi_eqn_wald_test(var->B, var->V, R, q, dfu,
2956 opt, prn);
2957 gretl_matrix_reuse(var->B, r, c);
2958 }
2959
2960 return err;
2961 }
2962
2963 static int
johansen_test_complete(GRETL_VAR * jvar,const DATASET * dset,gretlopt opt,PRN * prn)2964 johansen_test_complete (GRETL_VAR *jvar, const DATASET *dset,
2965 gretlopt opt, PRN *prn)
2966 {
2967 int (*jfun) (GRETL_VAR *, const DATASET *, gretlopt, PRN *);
2968 int err = 0;
2969
2970 gretl_error_clear();
2971
2972 jfun = get_plugin_function("johansen_coint_test");
2973
2974 if (jfun == NULL) {
2975 err = 1;
2976 } else {
2977 err = (* jfun) (jvar, dset, opt, prn);
2978 }
2979
2980 return err;
2981 }
2982
2983 static int
johansen_estimate_complete(GRETL_VAR * jvar,gretl_restriction * rset,const DATASET * dset,PRN * prn)2984 johansen_estimate_complete (GRETL_VAR *jvar, gretl_restriction *rset,
2985 const DATASET *dset, PRN *prn)
2986 {
2987 int (*jfun) (GRETL_VAR *, gretl_restriction *,
2988 const DATASET *, PRN *);
2989 int err = 0;
2990
2991 gretl_error_clear();
2992
2993 jfun = get_plugin_function("johansen_estimate");
2994
2995 if (jfun == NULL) {
2996 err = 1;
2997 } else {
2998 err = (* jfun) (jvar, rset, dset, prn);
2999 }
3000
3001 return err;
3002 }
3003
3004 /* N.B. we allow for the possibility that this allocation has
3005 already been done */
3006
allocate_johansen_extra_matrices(GRETL_VAR * v)3007 static int allocate_johansen_extra_matrices (GRETL_VAR *v)
3008 {
3009 if (v->jinfo->R0 != NULL &&
3010 v->jinfo->S00 != NULL &&
3011 v->jinfo->YY != NULL) {
3012 return 0;
3013 } else {
3014 int xc = v->X == NULL ? 0 : v->X->cols;
3015 int p0 = v->neqns;
3016 int p1 = p0 + n_restricted_terms(v);
3017 int p = p0 + p1;
3018
3019 clear_gretl_matrix_err();
3020
3021 if (v->jinfo->R0 == NULL) {
3022 v->jinfo->R0 = gretl_matrix_alloc(v->T, p0);
3023 v->jinfo->R1 = gretl_matrix_alloc(v->T, p1);
3024 }
3025
3026 if (v->jinfo->S00 == NULL) {
3027 v->jinfo->S00 = gretl_matrix_alloc(p0, p0);
3028 v->jinfo->S11 = gretl_matrix_alloc(p1, p1);
3029 v->jinfo->S01 = gretl_matrix_alloc(p0, p1);
3030 }
3031
3032 if (xc > 0 && v->ncoeff > 0 && v->jinfo->YY == NULL) {
3033 v->jinfo->YY = gretl_matrix_alloc(v->T, p);
3034 v->jinfo->RR = gretl_matrix_alloc(v->T, p);
3035 v->jinfo->BB = gretl_matrix_alloc(v->X->cols, p);
3036 }
3037
3038 return get_gretl_matrix_err();
3039 }
3040 }
3041
johansen_fill_S_matrices(GRETL_VAR * v)3042 static void johansen_fill_S_matrices (GRETL_VAR *v)
3043 {
3044 gretl_matrix_multiply_mod(v->jinfo->R0, GRETL_MOD_TRANSPOSE,
3045 v->jinfo->R0, GRETL_MOD_NONE,
3046 v->jinfo->S00, GRETL_MOD_NONE);
3047 gretl_matrix_multiply_mod(v->jinfo->R1, GRETL_MOD_TRANSPOSE,
3048 v->jinfo->R1, GRETL_MOD_NONE,
3049 v->jinfo->S11, GRETL_MOD_NONE);
3050 gretl_matrix_multiply_mod(v->jinfo->R0, GRETL_MOD_TRANSPOSE,
3051 v->jinfo->R1, GRETL_MOD_NONE,
3052 v->jinfo->S01, GRETL_MOD_NONE);
3053
3054 gretl_matrix_divide_by_scalar(v->jinfo->S00, v->T);
3055 gretl_matrix_divide_by_scalar(v->jinfo->S11, v->T);
3056 gretl_matrix_divide_by_scalar(v->jinfo->S01, v->T);
3057 }
3058
3059 /* We come here if there are no regressors at stage 1 of
3060 the Johansen procedure. This happens if the associated
3061 VAR is of order 1 (becomes order 0 on differencing)
3062 and there are no unrestricted exogenous variables.
3063 The residuals are then just the values of the
3064 variables themselves.
3065 */
3066
johansen_degenerate_stage_1(GRETL_VAR * v,const DATASET * dset)3067 static int johansen_degenerate_stage_1 (GRETL_VAR *v,
3068 const DATASET *dset)
3069 {
3070 const double **Z = (const double **) dset->Z;
3071 gretl_matrix *R0 = v->jinfo->R0;
3072 gretl_matrix *R1 = v->jinfo->R1;
3073 int i, vi, s, t, j = 0;
3074
3075 #if 0
3076 fprintf(stderr, "degenerate stage 1\n");
3077 #endif
3078
3079 for (i=0; i<v->neqns; i++) {
3080 vi = v->ylist[i+1];
3081 s = 0;
3082 for (t=v->t1; t<=v->t2; t++) {
3083 gretl_matrix_set(R0, s, j, Z[vi][t] - Z[vi][t-1]);
3084 gretl_matrix_set(R1, s, j, Z[vi][t-1]);
3085 s++;
3086 }
3087 j++;
3088 }
3089
3090 if (auto_restr(v)) {
3091 int trend = (jcode(v) == J_REST_TREND);
3092
3093 for (t=0; t<v->T; t++) {
3094 gretl_matrix_set(R1, t, j, (trend)? v->t1 + t : 1);
3095 }
3096 j++;
3097 }
3098
3099 if (v->rlist != NULL) {
3100 for (i=0; i<v->rlist[0]; i++) {
3101 vi = v->rlist[i+1];
3102 s = 0;
3103 for (t=v->t1; t<=v->t2; t++) {
3104 gretl_matrix_set(R1, s++, j, Z[vi][t]); /* was t-1 */
3105 }
3106 j++;
3107 }
3108 }
3109
3110 johansen_fill_S_matrices(v);
3111
3112 return 0;
3113 }
3114
print_stage_1_coeffs(GRETL_VAR * v,gretl_matrix * B,PRN * prn)3115 static void print_stage_1_coeffs (GRETL_VAR *v, gretl_matrix *B,
3116 PRN *prn)
3117 {
3118 gretl_matrix tmp;
3119
3120 gretl_matrix_init(&tmp);
3121
3122 tmp.rows = B->rows;
3123 tmp.cols = v->neqns;
3124 tmp.val = B->val;
3125
3126 gretl_matrix_print_to_prn(&tmp, "\nCoefficients, VAR in differences",
3127 prn);
3128
3129 tmp.cols += n_restricted_terms(v);
3130 tmp.val += v->neqns * tmp.rows;
3131
3132 gretl_matrix_print_to_prn(&tmp, "Coefficients, eqns in lagged levels",
3133 prn);
3134 }
3135
3136 #define JVAR_USE_SVD 1
3137
johansen_partition_residuals(GRETL_VAR * v)3138 static void johansen_partition_residuals (GRETL_VAR *v)
3139 {
3140 int n = v->neqns * v->T;
3141 int m = n + n_restricted_terms(v) * v->T;
3142 gretl_matrix *R = v->jinfo->RR;
3143
3144 memcpy(v->jinfo->R0->val, R->val, n * sizeof(double));
3145 memcpy(v->jinfo->R1->val, R->val + n, m * sizeof(double));
3146 }
3147
3148 /* For Johansen analysis: estimate VAR in differences along with the
3149 other auxiliary regressions required to compute the relevant
3150 matrices of residuals, for concentration of the log-likelihood.
3151 Then compute S00, S11, S01. See for example James Hamilton,
3152 "Time Series Analysis", section 20.2.
3153 */
3154
johansen_stage_1(GRETL_VAR * v,const DATASET * dset,gretlopt opt,PRN * prn)3155 int johansen_stage_1 (GRETL_VAR *v, const DATASET *dset,
3156 gretlopt opt, PRN *prn)
3157 {
3158 int err;
3159
3160 err = allocate_johansen_extra_matrices(v);
3161 if (err) {
3162 fprintf(stderr, "allocate_johansen_extra_matrices: err = %d\n", err);
3163 return err;
3164 }
3165
3166 #if 0
3167 fprintf(stderr, "johansen_stage_1: ncoeff = %d\n", v->ncoeff);
3168 #endif
3169
3170 if (v->ncoeff == 0 || v->jinfo->BB == NULL) {
3171 /* nothing to concentrate out */
3172 if (opt & OPT_V) {
3173 pputs(prn, "\nNo initial VAR estimation is required\n\n");
3174 }
3175 johansen_degenerate_stage_1(v, dset);
3176 } else {
3177 gretl_matrix *Y = v->jinfo->YY;
3178 gretl_matrix *B = v->jinfo->BB;
3179 gretl_matrix *R = v->jinfo->RR;
3180
3181 VECM_fill_Y(v, dset, Y);
3182
3183 #if 0
3184 gretl_matrix_print(v->X, "X in johansen_stage_1");
3185 gretl_matrix_print(Y, "Y in johansen_stage_1");
3186 #endif
3187
3188 #if JVAR_USE_SVD
3189 err = gretl_matrix_multi_SVD_ols(Y, v->X, B, R, NULL);
3190 #else
3191 err = gretl_matrix_multi_ols(Y, v->X, B, R, NULL);
3192 #endif
3193
3194 if (!err) {
3195 if (opt & OPT_V) {
3196 print_stage_1_coeffs(v, B, prn);
3197 }
3198 johansen_partition_residuals(v);
3199 johansen_fill_S_matrices(v);
3200 }
3201 }
3202
3203 return err;
3204 }
3205
opt_from_jcode(JohansenCode jc)3206 static gretlopt opt_from_jcode (JohansenCode jc)
3207 {
3208 gretlopt opt = OPT_NONE;
3209
3210 if (jc == J_NO_CONST) {
3211 opt = OPT_N;
3212 } else if (jc == J_UNREST_TREND) {
3213 opt = OPT_T;
3214 } else if (jc == J_REST_CONST) {
3215 opt = OPT_R;
3216 } else if (jc == J_REST_TREND) {
3217 opt = OPT_A;
3218 }
3219
3220 return opt;
3221 }
3222
jcode_from_opt(gretlopt opt)3223 static JohansenCode jcode_from_opt (gretlopt opt)
3224 {
3225 JohansenCode jc = J_UNREST_CONST;
3226
3227 if (opt & OPT_N) {
3228 jc = J_NO_CONST;
3229 } else if (opt & OPT_T) {
3230 jc = J_UNREST_TREND;
3231 } else if (opt & OPT_R) {
3232 jc = J_REST_CONST;
3233 } else if (opt & OPT_A) {
3234 jc = J_REST_TREND;
3235 }
3236
3237 return jc;
3238 }
3239
3240 static JohansenInfo *
johansen_info_new(GRETL_VAR * var,int rank,gretlopt opt)3241 johansen_info_new (GRETL_VAR *var, int rank, gretlopt opt)
3242 {
3243 JohansenInfo *jv = malloc(sizeof *jv);
3244
3245 if (jv == NULL) {
3246 return NULL;
3247 }
3248
3249 jv->ID = 0;
3250 jv->code = jcode_from_opt(opt);
3251 jv->rank = rank;
3252 jv->seasonals = 0;
3253
3254 jv->R0 = NULL;
3255 jv->R1 = NULL;
3256 jv->S00 = NULL;
3257 jv->S11 = NULL;
3258 jv->S01 = NULL;
3259 jv->evals = NULL;
3260 jv->Beta = NULL;
3261 jv->Alpha = NULL;
3262 jv->Bse = NULL;
3263 jv->Ase = NULL;
3264 jv->Bvar = NULL;
3265 jv->R = NULL;
3266 jv->q = NULL;
3267 jv->Ra = NULL;
3268 jv->qa = NULL;
3269 jv->YY = NULL;
3270 jv->RR = NULL;
3271 jv->BB = NULL;
3272
3273 jv->ll0 = jv->prior_ll = NADBL;
3274 jv->lrdf = jv->prior_df = 0;
3275
3276 if (rank > 0) {
3277 jv->Alpha = gretl_zero_matrix_new(var->neqns, rank);
3278 if (jv->Alpha == NULL) {
3279 free(jv);
3280 jv = NULL;
3281 }
3282 }
3283
3284 return jv;
3285 }
3286
coint_test_header(const GRETL_VAR * v,const DATASET * dset,PRN * prn)3287 static void coint_test_header (const GRETL_VAR *v,
3288 const DATASET *dset,
3289 PRN *prn)
3290 {
3291 char stobs[OBSLEN], endobs[OBSLEN];
3292
3293 pprintf(prn, "\n%s:\n", _("Johansen test"));
3294 pprintf(prn, "%s = %d\n", _("Number of equations"), v->neqns);
3295 pprintf(prn, "%s = %d\n", _("Lag order"), v->order + 1);
3296 pprintf(prn, "%s: %s - %s (T = %d)\n", _("Estimation period"),
3297 ntolabel(stobs, v->t1, dset),
3298 ntolabel(endobs, v->t2, dset), v->T);
3299
3300 }
3301
jvar_check_allocation(GRETL_VAR * v,const DATASET * dset)3302 static int jvar_check_allocation (GRETL_VAR *v, const DATASET *dset)
3303 {
3304 int err = VAR_add_models(v, dset);
3305
3306 if (!err) {
3307 err = VAR_allocate_companion_matrix(v);
3308 }
3309 if (!err) {
3310 err = VAR_allocate_cholesky_matrix(v);
3311 }
3312 if (!err) {
3313 err = VAR_allocate_residuals_matrix(v);
3314 }
3315
3316 return err;
3317 }
3318
3319 /* Driver function for Johansen analysis. An appropriately
3320 initialized "jvar" must have been set up already.
3321 */
3322
3323 static int
johansen_driver(GRETL_VAR * jvar,gretl_restriction * rset,const DATASET * dset,gretlopt opt,PRN * prn)3324 johansen_driver (GRETL_VAR *jvar,
3325 gretl_restriction *rset,
3326 const DATASET *dset,
3327 gretlopt opt, PRN *prn)
3328 {
3329 int r = jrank(jvar);
3330
3331 if (r == 0) {
3332 /* doing cointegration test */
3333 coint_test_header(jvar, dset, prn);
3334 }
3335
3336 jvar->err = johansen_stage_1(jvar, dset, opt, prn);
3337 if (jvar->err) {
3338 return jvar->err;
3339 }
3340
3341 if (opt & OPT_V) {
3342 print_johansen_sigmas(jvar->jinfo, prn);
3343 }
3344
3345 if (r > 0) {
3346 /* estimating VECM, not just doing cointegration test */
3347 jvar->err = jvar_check_allocation(jvar, dset);
3348 }
3349
3350 if (!jvar->err) {
3351 /* call the johansen plugin to finish the job */
3352 if (r > 0) {
3353 jvar->err = johansen_estimate_complete(jvar, rset, dset,
3354 prn);
3355 } else {
3356 jvar->err = johansen_test_complete(jvar, dset, opt, prn);
3357 }
3358 }
3359
3360 return jvar->err;
3361 }
3362
3363 static GRETL_VAR *
johansen_wrapper(int code,int order,int rank,const int * lags,const int * list,gretl_restriction * rset,const DATASET * dset,gretlopt opt,PRN * prn,int * err)3364 johansen_wrapper (int code, int order, int rank,
3365 const int *lags, const int *list,
3366 gretl_restriction *rset,
3367 const DATASET *dset,
3368 gretlopt opt, PRN *prn, int *err)
3369 {
3370 GRETL_VAR *jvar;
3371
3372 jvar = gretl_VAR_new(code, order - 1, rank, lags, list, dset,
3373 opt, err);
3374 if (jvar != NULL && !jvar->err) {
3375 *err = jvar->err = johansen_driver(jvar, rset, dset, opt, prn);
3376 }
3377
3378 return jvar;
3379 }
3380
3381 /**
3382 * johansen_test:
3383 * @order: lag order for test.
3384 * @list: list of variables to test for cointegration.
3385 * @dset: dataset struct.
3386 * @opt: %OPT_A: include constant plus restricted trend; %OPT_D:
3387 * include centered seasonals; %OPT_N: no constant; %OPT_R:
3388 * restricted constant; %OPT_T: constant and unrestricted trend
3389 * (note: default "case" is unrestricted constant).
3390 * %OPT_V: produce verbose results.
3391 * @prn: gretl printing struct.
3392 *
3393 * Carries out the Johansen test for cointegration.
3394 *
3395 * Returns: pointer to struct containing information on
3396 * the test.
3397 */
3398
johansen_test(int order,const int * list,const DATASET * dset,gretlopt opt,PRN * prn)3399 GRETL_VAR *johansen_test (int order, const int *list,
3400 const DATASET *dset,
3401 gretlopt opt, PRN *prn)
3402 {
3403 int err = 0;
3404
3405 return johansen_wrapper(VECM_CTEST, order, 0, NULL, list, NULL,
3406 dset, opt, prn, &err);
3407 }
3408
3409 /**
3410 * johansen_test_simple:
3411 * @order: lag order for test.
3412 * @list: list of variables to test for cointegration.
3413 * @dset: dataset struct.
3414 * @opt: %OPT_A: include constant plus restricted trend; %OPT_D:
3415 * include centered seasonals; %OPT_N: no constant; %OPT_R:
3416 * restricted constant; %OPT_T: constant and unrestricted trend
3417 * (note: default "case" is unrestricted constant);
3418 * %OPT_V: produce verbose results; %OPT_Q: just print the tests;
3419 * %OPT_S: don't print anything.
3420 * @prn: gretl printing struct.
3421 *
3422 * Carries out the Johansen test for cointegration and prints the
3423 * results (but unlike johansen_test(), does not return the
3424 * allocated results in VAR form).
3425 *
3426 * Returns: 0 on success, non-zero code on error.
3427 */
3428
johansen_test_simple(int order,const int * list,const DATASET * dset,gretlopt opt,PRN * prn)3429 int johansen_test_simple (int order, const int *list,
3430 const DATASET *dset,
3431 gretlopt opt, PRN *prn)
3432 {
3433 GRETL_VAR *jvar = NULL;
3434 int err = 0;
3435
3436 jvar = johansen_wrapper(VECM_CTEST, order, 0, NULL, list, NULL,
3437 dset, opt, (opt & OPT_S)? NULL : prn,
3438 &err);
3439
3440 if (jvar != NULL) {
3441 gretl_VAR_free(jvar);
3442 }
3443
3444 return err;
3445 }
3446
3447 /**
3448 * gretl_VECM:
3449 * @order: lag order.
3450 * @rank: cointegration rank.
3451 * @list: list of endogenous variables, possibly plus
3452 * exogenous variables.
3453 * @dset: dataset struct.
3454 * @opt: may include OPT_N ("nc"), OPT_R ("rc"), OPT_A ("crt"),
3455 * OPT_T ("ct"), OPT_D (include seasonals), OPT_F (show variance
3456 * decompositions), OPT_I (show impulse responses), OPT_V
3457 * (verbose operation), OPT_Q (quiet), OPT_S (silent).
3458 * @prn: gretl printing struct.
3459 * @err: location to receive error code.
3460 *
3461 * Returns: pointer to struct containing a VECM system,
3462 * or %NULL on failure.
3463 */
3464
gretl_VECM(int order,int rank,int * list,const DATASET * dset,gretlopt opt,PRN * prn,int * err)3465 GRETL_VAR *gretl_VECM (int order, int rank, int *list,
3466 const DATASET *dset,
3467 gretlopt opt, PRN *prn, int *err)
3468 {
3469 GRETL_VAR *jvar = NULL;
3470 int *lags = NULL;
3471
3472 if (rank <= 0) {
3473 gretl_errmsg_sprintf(_("vecm: rank %d is out of bounds"), rank);
3474 *err = E_DATA;
3475 return NULL;
3476 }
3477
3478 jvar = johansen_wrapper(VECM_ESTIMATE, order, rank, lags, list,
3479 NULL, dset, opt, prn, err);
3480
3481 if (jvar != NULL && !jvar->err) {
3482 gretl_VAR_print(jvar, dset, opt, prn);
3483 }
3484
3485 return jvar;
3486 }
3487
VAR_list_composite(const int * ylist,const int * xlist,const int * rlist)3488 int *VAR_list_composite (const int *ylist, const int *xlist,
3489 const int *rlist)
3490 {
3491 int *big = NULL;
3492 int i, k, n = ylist[0];
3493
3494 if (xlist != NULL && xlist[0] > 0) {
3495 n += xlist[0] + 1;
3496 }
3497
3498 if (rlist != NULL && rlist[0] > 0) {
3499 n += rlist[0] + 1;
3500 if (xlist == NULL || xlist[0] == 0) {
3501 /* extra separator needed */
3502 n++;
3503 }
3504 }
3505
3506 big = gretl_list_new(n);
3507 if (big == NULL) {
3508 return NULL;
3509 }
3510
3511 k = 1;
3512
3513 for (i=1; i<=ylist[0]; i++) {
3514 big[k++] = ylist[i];
3515 }
3516
3517 if (xlist != NULL && xlist[0] > 0) {
3518 big[k++] = LISTSEP;
3519 for (i=1; i<=xlist[0]; i++) {
3520 big[k++] = xlist[i];
3521 }
3522 }
3523
3524 if (rlist != NULL && rlist[0] > 0) {
3525 if (xlist == NULL || xlist[0] == 0) {
3526 /* placeholder for empty xlist */
3527 big[k++] = LISTSEP;
3528 }
3529 big[k++] = LISTSEP;
3530 for (i=1; i<=rlist[0]; i++) {
3531 big[k++] = rlist[i];
3532 }
3533 }
3534
3535 return big;
3536 }
3537
rebuild_full_VAR_list(const GRETL_VAR * var)3538 static int *rebuild_full_VAR_list (const GRETL_VAR *var)
3539 {
3540 int *list;
3541
3542 if (var->xlist == NULL && var->rlist == NULL) {
3543 list = gretl_list_copy(var->ylist);
3544 } else {
3545 list = VAR_list_composite(var->ylist, var->xlist,
3546 var->rlist);
3547 }
3548
3549 return list;
3550 }
3551
3552 /**
3553 * real_gretl_restricted_vecm:
3554 * @orig: orginal VECM model.
3555 * @rset: restriction information.
3556 * @dset: dataset struct.
3557 * @prn: gretl printing struct.
3558 * @err: location to receive error code.
3559 *
3560 * Returns: pointer to struct containing information on
3561 * the newly restricted VECM system or %NULL on failure.
3562 */
3563
3564 GRETL_VAR *
real_gretl_restricted_vecm(GRETL_VAR * orig,gretl_restriction * rset,const DATASET * dset,PRN * prn,int * err)3565 real_gretl_restricted_vecm (GRETL_VAR *orig,
3566 gretl_restriction *rset,
3567 const DATASET *dset,
3568 PRN *prn, int *err)
3569 {
3570 GRETL_VAR *jvar = NULL;
3571 gretlopt jopt = OPT_NONE;
3572 int *list = NULL;
3573
3574 if (orig == NULL || orig->jinfo == NULL || rset == NULL) {
3575 *err = E_DATA;
3576 return NULL;
3577 }
3578
3579 list = rebuild_full_VAR_list(orig);
3580 if (list == NULL) {
3581 *err = E_ALLOC;
3582 return NULL;
3583 }
3584
3585 jopt |= opt_from_jcode(orig->jinfo->code);
3586
3587 if (orig->jinfo->seasonals > 0) {
3588 jopt |= OPT_D;
3589 }
3590
3591 jvar = johansen_wrapper(VECM_ESTIMATE, orig->order + 1,
3592 orig->jinfo->rank, orig->lags,
3593 list, rset, dset,
3594 jopt, prn, err);
3595
3596 if (jvar != NULL && !jvar->err) {
3597 gretlopt ropt, prnopt = OPT_NONE;
3598 int df;
3599
3600 df = jvar->jinfo->lrdf - orig->jinfo->lrdf;
3601
3602 if (df > 0) {
3603 double x = 2 * (orig->ll - jvar->ll);
3604 double pv = chisq_cdf_comp(df, x);
3605
3606 rset_add_results(rset, x, pv, jvar->ll);
3607 rset_record_LR_result(rset);
3608 }
3609
3610 jvar->jinfo->prior_ll = orig->ll;
3611 jvar->jinfo->prior_df = orig->jinfo->lrdf;
3612
3613 ropt = gretl_restriction_get_options(rset);
3614 if (ropt & OPT_Q) {
3615 prnopt = OPT_Q;
3616 }
3617
3618 if (!(ropt & OPT_S)) {
3619 /* FIXME OPT_I, OPT_F, impulses and decomp? */
3620 gretl_VAR_print(jvar, dset, prnopt, prn);
3621 }
3622 }
3623
3624 free(list);
3625
3626 return jvar;
3627 }
3628
gretl_VAR_set_name(GRETL_VAR * var,const char * name)3629 void gretl_VAR_set_name (GRETL_VAR *var, const char *name)
3630 {
3631 if (name == var->name) {
3632 return;
3633 }
3634
3635 if (var->name == NULL) {
3636 var->name = malloc(MAXSAVENAME);
3637 }
3638
3639 if (var->name != NULL) {
3640 *var->name = '\0';
3641 strncat(var->name, name, MAXSAVENAME - 1);
3642 }
3643 }
3644
gretl_VAR_get_name(const GRETL_VAR * var)3645 const char *gretl_VAR_get_name (const GRETL_VAR *var)
3646 {
3647 return var->name;
3648 }
3649
gretl_VAR_get_resid_series(GRETL_VAR * var,int eqnum,int * err)3650 double *gretl_VAR_get_resid_series (GRETL_VAR *var, int eqnum,
3651 int *err)
3652 {
3653 MODEL *pmod;
3654 double *u = NULL;
3655
3656 if (var->models == NULL || eqnum < 0 || eqnum >= var->neqns) {
3657 *err = E_BADSTAT;
3658 return NULL;
3659 }
3660
3661 pmod = var->models[eqnum];
3662 u = copyvec(pmod->uhat, pmod->full_n);
3663
3664 if (u == NULL) {
3665 *err = E_ALLOC;
3666 }
3667
3668 return u;
3669 }
3670
gretl_VAR_set_ordering(GRETL_VAR * var,gretl_matrix * ord)3671 int gretl_VAR_set_ordering (GRETL_VAR *var, gretl_matrix *ord)
3672 {
3673 gretl_matrix_free(var->ord);
3674 var->ord = ord;
3675 return 0;
3676 }
3677
gretl_VAR_do_irf(GRETL_VAR * var,const char * line,const DATASET * dset)3678 int gretl_VAR_do_irf (GRETL_VAR *var, const char *line,
3679 const DATASET *dset)
3680 {
3681 int targ = -1, shock = 1;
3682 int h = 20, boot = 0;
3683 double alpha = 0.10;
3684 int err = 0;
3685 char *s;
3686
3687 s = strstr(line, "--targ=");
3688 if (s != NULL) {
3689 targ = atoi(s + 7) - 1;
3690 }
3691
3692 s = strstr(line, "--shock=");
3693 if (s != NULL) {
3694 shock = atoi(s + 8) - 1;
3695 }
3696
3697 s = strstr(line, "--horizon=");
3698 if (s != NULL) {
3699 h = atoi(s + 10);
3700 }
3701
3702 s = strstr(line, "--alpha=");
3703 if (s != NULL) {
3704 alpha = dot_atof(s + 8);
3705 }
3706
3707 if (strstr(line, "--bootstrap") != NULL) {
3708 boot = 1;
3709 }
3710
3711 #if 0
3712 fprintf(stderr, "targ=%d, shock=%d, h=%d, alpha=%g, boot=%d\n",
3713 targ, shock, h, alpha, boot);
3714 #endif
3715
3716 if (targ < 0 || shock < 0 || h <= 0 ||
3717 alpha < .01 || alpha > 0.5) {
3718 err = E_INVARG;
3719 } else if (boot) {
3720 err = gretl_VAR_plot_impulse_response(var, targ, shock,
3721 h, alpha, dset,
3722 OPT_NONE);
3723 } else {
3724 err = gretl_VAR_plot_impulse_response(var, targ, shock,
3725 h, 0, dset,
3726 OPT_NONE);
3727 }
3728
3729 return err;
3730 }
3731
gretl_VAR_get_highest_variable(const GRETL_VAR * var)3732 int gretl_VAR_get_highest_variable (const GRETL_VAR *var)
3733 {
3734 int i, vi, vmax = 0;
3735
3736 if (var->ylist != NULL) {
3737 for (i=1; i<=var->ylist[0]; i++) {
3738 vi = var->ylist[i];
3739 if (vi > vmax) {
3740 vmax = vi;
3741 }
3742 }
3743 }
3744
3745 if (var->xlist != NULL) {
3746 for (i=1; i<=var->xlist[0]; i++) {
3747 vi = var->xlist[i];
3748 if (vi > vmax) {
3749 vmax = vi;
3750 }
3751 }
3752 }
3753
3754 if (var->rlist != NULL) {
3755 for (i=1; i<=var->rlist[0]; i++) {
3756 vi = var->rlist[i];
3757 if (vi > vmax) {
3758 vmax = vi;
3759 }
3760 }
3761 }
3762
3763 return vmax;
3764 }
3765
gretl_VECM_id(GRETL_VAR * vecm)3766 int gretl_VECM_id (GRETL_VAR *vecm)
3767 {
3768 static int nvecm;
3769
3770 if (vecm->jinfo->ID == 0) {
3771 vecm->jinfo->ID = ++nvecm;
3772 }
3773
3774 return vecm->jinfo->ID;
3775 }
3776
gretl_VECM_n_beta(const GRETL_VAR * vecm)3777 int gretl_VECM_n_beta (const GRETL_VAR *vecm)
3778 {
3779 int nb = 0;
3780
3781 if (vecm->jinfo != NULL && vecm->jinfo->Beta != NULL) {
3782 nb = gretl_matrix_rows(vecm->jinfo->Beta);
3783 }
3784
3785 return nb;
3786 }
3787
gretl_VECM_n_alpha(const GRETL_VAR * vecm)3788 int gretl_VECM_n_alpha (const GRETL_VAR *vecm)
3789 {
3790 int na = 0;
3791
3792 if (vecm->jinfo != NULL && vecm->jinfo->Alpha != NULL) {
3793 na = gretl_matrix_rows(vecm->jinfo->Alpha);
3794 }
3795
3796 return na;
3797 }
3798
gretl_VECM_rank(const GRETL_VAR * vecm)3799 int gretl_VECM_rank (const GRETL_VAR *vecm)
3800 {
3801 int r = 0;
3802
3803 if (vecm->jinfo != NULL) {
3804 r = vecm->jinfo->rank;
3805 }
3806
3807 return r;
3808 }
3809
beta_restricted_VECM(const GRETL_VAR * vecm)3810 int beta_restricted_VECM (const GRETL_VAR *vecm)
3811 {
3812 if (vecm->jinfo != NULL && vecm->jinfo->R != NULL) {
3813 return 1;
3814 }
3815
3816 return 0;
3817 }
3818
alpha_restricted_VECM(const GRETL_VAR * vecm)3819 int alpha_restricted_VECM (const GRETL_VAR *vecm)
3820 {
3821 if (vecm->jinfo != NULL && vecm->jinfo->Ra != NULL) {
3822 return 1;
3823 }
3824
3825 return 0;
3826 }
3827
restricted_VECM(const GRETL_VAR * vecm)3828 int restricted_VECM (const GRETL_VAR *vecm)
3829 {
3830 if (vecm->jinfo != NULL &&
3831 (vecm->jinfo->R != NULL || vecm->jinfo->Ra != NULL)) {
3832 return 1;
3833 }
3834
3835 return 0;
3836 }
3837
gretl_VECM_R_matrix(const GRETL_VAR * vecm)3838 const gretl_matrix *gretl_VECM_R_matrix (const GRETL_VAR *vecm)
3839 {
3840 return (vecm->jinfo != NULL)? vecm->jinfo->R : NULL;
3841 }
3842
gretl_VECM_q_matrix(const GRETL_VAR * vecm)3843 const gretl_matrix *gretl_VECM_q_matrix (const GRETL_VAR *vecm)
3844 {
3845 return (vecm->jinfo != NULL)? vecm->jinfo->q : NULL;
3846 }
3847
gretl_VECM_Ra_matrix(const GRETL_VAR * vecm)3848 const gretl_matrix *gretl_VECM_Ra_matrix (const GRETL_VAR *vecm)
3849 {
3850 return (vecm->jinfo != NULL)? vecm->jinfo->Ra : NULL;
3851 }
3852
gretl_VECM_qa_matrix(const GRETL_VAR * vecm)3853 const gretl_matrix *gretl_VECM_qa_matrix (const GRETL_VAR *vecm)
3854 {
3855 return (vecm->jinfo != NULL)? vecm->jinfo->qa : NULL;
3856 }
3857
gretl_VAR_get_series(const GRETL_VAR * var,const DATASET * dset,int idx,const char * key,int * err)3858 double *gretl_VAR_get_series (const GRETL_VAR *var, const DATASET *dset,
3859 int idx, const char *key, int *err)
3860 {
3861 double *x = NULL;
3862 const char *msel;
3863 int t, col = 0;
3864
3865 if (var == NULL || idx != M_UHAT) { /* FIXME generalize this */
3866 *err = E_BADSTAT;
3867 return NULL;
3868 }
3869
3870 msel = strchr(key, '[');
3871 if (msel == NULL || sscanf(msel, "[,%d]", &col) != 1) {
3872 *err = E_PARSE;
3873 } else if (col <= 0 || col > var->neqns) {
3874 *err = E_DATA;
3875 }
3876
3877 if (!*err) {
3878 x = malloc(dset->n * sizeof *x);
3879 if (x == NULL) {
3880 *err = E_ALLOC;
3881 }
3882 }
3883
3884 if (!*err) {
3885 const MODEL *pmod = var->models[col-1];
3886
3887 if (pmod == NULL || pmod->full_n != dset->n) {
3888 *err = E_DATA;
3889 free(x);
3890 x = NULL;
3891 } else {
3892 for (t=0; t<dset->n; t++) {
3893 x[t] = pmod->uhat[t];
3894 }
3895 }
3896 }
3897
3898 return x;
3899 }
3900
3901 /* get coefficients or standard errors from the models
3902 in the VAR */
3903
3904 static gretl_matrix *
VAR_matrix_from_models(const GRETL_VAR * var,int idx,int * err)3905 VAR_matrix_from_models (const GRETL_VAR *var, int idx, int *err)
3906 {
3907 gretl_matrix *m = NULL;
3908 const MODEL *pmod;
3909 double x;
3910 int i, j;
3911
3912 if (idx == M_VECG) {
3913 if (var->ci != VECM || var->order == 0) {
3914 *err = E_BADSTAT;
3915 return NULL;
3916 }
3917 }
3918
3919 if (idx == M_VECG) {
3920 m = gretl_matrix_alloc(var->neqns, var->neqns * var->order);
3921 } else {
3922 m = gretl_matrix_alloc(var->models[0]->ncoeff, var->neqns);
3923 }
3924
3925 if (m == NULL) {
3926 *err = E_ALLOC;
3927 return NULL;
3928 }
3929
3930 if (idx == M_VECG) {
3931 int k, mcol, cnum;
3932
3933 for (j=0; j<var->neqns; j++) {
3934 pmod = var->models[j];
3935 mcol = 0;
3936 for (i=0; i<var->order; i++) {
3937 cnum = pmod->ifc + i;
3938 for (k=0; k<var->neqns; k++) {
3939 x = pmod->coeff[cnum];
3940 gretl_matrix_set(m, j, mcol++, x);
3941 cnum += var->order;
3942 }
3943 }
3944 }
3945 } else {
3946 for (j=0; j<var->neqns; j++) {
3947 pmod = var->models[j];
3948 for (i=0; i<pmod->ncoeff; i++) {
3949 if (idx == M_COEFF) {
3950 x = pmod->coeff[i];
3951 } else {
3952 x = pmod->sderr[i];
3953 }
3954 gretl_matrix_set(m, i, j, x);
3955 }
3956 }
3957 }
3958
3959 return m;
3960 }
3961
alt_VECM_get_EC_matrix(const GRETL_VAR * vecm,const DATASET * dset,int * err)3962 static gretl_matrix *alt_VECM_get_EC_matrix (const GRETL_VAR *vecm,
3963 const DATASET *dset,
3964 int *err)
3965 {
3966 const gretl_matrix *B = vecm->jinfo->Beta;
3967 int r = jrank(vecm);
3968 gretl_matrix *EC = NULL;
3969 double xti, xj;
3970 int s, t, T;
3971 int i, j, k;
3972
3973 if (dset == NULL || dset->Z == NULL) {
3974 *err = E_BADSTAT;
3975 return NULL;
3976 }
3977
3978 for (i=1; i<=vecm->ylist[0]; i++) {
3979 if (vecm->ylist[i] >= dset->v) {
3980 *err = E_DATA;
3981 return NULL;
3982 }
3983 }
3984
3985 T = vecm->t2 - vecm->t1 + 1;
3986
3987 EC = gretl_matrix_alloc(T, r);
3988 if (EC == NULL) {
3989 *err = E_ALLOC;
3990 return NULL;
3991 }
3992
3993 s = 0;
3994
3995 for (t=vecm->t1; t<=vecm->t2; t++) {
3996 for (j=0; j<r; j++) {
3997 xj = 0.0;
3998 /* beta * X(t-1) */
3999 k = 0;
4000 for (i=0; i<vecm->neqns; i++) {
4001 xti = dset->Z[vecm->ylist[i+1]][t-1];
4002 if (na(xti)) {
4003 xj = NADBL;
4004 break;
4005 }
4006 xj += xti * gretl_matrix_get(B, k++, j);
4007 }
4008
4009 /* restricted const or trend */
4010 if (auto_restr(vecm) && !na(xj)) {
4011 xti = gretl_matrix_get(B, k++, j);
4012 if (jcode(vecm) == J_REST_TREND) {
4013 xti *= t;
4014 }
4015 xj += xti;
4016 }
4017
4018 /* restricted exog vars */
4019 if (vecm->rlist != NULL && !na(xj)) {
4020 for (i=0; i<vecm->rlist[0]; i++) {
4021 xti = dset->Z[vecm->rlist[i+1]][t-1];
4022 if (na(xti)) {
4023 xj = NADBL;
4024 break;
4025 }
4026 xj += xti * gretl_matrix_get(B, k++, j);
4027 }
4028 }
4029
4030 if (na(xj)) {
4031 gretl_matrix_set(EC, s, j, NADBL);
4032 } else {
4033 gretl_matrix_set(EC, s, j, xj);
4034 }
4035 }
4036 s++;
4037 }
4038
4039 gretl_matrix_set_t1(EC, vecm->t1);
4040 gretl_matrix_set_t2(EC, vecm->t2);
4041
4042 return EC;
4043 }
4044
VECM_get_EC_matrix(const GRETL_VAR * v,const DATASET * dset,int * err)4045 gretl_matrix *VECM_get_EC_matrix (const GRETL_VAR *v,
4046 const DATASET *dset,
4047 int *err)
4048 {
4049 gretl_matrix *EC = NULL;
4050 double x;
4051 int rank, k, k0;
4052 int j, t, T;
4053
4054 rank = jrank(v);
4055 if (rank == 0) {
4056 *err = E_BADSTAT;
4057 return NULL;
4058 }
4059
4060 if (v->X == NULL) {
4061 fprintf(stderr, "VECM_get_EC_matrix: v->X is NULL\n");
4062 *err = E_BADSTAT;
4063 return NULL;
4064 } else if (v->X->cols < v->ncoeff) {
4065 fprintf(stderr, "VECM_get_EC_matrix: v->X is short of cols\n");
4066 return alt_VECM_get_EC_matrix(v, dset, err);
4067 }
4068
4069 T = v->X->rows;
4070
4071 EC = gretl_matrix_alloc(T, rank);
4072 if (EC == NULL) {
4073 *err = E_ALLOC;
4074 return NULL;
4075 }
4076
4077 k0 = v->ncoeff - rank;
4078
4079 for (j=0, k=k0; j<rank; j++, k++) {
4080 for (t=0; t<T; t++) {
4081 x = gretl_matrix_get(v->X, t, k);
4082 gretl_matrix_set(EC, t, j, x);
4083 }
4084 }
4085
4086 gretl_matrix_set_t1(EC, v->t1);
4087 gretl_matrix_set_t2(EC, v->t2);
4088
4089 return EC;
4090 }
4091
4092 #define vecm_matrix(i) (i == M_JALPHA || i == M_JBETA || \
4093 i == M_JVBETA || i == M_JS00 || \
4094 i == M_JS11 || i == M_JS01 || \
4095 i == M_EC || i == M_EVALS)
4096
gretl_VAR_get_matrix(const GRETL_VAR * var,int idx,int * err)4097 gretl_matrix *gretl_VAR_get_matrix (const GRETL_VAR *var, int idx,
4098 int *err)
4099 {
4100 const gretl_matrix *src = NULL;
4101 gretl_matrix *M = NULL;
4102 int copy = 1;
4103
4104 if (var == NULL) {
4105 *err = E_BADSTAT;
4106 return NULL;
4107 }
4108
4109 if (idx == M_UHAT) {
4110 src = gretl_VAR_get_residual_matrix(var);
4111 } else if (idx == M_YHAT) {
4112 M = gretl_VAR_get_fitted_matrix(var);
4113 copy = 0;
4114 } else if (idx == M_COMPAN) {
4115 src = var->A;
4116 } else if (idx == M_COEFF || idx == M_SE || idx == M_VECG) {
4117 M = VAR_matrix_from_models(var, idx, err);
4118 copy = 0;
4119 } else if (idx == M_XTXINV) {
4120 src = var->XTX;
4121 } else if (idx == M_SIGMA) {
4122 src = var->S;
4123 } else if (idx == M_PMANTEAU) {
4124 M = gretl_VAR_get_portmanteau_test(var, err);
4125 copy = 0;
4126 } else if (vecm_matrix(idx)) {
4127 if (var->jinfo != NULL) {
4128 switch (idx) {
4129 case M_EVALS:
4130 src = var->jinfo->evals;
4131 break;
4132 case M_JALPHA:
4133 src = var->jinfo->Alpha;
4134 break;
4135 case M_JBETA:
4136 src = var->jinfo->Beta;
4137 break;
4138 case M_JVBETA:
4139 src = var->jinfo->Bvar;
4140 break;
4141 case M_JS00:
4142 src = var->jinfo->S00;
4143 break;
4144 case M_JS11:
4145 src = var->jinfo->S11;
4146 break;
4147 case M_JS01:
4148 src = var->jinfo->S01;
4149 break;
4150 case M_EC:
4151 M = VECM_get_EC_matrix(var, NULL, err);
4152 copy = 0;
4153 break;
4154 }
4155 }
4156 }
4157
4158 if (copy) {
4159 if (src == NULL) {
4160 *err = E_BADSTAT;
4161 } else {
4162 M = gretl_matrix_copy(src);
4163 if (M == NULL) {
4164 *err = E_ALLOC;
4165 }
4166 }
4167 }
4168
4169 return M;
4170 }
4171
4172 /* retrieve EC (j >= 0 && j < rank) as a full-length series */
4173
gretl_VECM_get_EC(GRETL_VAR * vecm,int j,const DATASET * dset,int * err)4174 double *gretl_VECM_get_EC (GRETL_VAR *vecm, int j,
4175 const DATASET *dset, int *err)
4176 {
4177 const gretl_matrix *B = vecm->jinfo->Beta;
4178 int r = jrank(vecm);
4179 double *x = NULL;
4180 double xti;
4181 int i, k, t, t0;
4182
4183 if (j < 0 || j >= r) {
4184 *err = E_DATA;
4185 return NULL;
4186 }
4187
4188 for (i=1; i<=vecm->ylist[0]; i++) {
4189 if (vecm->ylist[i] >= dset->v) {
4190 *err = E_DATA;
4191 return NULL;
4192 }
4193 }
4194
4195 x = malloc(dset->n * sizeof *x);
4196 if (x == NULL) {
4197 *err = E_ALLOC;
4198 return NULL;
4199 }
4200
4201 t0 = (dset->t1 >= 1)? dset->t1 : 1;
4202
4203 for (t=0; t<dset->n; t++) {
4204 if (t < t0 || t > dset->t2) {
4205 x[t] = NADBL;
4206 continue;
4207 }
4208 x[t] = 0.0;
4209
4210 k = 0;
4211
4212 /* beta * X(t-1) */
4213 for (i=0; i<vecm->neqns; i++) {
4214 xti = dset->Z[vecm->ylist[i+1]][t-1];
4215 if (na(xti)) {
4216 x[t] = NADBL;
4217 break;
4218 }
4219 x[t] += xti * gretl_matrix_get(B, k++, j);
4220 }
4221
4222 /* restricted const or trend */
4223 if (auto_restr(vecm) && !na(x[t])) {
4224 xti = gretl_matrix_get(B, k++, j);
4225 if (jcode(vecm) == J_REST_TREND) {
4226 xti *= t;
4227 }
4228 x[t] += xti;
4229 }
4230
4231 /* restricted exog vars */
4232 if (vecm->rlist != NULL && !na(x[t])) {
4233 for (i=0; i<vecm->rlist[0]; i++) {
4234 xti = dset->Z[vecm->rlist[i+1]][t-1];
4235 if (na(xti)) {
4236 x[t] = NADBL;
4237 break;
4238 }
4239 x[t] += xti * gretl_matrix_get(B, k++, j);
4240 }
4241 }
4242 }
4243
4244 return x;
4245 }
4246
gretl_VAR_rebuilder_new(void)4247 static GRETL_VAR *gretl_VAR_rebuilder_new (void)
4248 {
4249 GRETL_VAR *var = malloc(sizeof *var);
4250
4251 if (var == NULL) {
4252 return NULL;
4253 }
4254
4255 gretl_VAR_clear(var);
4256
4257 return var;
4258 }
4259
VAR_retrieve_jinfo(xmlNodePtr node,xmlDocPtr doc,GRETL_VAR * var)4260 static int VAR_retrieve_jinfo (xmlNodePtr node, xmlDocPtr doc,
4261 GRETL_VAR *var)
4262 {
4263 xmlNodePtr cur = node->xmlChildrenNode;
4264 JohansenInfo *jinfo;
4265 int ID, code, rank;
4266 int seas;
4267 int got = 0;
4268 int err = 0;
4269
4270 got += gretl_xml_get_prop_as_int(node, "ID", &ID);
4271 got += gretl_xml_get_prop_as_int(node, "code", &code);
4272 got += gretl_xml_get_prop_as_int(node, "rank", &rank);
4273 got += gretl_xml_get_prop_as_int(node, "seasonals", &seas);
4274
4275 if (got != 4) {
4276 return E_DATA;
4277 }
4278
4279 jinfo = johansen_info_new(var, rank, OPT_NONE);
4280 if (jinfo == NULL) {
4281 return E_ALLOC;
4282 }
4283
4284 jinfo->ID = ID;
4285 jinfo->code = code;
4286 jinfo->rank = rank;
4287 jinfo->seasonals = seas;
4288
4289 gretl_xml_get_prop_as_double(node, "ll0", &jinfo->ll0);
4290 gretl_xml_get_prop_as_int(node, "bdf", &jinfo->lrdf);
4291 gretl_xml_get_prop_as_double(node, "oldll", &jinfo->prior_ll);
4292 gretl_xml_get_prop_as_int(node, "olddf", &jinfo->prior_df);
4293
4294 cur = node->xmlChildrenNode;
4295
4296 while (cur != NULL && !err) {
4297 if (!xmlStrcmp(cur->name, (XUC) "gretl-matrix")) {
4298 char *mname;
4299
4300 gretl_xml_get_prop_as_string(cur, "name", &mname);
4301 if (mname == NULL) {
4302 err = E_DATA;
4303 } else {
4304 if (!strcmp(mname, "u")) {
4305 jinfo->R0 = gretl_xml_get_matrix(cur, doc, &err);
4306 } else if (!strcmp(mname, "v")) {
4307 jinfo->R1 = gretl_xml_get_matrix(cur, doc, &err);
4308 } else if (!strcmp(mname, "Suu")) {
4309 jinfo->S00 = gretl_xml_get_matrix(cur, doc, &err);
4310 } else if (!strcmp(mname, "Svv")) {
4311 jinfo->S11 = gretl_xml_get_matrix(cur, doc, &err);
4312 } else if (!strcmp(mname, "Suv")) {
4313 jinfo->S01 = gretl_xml_get_matrix(cur, doc, &err);
4314 } else if (!strcmp(mname, "evals")) {
4315 jinfo->evals = gretl_xml_get_matrix(cur, doc, &err);
4316 } else if (!strcmp(mname, "Beta")) {
4317 jinfo->Beta = gretl_xml_get_matrix(cur, doc, &err);
4318 } else if (!strcmp(mname, "Alpha")) {
4319 jinfo->Alpha = gretl_xml_get_matrix(cur, doc, &err);
4320 } else if (!strcmp(mname, "Bvar")) {
4321 jinfo->Bvar = gretl_xml_get_matrix(cur, doc, &err);
4322 } else if (!strcmp(mname, "Bse")) {
4323 jinfo->Bse = gretl_xml_get_matrix(cur, doc, &err);
4324 } else if (!strcmp(mname, "R")) {
4325 jinfo->R = gretl_xml_get_matrix(cur, doc, &err);
4326 } else if (!strcmp(mname, "q")) {
4327 jinfo->q = gretl_xml_get_matrix(cur, doc, &err);
4328 } else if (!strcmp(mname, "Ra")) {
4329 jinfo->Ra = gretl_xml_get_matrix(cur, doc, &err);
4330 } else if (!strcmp(mname, "qa")) {
4331 jinfo->qa = gretl_xml_get_matrix(cur, doc, &err);
4332 }
4333 free(mname);
4334 }
4335 }
4336 cur = cur->next;
4337 }
4338
4339 if (err) {
4340 johansen_info_free(jinfo);
4341 } else {
4342 var->jinfo = jinfo;
4343 }
4344
4345 fprintf(stderr, "VAR_retrieve_jinfo: err = %d\n", err);
4346
4347 return err;
4348 }
4349
VAR_retrieve_equations(xmlNodePtr node,xmlDocPtr doc,MODEL ** models,int neqns,const DATASET * dset)4350 static int VAR_retrieve_equations (xmlNodePtr node, xmlDocPtr doc,
4351 MODEL **models, int neqns,
4352 const DATASET *dset)
4353 {
4354 xmlNodePtr cur = node->xmlChildrenNode;
4355 int i = 0, err = 0;
4356
4357 while (cur != NULL && !err) {
4358 MODEL *pmod;
4359
4360 if (!xmlStrcmp(cur->name, (XUC) "gretl-model")) {
4361 pmod = gretl_model_from_XML(cur, doc, dset, &err);
4362 if (pmod != NULL) {
4363 models[i++] = pmod;
4364 }
4365 }
4366 cur = cur->next;
4367 }
4368
4369 fprintf(stderr, "VAR_retrieve_equations: got %d models (neqns = %d), err = %d\n",
4370 i, neqns, err);
4371
4372 return err;
4373 }
4374
4375 /* for backward compatibility */
4376
make_VAR_global_lists(GRETL_VAR * var)4377 static int make_VAR_global_lists (GRETL_VAR *var)
4378 {
4379 MODEL *pmod = var->models[0];
4380 int *list = pmod->list;
4381 int p = var->order;
4382 int n = var->neqns;
4383 int nx, np = n * p;
4384 int i, j, ifc;
4385
4386 if (list == NULL || list[0] < 3) {
4387 return E_DATA;
4388 }
4389
4390 /* the _endogenous_ vars start in position 2 or 3 (3 if a constant
4391 is included), and there are (order * neqns) such terms */
4392
4393 ifc = (list[2] == 0);
4394 nx = list[0] - 1 - ifc - np;
4395
4396 var->ylist = gretl_list_new(n);
4397 if (var->ylist == NULL) {
4398 return E_ALLOC;
4399 }
4400
4401 if (nx > 0) {
4402 var->xlist = gretl_list_new(nx);
4403 if (var->xlist == NULL) {
4404 free(var->ylist);
4405 var->ylist = NULL;
4406 return E_ALLOC;
4407 }
4408 }
4409
4410 for (i=0; i<n; i++) {
4411 pmod = var->models[i];
4412 var->ylist[i+1] = pmod->list[1];
4413 }
4414
4415 j = 2 + ifc + np;
4416 for (i=0; i<nx; i++) {
4417 var->xlist[i+1] = list[j++];
4418 }
4419
4420 return 0;
4421 }
4422
rebuild_VAR_matrices(GRETL_VAR * var)4423 static int rebuild_VAR_matrices (GRETL_VAR *var)
4424 {
4425 MODEL *pmod;
4426 double x;
4427 int gotA = (var->A != NULL);
4428 int gotX = (var->X != NULL);
4429 int j, i;
4430 int err = 0;
4431
4432 if (var->E == NULL) {
4433 err = VAR_allocate_residuals_matrix(var);
4434 }
4435
4436 if (!err && var->A == NULL) {
4437 err = VAR_allocate_companion_matrix(var);
4438 }
4439
4440 if (!err && gotX && var->XTX == NULL) {
4441 err = VAR_add_XTX_matrix(var);
4442 }
4443
4444 if (!err && var->C == NULL) {
4445 err = VAR_allocate_cholesky_matrix(var);
4446 }
4447
4448 if (!err && var->B == NULL) {
4449 var->B = gretl_matrix_alloc(var->models[0]->ncoeff,
4450 var->neqns);
4451 if (var->B == NULL) {
4452 err = E_ALLOC;
4453 }
4454 }
4455
4456 for (j=0; j<var->neqns && !err; j++) {
4457 pmod = var->models[j];
4458 for (i=0; i<pmod->ncoeff; i++) {
4459 x = pmod->coeff[i];
4460 gretl_matrix_set(var->B, i, j, x);
4461 }
4462 for (i=0; i<var->T; i++) {
4463 x = pmod->uhat[pmod->t1 + i];
4464 gretl_matrix_set(var->E, i, j, x);
4465 }
4466 }
4467
4468 if (!err && !gotA) {
4469 /* note: for VECMs, A should be retrieved from the session
4470 file, and gotA should be non-zero
4471 */
4472 VAR_write_A_matrix(var, GRETL_MOD_NONE);
4473 }
4474
4475 if (!err && !gotX) {
4476 fprintf(stderr, "Can't we rebuild VAR->X somehow?\n");
4477 }
4478
4479 return err;
4480 }
4481
gretl_VAR_from_XML(xmlNodePtr node,xmlDocPtr doc,const DATASET * dset,int * err)4482 GRETL_VAR *gretl_VAR_from_XML (xmlNodePtr node, xmlDocPtr doc,
4483 const DATASET *dset,
4484 int *err)
4485 {
4486 GRETL_VAR *var;
4487 MODEL *pmod;
4488 xmlNodePtr cur;
4489 char *vname;
4490 int i, n, got = 0;
4491
4492 var = gretl_VAR_rebuilder_new();
4493 if (var == NULL) {
4494 *err = E_ALLOC;
4495 return NULL;
4496 }
4497
4498 got += gretl_xml_get_prop_as_int(node, "ecm", &var->ci);
4499 got += gretl_xml_get_prop_as_int(node, "neqns", &var->neqns);
4500 got += gretl_xml_get_prop_as_int(node, "order", &var->order);
4501
4502 if (got < 3) {
4503 *err = E_DATA;
4504 goto bailout;
4505 }
4506
4507 var->ci = (var->ci == 0)? VAR : VECM;
4508
4509 gretl_xml_get_prop_as_string(node, "name", &vname);
4510 if (vname != NULL) {
4511 gretl_VAR_set_name(var, vname);
4512 free(vname);
4513 }
4514
4515 /* these are not show-stoppers */
4516 gretl_xml_get_prop_as_int(node, "robust", &var->robust);
4517 gretl_xml_get_prop_as_int(node, "detflags", &var->detflags);
4518 gretl_xml_get_prop_as_int(node, "LBs", &var->LBs);
4519 gretl_xml_get_prop_as_double(node, "LB", &var->LB);
4520
4521 var->models = malloc(var->neqns * sizeof *var->models);
4522 if (var->models == NULL) {
4523 *err = E_ALLOC;
4524 goto bailout;
4525 }
4526
4527 for (i=0; i<var->neqns; i++) {
4528 var->models[i] = NULL;
4529 }
4530
4531 gretl_push_c_numeric_locale();
4532
4533 cur = node->xmlChildrenNode;
4534
4535 while (cur != NULL && !*err) {
4536 if (!xmlStrcmp(cur->name, (XUC) "lags")) {
4537 var->lags = gretl_xml_get_list(cur, doc, err);
4538 } else if (!xmlStrcmp(cur->name, (XUC) "ylist")) {
4539 var->ylist = gretl_xml_get_list(cur, doc, err);
4540 } else if (!xmlStrcmp(cur->name, (XUC) "xlist")) {
4541 var->xlist = gretl_xml_get_list(cur, doc, err);
4542 } else if (!xmlStrcmp(cur->name, (XUC) "rlist")) {
4543 var->rlist = gretl_xml_get_list(cur, doc, err);
4544 } else if (!xmlStrcmp(cur->name, (XUC) "Fvals")) {
4545 var->Fvals = gretl_xml_get_double_array(cur, doc, &n, err);
4546 } else if (!xmlStrcmp(cur->name, (XUC) "Ivals")) {
4547 var->Ivals = gretl_xml_get_double_array(cur, doc, &n, err);
4548 } else if (!xmlStrcmp(cur->name, (XUC) "equations")) {
4549 *err = VAR_retrieve_equations(cur, doc, var->models, var->neqns, dset);
4550 } else if (!xmlStrcmp(cur->name, (XUC) "gretl-johansen")) {
4551 *err = VAR_retrieve_jinfo(cur, doc, var);
4552 } else if (!xmlStrcmp(cur->name, (XUC) "gretl-matrix")) {
4553 char *mname;
4554
4555 gretl_xml_get_prop_as_string(cur, "name", &mname);
4556 if (mname == NULL) {
4557 *err = E_DATA;
4558 } else {
4559 if (!strcmp(mname, "A")) {
4560 var->A = gretl_xml_get_matrix(cur, doc, err);
4561 } else if (!strcmp(mname, "X")) {
4562 var->X = gretl_xml_get_matrix(cur, doc, err);
4563 } else if (!strcmp(mname, "Y")) {
4564 var->Y = gretl_xml_get_matrix(cur, doc, err);
4565 } else if (!strcmp(mname, "ord")) {
4566 var->ord = gretl_xml_get_matrix(cur, doc, err);
4567 }
4568 free(mname);
4569 }
4570 }
4571 cur = cur->next;
4572 }
4573
4574 gretl_pop_c_numeric_locale();
4575
4576 if (*err) {
4577 goto bailout;
4578 }
4579
4580 pmod = var->models[0];
4581
4582 /* Note: get "global" info from the first model (equation) */
4583 var->ncoeff = pmod->ncoeff;
4584 var->ifc = pmod->ifc;
4585 var->t1 = pmod->t1;
4586 var->t2 = pmod->t2;
4587 var->T = var->t2 - var->t1 + 1;
4588 var->df = var->T - var->ncoeff;
4589
4590 if (var->ylist == NULL) {
4591 *err = make_VAR_global_lists(var);
4592 }
4593
4594 if (!*err) {
4595 *err = rebuild_VAR_matrices(var);
4596 }
4597
4598 if (!*err) {
4599 *err = VAR_add_stats(var, 0);
4600 }
4601
4602 if (!*err) {
4603 *err = gretl_VAR_do_error_decomp(var->S, var->C, NULL);
4604 }
4605
4606 /* FIXME vecm tests on beta/alpha */
4607
4608 bailout:
4609
4610 if (*err) {
4611 gretl_VAR_free(var);
4612 var = NULL;
4613 }
4614
4615 return var;
4616 }
4617
johansen_serialize(JohansenInfo * j,PRN * prn)4618 static void johansen_serialize (JohansenInfo *j, PRN *prn)
4619 {
4620 pprintf(prn, "<gretl-johansen ID=\"%d\" code=\"%d\" rank=\"%d\" ",
4621 j->ID, j->code, j->rank);
4622 pprintf(prn, "seasonals=\"%d\" ", j->seasonals);
4623
4624 if (j->lrdf > 0 && !na(j->ll0)) {
4625 gretl_xml_put_double("ll0", j->ll0, prn);
4626 gretl_xml_put_int("bdf", j->lrdf, prn);
4627 }
4628
4629 if (j->prior_df > 0 && !na(j->prior_ll)) {
4630 gretl_xml_put_double("oldll", j->prior_ll, prn);
4631 gretl_xml_put_int("olddf", j->prior_df, prn);
4632 }
4633
4634 pputs(prn, ">\n");
4635
4636 gretl_matrix_serialize(j->R0, "u", prn);
4637 gretl_matrix_serialize(j->R1, "v", prn);
4638 gretl_matrix_serialize(j->S00, "Suu", prn);
4639 gretl_matrix_serialize(j->S11, "Svv", prn);
4640 gretl_matrix_serialize(j->S01, "Suv", prn);
4641 gretl_matrix_serialize(j->evals, "evals", prn);
4642 gretl_matrix_serialize(j->Beta, "Beta", prn);
4643 gretl_matrix_serialize(j->Alpha, "Alpha", prn);
4644 gretl_matrix_serialize(j->Bvar, "Bvar", prn);
4645 gretl_matrix_serialize(j->Bse, "Bse", prn);
4646 gretl_matrix_serialize(j->R, "R", prn);
4647 gretl_matrix_serialize(j->q, "q", prn);
4648 gretl_matrix_serialize(j->Ra, "Ra", prn);
4649 gretl_matrix_serialize(j->qa, "qa", prn);
4650
4651 pputs(prn, "</gretl-johansen>\n");
4652 }
4653
4654 /* Retrieve enough VECM-related info from @b to carry
4655 out an IRF bootstrap: this requires more testing.
4656 */
4657
retrieve_johansen_basics(GRETL_VAR * var,gretl_bundle * b)4658 static int retrieve_johansen_basics (GRETL_VAR *var,
4659 gretl_bundle *b)
4660 {
4661 int err = 0;
4662
4663 var->jinfo = calloc(1, sizeof *var->jinfo);
4664
4665 if (var->jinfo == NULL) {
4666 err = E_ALLOC;
4667 } else {
4668 int i, e[10] = {0};
4669
4670 var->jinfo->code = gretl_bundle_get_int(b, "code", &e[0]);
4671 var->jinfo->rank = gretl_bundle_get_int(b, "rank", &e[1]);
4672 var->jinfo->seasonals = gretl_bundle_get_int(b, "seasonals", &e[2]);
4673 var->jinfo->R0 = gretl_bundle_get_matrix(b, "u", &e[3]);
4674 var->jinfo->R1 = gretl_bundle_get_matrix(b, "v", &e[4]);
4675 var->jinfo->S00 = gretl_bundle_get_matrix(b, "Suu", &e[5]);
4676 var->jinfo->S11 = gretl_bundle_get_matrix(b, "Svv", &e[6]);
4677 var->jinfo->S01 = gretl_bundle_get_matrix(b, "Suv", &e[7]);
4678 var->jinfo->Beta = gretl_bundle_get_matrix(b, "Beta", &e[8]);
4679 var->jinfo->Alpha = gretl_bundle_get_matrix(b, "Alpha", &e[9]);
4680
4681 for (i=0; i<10; i++) {
4682 if (e[i]) {
4683 err = e[i];
4684 break;
4685 }
4686 }
4687
4688 if (!err) {
4689 /* restriction matrices: may not be present */
4690 var->jinfo->R = gretl_bundle_get_matrix(b, "R", NULL);
4691 var->jinfo->q = gretl_bundle_get_matrix(b, "q", NULL);
4692 var->jinfo->Ra = gretl_bundle_get_matrix(b, "Ra", NULL);
4693 var->jinfo->qa = gretl_bundle_get_matrix(b, "qa", NULL);
4694 }
4695 }
4696
4697 if (err && var->jinfo != NULL) {
4698 free(var->jinfo);
4699 var->jinfo = NULL;
4700 }
4701
4702 return err;
4703 }
4704
johansen_bundlize(JohansenInfo * j)4705 static gretl_bundle *johansen_bundlize (JohansenInfo *j)
4706 {
4707 gretl_bundle *b = gretl_bundle_new();
4708
4709 if (b == NULL) {
4710 return NULL;
4711 }
4712
4713 /* "j->ID"? */
4714 gretl_bundle_set_int(b, "code", j->code);
4715 gretl_bundle_set_int(b, "rank", j->rank);
4716 gretl_bundle_set_int(b, "seasonals", j->seasonals);
4717
4718 if (j->lrdf > 0 && !na(j->ll0)) {
4719 gretl_bundle_set_scalar(b, "ll0", j->ll0);
4720 gretl_bundle_set_scalar(b, "bdf", j->lrdf);
4721 }
4722
4723 gretl_bundle_set_matrix(b, "u", j->R0);
4724 gretl_bundle_set_matrix(b, "v", j->R1);
4725 gretl_bundle_set_matrix(b, "Suu", j->S00);
4726 gretl_bundle_set_matrix(b, "Svv", j->S11);
4727 gretl_bundle_set_matrix(b, "Suv", j->S01);
4728 gretl_bundle_set_matrix(b, "evals", j->evals);
4729 gretl_bundle_set_matrix(b, "Beta", j->Beta);
4730 gretl_bundle_set_matrix(b, "Alpha", j->Alpha);
4731 gretl_bundle_set_matrix(b, "Bvar", j->Bvar);
4732 gretl_bundle_set_matrix(b, "Bse", j->Bse);
4733
4734 if (j->R != NULL) {
4735 gretl_bundle_set_matrix(b, "R", j->R);
4736 }
4737 if (j->q != NULL) {
4738 gretl_bundle_set_matrix(b, "q", j->q);
4739 }
4740 if (j->Ra != NULL) {
4741 gretl_bundle_set_matrix(b, "Ra", j->Ra);
4742 }
4743 if (j->qa != NULL) {
4744 gretl_bundle_set_matrix(b, "qa", j->qa);
4745 }
4746
4747 return b;
4748 }
4749
gretl_VAR_serialize(const GRETL_VAR * var,SavedObjectFlags flags,PRN * prn)4750 int gretl_VAR_serialize (const GRETL_VAR *var, SavedObjectFlags flags,
4751 PRN *prn)
4752 {
4753 char *xmlname = NULL;
4754 int g = var->neqns;
4755 int m = g * g + g;
4756 int i, err = 0;
4757
4758 if (var->name == NULL || *var->name == '\0') {
4759 xmlname = gretl_strdup("none");
4760 } else {
4761 xmlname = gretl_xml_encode(var->name);
4762 }
4763
4764 pprintf(prn, "<gretl-VAR name=\"%s\" saveflags=\"%d\" ", xmlname, (int) flags);
4765 free(xmlname);
4766
4767 pprintf(prn, "ecm=\"%d\" neqns=\"%d\" order=\"%d\" detflags=\"%d\" ",
4768 (var->ci == VECM), var->neqns, var->order, var->detflags);
4769
4770 if (var->robust) {
4771 gretl_xml_put_int("robust", var->robust, prn);
4772 }
4773
4774 if (var->LBs > 0 && !na(var->LB)) {
4775 /* Portmanteau test */
4776 gretl_xml_put_double("LB", var->LB, prn);
4777 gretl_xml_put_int("LBs", var->LBs, prn);
4778 }
4779
4780 pputs(prn, ">\n");
4781
4782 gretl_xml_put_tagged_list("lags", var->lags, prn);
4783 gretl_xml_put_tagged_list("ylist", var->ylist, prn);
4784 gretl_xml_put_tagged_list("xlist", var->xlist, prn);
4785 gretl_xml_put_tagged_list("rlist", var->rlist, prn);
4786
4787 gretl_push_c_numeric_locale();
4788
4789 if (var->Fvals != NULL) {
4790 gretl_xml_put_double_array("Fvals", var->Fvals, m, prn);
4791 }
4792
4793 if (var->Ivals != NULL) {
4794 gretl_xml_put_double_array("Ivals", var->Ivals, N_IVALS, prn);
4795 }
4796
4797 if (var->X != NULL && var->Y != NULL) {
4798 /* could be fiddly to reconstruct, needed for IRF bootstrap */
4799 gretl_matrix_serialize(var->X, "X", prn);
4800 gretl_matrix_serialize(var->Y, "Y", prn);
4801 }
4802
4803 if (var->ord != NULL) {
4804 gretl_matrix_serialize(var->ord, "ord", prn);
4805 }
4806
4807 if (var->ci == VECM) {
4808 /* this is hard to reconstruct for VECMs */
4809 gretl_matrix_serialize(var->A, "A", prn);
4810 }
4811
4812 gretl_pop_c_numeric_locale();
4813
4814 pputs(prn, "<equations>\n");
4815
4816 for (i=0; i<var->neqns; i++) {
4817 gretl_model_serialize(var->models[i], 0, prn);
4818 }
4819
4820 pputs(prn, "</equations>\n");
4821
4822 if (var->jinfo != NULL) {
4823 johansen_serialize(var->jinfo, prn);
4824 }
4825
4826 pputs(prn, "</gretl-VAR>\n");
4827
4828 return err;
4829 }
4830
4831 /* Clean-up of temporary "jinfo" used in IRF
4832 bootstrapping. We can't use the regular destructor,
4833 johansen_info_free(), because the primary matrix
4834 pointers on var->jinfo are in this case borrowed
4835 from a $system bundle; but we do need to free any
4836 "extra" matrices that have been added.
4837 */
4838
free_temp_jinfo(GRETL_VAR * var)4839 static void free_temp_jinfo (GRETL_VAR *var)
4840 {
4841 /* handle extra matrices */
4842 gretl_matrix_free(var->jinfo->YY);
4843 gretl_matrix_free(var->jinfo->RR);
4844 gretl_matrix_free(var->jinfo->BB);
4845
4846 free(var->jinfo);
4847 }
4848
4849 /* When destroying a VAR that was reconstituted by reading
4850 from a $system bundle, we must be careful not to free
4851 any elements which were just borrowed from the bundle;
4852 we set the associated pointers to NULL.
4853 */
4854
destroy_VAR_from_bundle(GRETL_VAR * var)4855 static void destroy_VAR_from_bundle (GRETL_VAR *var)
4856 {
4857 gretl_matrix **mm[] = {
4858 &var->A, &var->C, &var->X, &var->Y,
4859 &var->XTX, &var->B, &var->S, &var->E,
4860 NULL
4861 };
4862 int i, imin = 0;
4863
4864 /* var->A is NOT borrowed */
4865 imin = 1;
4866
4867 /* nullify borrowed matrices */
4868 for (i=imin; mm[i] != NULL; i++) {
4869 *mm[i] = NULL;
4870 }
4871
4872 /* nullify borrowed lists */
4873 var->xlist = var->ylist = var->rlist = NULL;
4874
4875 if (var->jinfo != NULL) {
4876 free_temp_jinfo(var);
4877 var->jinfo = NULL;
4878 }
4879
4880 gretl_VAR_free(var);
4881 }
4882
4883 /* Reconstitute a VAR or VECM (partially) from the bundle @b,
4884 for the purposes of generating an FEVD or IRF. How far we
4885 need to go in rebuilding the GRETL_VAR struct depends on
4886 the purpose, and in the IRF case, whether bootstrapping is
4887 required. For FEVD we basically just need the A and C
4888 matrices from the bundle.
4889 */
4890
build_VAR_from_bundle(gretl_bundle * b,int irf,int boot,int * err)4891 static GRETL_VAR *build_VAR_from_bundle (gretl_bundle *b,
4892 int irf, int boot,
4893 int *err)
4894 {
4895 GRETL_VAR *var = malloc(sizeof *var);
4896 int i, ierr[7];
4897
4898 if (var == NULL) {
4899 *err = E_ALLOC;
4900 return NULL;
4901 }
4902
4903 gretl_VAR_clear(var);
4904
4905 if (gretl_bundle_get_int(b, "ecm", NULL)) {
4906 var->ci = VECM;
4907 } else {
4908 var->ci = VAR;
4909 }
4910
4911 /* get integer dimensions */
4912 var->neqns = gretl_bundle_get_int(b, "neqns", &ierr[0]);
4913 var->order = gretl_bundle_get_int(b, "order", &ierr[1]);
4914 var->ncoeff = gretl_bundle_get_int(b, "ncoeff", &ierr[2]);
4915 var->t1 = gretl_bundle_get_int(b, "t1", &ierr[3]);
4916 var->t2 = gretl_bundle_get_int(b, "t2", &ierr[4]);
4917 var->df = gretl_bundle_get_int(b, "df", &ierr[6]);
4918 var->T = gretl_bundle_get_int(b, "T", &ierr[5]);
4919 for (i=0; i<7; i++) {
4920 if (ierr[i]) {
4921 *err = ierr[i];
4922 break;
4923 }
4924 }
4925
4926 if (!*err) {
4927 /* convert sample range to zero-based */
4928 var->t1 -= 1;
4929 var->t2 -= 1;
4930 var->ifc = gretl_bundle_get_int(b, "ifc", NULL);
4931 /* note: borrowing of lists from bundle */
4932 var->ylist = gretl_bundle_get_list(b, "ylist", NULL);
4933 var->xlist = gretl_bundle_get_list(b, "xlist", NULL);
4934 var->rlist = gretl_bundle_get_list(b, "rlist", NULL);
4935 }
4936
4937 if (!*err) {
4938 /* note: borrowing of matrices */
4939 gretl_matrix **mm[] = {
4940 &var->A, &var->C, &var->X, &var->Y,
4941 &var->XTX, &var->B, &var->S, &var->E
4942 };
4943 const char *keys[] = {
4944 "A", "C", "X", "Y", "xtxinv",
4945 "coeff", "sigma", "uhat"
4946 };
4947 int n = 2 + 2*irf + 4*boot;
4948
4949 for (i=0; i<n && !*err; i++) {
4950 *mm[i] = gretl_bundle_get_matrix(b, keys[i], err);
4951 }
4952 if (!*err) {
4953 /* matrix @A from the bundle must be converted into
4954 the full companion matrix
4955 */
4956 gretl_matrix *A = companionize(var->A, err);
4957
4958 var->A = A;
4959 }
4960 if (!*err && var->ci == VECM && boot) {
4961 /* not fully ready, won't be reached at present? */
4962 gretl_bundle *jb;
4963
4964 jb = gretl_bundle_get_bundle(b, "vecm_info", err);
4965 if (!*err) {
4966 *err = retrieve_johansen_basics(var, jb);
4967 }
4968 }
4969 }
4970
4971 if (!*err && var->neqns == 0) {
4972 /* maybe we have a minimal "home-made" bundle? */
4973 if (var->C != NULL && var->C->cols > 0) {
4974 var->neqns = var->C->cols;
4975 } else {
4976 *err = E_INVARG;
4977 }
4978 }
4979
4980 if (*err) {
4981 destroy_VAR_from_bundle(var);
4982 var = NULL;
4983 }
4984
4985 return var;
4986 }
4987
gretl_FEVD_from_bundle(gretl_bundle * b,int targ,int shock,const DATASET * dset,int * err)4988 gretl_matrix *gretl_FEVD_from_bundle (gretl_bundle *b,
4989 int targ, int shock,
4990 const DATASET *dset,
4991 int *err)
4992 {
4993 GRETL_VAR *var = build_VAR_from_bundle(b, 0, 0, err);
4994 gretl_matrix *ret = NULL;
4995
4996 if (var != NULL) {
4997 ret = gretl_VAR_get_FEVD_matrix(var, targ, shock, 0, dset, err);
4998 destroy_VAR_from_bundle(var);
4999 }
5000
5001 return ret;
5002 }
5003
gretl_IRF_from_bundle(gretl_bundle * b,int targ,int shock,double alpha,const DATASET * dset,int * err)5004 gretl_matrix *gretl_IRF_from_bundle (gretl_bundle *b,
5005 int targ, int shock,
5006 double alpha,
5007 const DATASET *dset,
5008 int *err)
5009 {
5010 gretl_matrix *ret = NULL;
5011 GRETL_VAR *var = NULL;
5012 int boot = 0;
5013
5014 if (alpha != 0 && (alpha < 0.01 || alpha > 0.6)) {
5015 *err = E_INVARG;
5016 return NULL;
5017 } else if (alpha != 0) {
5018 boot = 1;
5019 }
5020
5021 var = build_VAR_from_bundle(b, 1, boot, err);
5022
5023 if (!*err && var != NULL) {
5024 ret = gretl_VAR_get_impulse_response(var, targ, shock,
5025 0, alpha, dset, err);
5026 }
5027
5028 if (var != NULL) {
5029 destroy_VAR_from_bundle(var);
5030 }
5031
5032 return ret;
5033 }
5034
make_detflags_matrix(const GRETL_VAR * var)5035 static gretl_matrix *make_detflags_matrix (const GRETL_VAR *var)
5036 {
5037 gretl_matrix *d = gretl_zero_matrix_new(1, 3);
5038
5039 if (var->ifc || (var->detflags & DET_CONST)) {
5040 d->val[0] = 1;
5041 }
5042 if (var->detflags & DET_TREND) {
5043 d->val[1] = 1;
5044 }
5045 if (var->detflags & DET_SEAS) {
5046 d->val[2] = 1;
5047 }
5048
5049 return d;
5050 }
5051
gretl_VAR_bundlize(const GRETL_VAR * var,DATASET * dset,gretl_bundle * b)5052 int gretl_VAR_bundlize (const GRETL_VAR *var,
5053 DATASET *dset,
5054 gretl_bundle *b)
5055 {
5056 gretl_matrix *d;
5057 int err = 0;
5058
5059 gretl_bundle_set_string(b, "command", gretl_command_word(var->ci));
5060
5061 if (var->name != NULL) {
5062 gretl_bundle_set_string(b, "name", var->name);
5063 }
5064 gretl_bundle_set_int(b, "ecm", var->ci == VECM);
5065 gretl_bundle_set_int(b, "neqns", var->neqns);
5066 gretl_bundle_set_int(b, "ncoeff", var->ncoeff);
5067 gretl_bundle_set_int(b, "order", var->order);
5068 gretl_bundle_set_int(b, "robust", var->robust);
5069 gretl_bundle_set_int(b, "t1", var->t1 + 1);
5070 gretl_bundle_set_int(b, "t2", var->t2 + 1);
5071 gretl_bundle_set_int(b, "df", var->df);
5072 gretl_bundle_set_int(b, "T", var->T);
5073 gretl_bundle_set_int(b, "ifc", var->ifc);
5074
5075 if (var->models != NULL && var->models[0] != NULL) {
5076 gretl_bundle_set_int(b, "sample_t1", var->models[0]->smpl.t1 + 1);
5077 gretl_bundle_set_int(b, "sample_t2", var->models[0]->smpl.t2 + 1);
5078 }
5079
5080 gretl_bundle_set_scalar(b, "lnl", var->ll);
5081 gretl_bundle_set_scalar(b, "ldet", var->ldet);
5082
5083 if (var->LBs > 0 && !na(var->LB)) {
5084 /* Portmanteau test */
5085 gretl_bundle_set_scalar(b, "Ljung_Box", var->LB);
5086 gretl_bundle_set_scalar(b, "LB_order", var->LBs);
5087 }
5088
5089 /* deterministic terms: pack in matrix */
5090 d = make_detflags_matrix(var);
5091 gretl_bundle_donate_data(b, "detflags", d,
5092 GRETL_TYPE_MATRIX, 0);
5093
5094 /* lists: lags, ylist, xlist, rlist */
5095 if (var->lags != NULL) {
5096 gretl_matrix *v = gretl_list_to_vector(var->lags, &err);
5097
5098 if (!err) {
5099 gretl_bundle_donate_data(b, "lags", v,
5100 GRETL_TYPE_MATRIX, 0);
5101 gretl_bundle_set_note(b, "lags", "gappy lags vector");
5102 }
5103 }
5104 if (var->ylist != NULL) {
5105 gretl_bundle_set_list(b, "ylist", var->ylist);
5106 }
5107 if (var->xlist != NULL) {
5108 gretl_bundle_set_list(b, "xlist", var->xlist);
5109 }
5110 if (var->rlist != NULL) {
5111 gretl_bundle_set_list(b, "rlist", var->rlist);
5112 }
5113
5114 /* doubles arrays: Fvals, Ivals? */
5115
5116 if (var->A != NULL) {
5117 /* companion matrix */
5118 if (var->A->rows == var->A->cols) {
5119 gretl_bundle_set_matrix(b, "A", var->A);
5120 } else {
5121 /* A needs trimming */
5122 gretl_matrix *A = decompanionize(var->A, var->neqns,
5123 GRETL_MOD_NONE);
5124
5125 if (A != NULL) {
5126 gretl_bundle_donate_data(b, "A", A, GRETL_TYPE_MATRIX, 0);
5127 }
5128 }
5129 }
5130 if (var->C != NULL) {
5131 gretl_bundle_set_matrix(b, "C", var->C);
5132 }
5133 if (var->B != NULL) {
5134 gretl_bundle_set_matrix(b, "coeff", var->B);
5135 }
5136 if (var->S != NULL) {
5137 gretl_bundle_set_matrix(b, "sigma", var->S);
5138 }
5139 if (var->XTX != NULL) {
5140 gretl_bundle_set_matrix(b, "xtxinv", var->XTX);
5141 }
5142 if (var->E != NULL) {
5143 gretl_bundle_set_matrix(b, "uhat", var->E);
5144 }
5145 if (var->X != NULL && var->Y != NULL) {
5146 gretl_bundle_set_matrix(b, "X", var->X);
5147 gretl_bundle_set_matrix(b, "Y", var->Y);
5148 }
5149 if (var->ord != NULL) {
5150 gretl_bundle_set_matrix(b, "ord", var->ord);
5151 }
5152
5153 #if 0
5154 /* bundle up the individual MODEL structs (?) */
5155 gretl_array *models;
5156 int i;
5157
5158 models = gretl_array_new(GRETL_TYPE_BUNDLES, var->neqns, &err);
5159 for (i=0; i<var->neqns && !err; i++) {
5160 gretl_bundle *mbi;
5161
5162 mbi = bundle_from_model(var->models[i], dset, &err);
5163 if (!err) {
5164 gretl_array_set_bundle(models, i, mbi, 0);
5165 }
5166 }
5167 if (!err) {
5168 err = gretl_bundle_set_data(b, "models", models,
5169 GRETL_TYPE_ARRAY, 0);
5170 } else if (models != NULL) {
5171 gretl_array_destroy(models);
5172 }
5173 #endif
5174
5175 if (var->jinfo != NULL) {
5176 /* VECM specific info */
5177 gretl_bundle *jb = johansen_bundlize(var->jinfo);
5178
5179 if (jb != NULL) {
5180 err = gretl_bundle_donate_data(b, "vecm_info", jb,
5181 GRETL_TYPE_BUNDLE, 0);
5182 }
5183 }
5184
5185 return err;
5186 }
5187