"Fossies" - the Fresh Open Source Software Archive

Member "gretl-2020e/lib/src/var.c" (21 Nov 2020, 107127 Bytes) of package /linux/misc/gretl-2020e.tar.xz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "var.c" see the Fossies "Dox" file reference documentation and the last Fossies "Diffs" side-by-side code changes report: 2020d_vs_2020e.

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