"Fossies" - the Fresh Open Source Software Archive

Member "gretl-2020e/plugin/sysest.c" (7 Oct 2019, 30830 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 "sysest.c" see the Fossies "Dox" file reference documentation.

    1 /*
    2  *  gretl -- Gnu Regression, Econometrics and Time-series Library
    3  *  Copyright (C) 2001 Allin Cottrell and Riccardo "Jack" Lucchetti
    4  *
    5  *  This program is free software: you can redistribute it and/or modify
    6  *  it under the terms of the GNU General Public License as published by
    7  *  the Free Software Foundation, either version 3 of the License, or
    8  *  (at your option) any later version.
    9  *
   10  *  This program is distributed in the hope that it will be useful,
   11  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
   12  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   13  *  GNU General Public License for more details.
   14  *
   15  *  You should have received a copy of the GNU General Public License
   16  *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
   17  *
   18  */
   19 
   20 #include "libgretl.h"
   21 #include "version.h"
   22 #include "gretl_matrix.h"
   23 #include "libset.h"
   24 #include "system.h"
   25 #include "tsls.h"
   26 #include "sysml.h"
   27 
   28 #define SDEBUG 0
   29 
   30 #define sys_ols_ok(s) (s->method == SYS_METHOD_SUR || \
   31                        s->method == SYS_METHOD_OLS || \
   32                        s->method == SYS_METHOD_WLS)
   33 
   34 /* Insert the elements of sub-matrix @M, multiplied by @scale, in the
   35    appropriate position within the big matrix @X. We're building
   36    a symmetric matrix here so we insert off-diagonal elements both
   37    above and below the diagonal.
   38 */
   39 
   40 static void
   41 insert_sys_X_block (gretl_matrix *X, const gretl_matrix *M,
   42             int startrow, int startcol, double scale)
   43 {
   44     int r, c, i, j;
   45     double mij;
   46 
   47     for (i=0; i<M->rows; i++) {
   48     r = startrow + i;
   49     for (j=0; j<M->cols; j++) {
   50         c = startcol + j;
   51         mij = gretl_matrix_get(M, i, j);
   52         gretl_matrix_set(X, r, c, mij * scale);
   53     }
   54     }
   55 
   56     if (startrow != startcol) {
   57     for (i=0; i<M->rows; i++) {
   58         c = startrow + i;
   59         for (j=0; j<M->cols; j++) {
   60         r = startcol + j;
   61         mij = gretl_matrix_get(M, i, j);
   62         gretl_matrix_set(X, r, c, mij * scale);
   63         }
   64     }
   65     }
   66 }
   67 
   68 /* Retrieve the data placed on @pmod in the course of 2sls estimation:
   69    for endogenous regressors these values are the first-stage fitted
   70    values, otherwise they're just the original data.  Note that @endog
   71    is a mask with nonzero values identifying the endogenous
   72    regressors, and the special array @X contains a column for each
   73    endogenous variable.
   74 */
   75 
   76 double *model_get_Xi (const MODEL *pmod, DATASET *dset, int i)
   77 {
   78     gretl_matrix *endog = gretl_model_get_data(pmod, "endog");
   79     double *xi = NULL;
   80 
   81     if (endog == NULL || endog->val[i] == 0) {
   82     /* use the original data */
   83     xi = dset->Z[pmod->list[i+2]];
   84     } else {
   85     double **X = gretl_model_get_data(pmod, "tslsX");
   86 
   87     if (X != NULL) {
   88         /* find and return the correct column of X */
   89         int j, k = 0;
   90 
   91         for (j=0; j<i; j++) {
   92         if (endog->val[j] != 0) k++;
   93         }
   94         xi = X[k];
   95     }
   96     }
   97 
   98     return xi;
   99 }
  100 
  101 /* retrieve the special k-class data wanted for LIML estimation: these
  102    represent a further transformation of the data placed on the model
  103    via 2sls estimation.
  104 */
  105 
  106 static int make_liml_X_block (gretl_matrix *X, const MODEL *pmod,
  107                   DATASET *dset, int t1)
  108 {
  109     const double *Xi;
  110     int i, t, err = 0;
  111 
  112     X->cols = pmod->ncoeff;
  113 
  114     for (i=0; i<X->cols && !err; i++) {
  115     Xi = model_get_Xi(pmod, dset, i);
  116     if (Xi == NULL) {
  117         err = E_DATA;
  118     } else {
  119         for (t=0; t<X->rows; t++) {
  120         gretl_matrix_set(X, t, i, Xi[t+t1]);
  121         }
  122     }
  123     }
  124 
  125     return err;
  126 }
  127 
  128 /* construct the X data block pertaining to a specific equation, using
  129    either the original data or fitted values from regression on a set
  130    of instruments
  131 */
  132 
  133 static int
  134 make_sys_X_block (gretl_matrix *X, const MODEL *pmod,
  135           DATASET *dset, int t1, int method)
  136 {
  137     const double *Xi;
  138     int i, t, err = 0;
  139 
  140     X->cols = pmod->ncoeff;
  141 
  142     for (i=0; i<X->cols && !err; i++) {
  143     if (method == SYS_METHOD_3SLS ||
  144         method == SYS_METHOD_FIML ||
  145         method == SYS_METHOD_TSLS) {
  146         Xi = model_get_Xi(pmod, dset, i);
  147     } else {
  148         Xi = dset->Z[pmod->list[i+2]];
  149     }
  150     if (Xi == NULL) {
  151         err = E_DATA;
  152     } else {
  153         for (t=0; t<X->rows; t++) {
  154         gretl_matrix_set(X, t, i, Xi[t+t1]);
  155         }
  156     }
  157     }
  158 
  159     return err;
  160 }
  161 
  162 /* Populate the cross-equation covariance matrix, @S, based on
  163    the per-equation residuals: if @do_diag is non-zero we also
  164    want to compute the Breusch-Pagan test for diagonality
  165    of the matrix. For the robust variant see Halunga, Orme and
  166    Yamagata, "A heteroskasticity robust Breusch-Pagan test for
  167    contemporaneous correlation...", Journal of Econometrics,
  168    198 (2017) pp. 209-230.
  169 */
  170 
  171 static int
  172 gls_sigma_from_uhat (equation_system *sys, gretl_matrix *S,
  173              int do_diag)
  174 {
  175     int geomean = system_vcv_geomean(sys);
  176     double *rdenom = NULL;
  177     double eti, etj, sij, dij;
  178     int m = sys->neqns;
  179     int robust = 0;
  180     int i, j, t, k;
  181 
  182     if (do_diag && (sys->flags & SYSTEM_ROBUST)) {
  183     /* denominator of robust B-P test */
  184     rdenom = malloc((m * m - m) / 2 * sizeof *rdenom);
  185     if (rdenom != NULL) {
  186         robust = 1;
  187     }
  188     }
  189 
  190     k = 0;
  191     for (i=0; i<m; i++) {
  192     for (j=i; j<m; j++) {
  193         sij = dij = 0.0;
  194         for (t=0; t<sys->T; t++) {
  195         eti = gretl_matrix_get(sys->E, t, i);
  196         etj = gretl_matrix_get(sys->E, t, j);
  197         /* increment sum of cross-products */
  198         sij += eti * etj;
  199         if (i != j && robust) {
  200             /* also sum of cross-products of squares */
  201             dij += eti * etj * eti * etj;
  202         }
  203         }
  204         if (i != j && robust) {
  205         rdenom[k++] = dij;
  206         }
  207         gretl_matrix_set(S, i, j, sij);
  208         if (j != i) {
  209         gretl_matrix_set(S, j, i, sij);
  210         }
  211     }
  212     }
  213 
  214     if (do_diag) {
  215     /* compute B-P test statistic */
  216     double sii, sjj;
  217 
  218     sys->diag_test = 0.0;
  219 
  220     k = 0;
  221     for (i=0; i<m-1; i++) {
  222         sii = gretl_matrix_get(S, i, i);
  223         for (j=i+1; j<m; j++) {
  224         sij = gretl_matrix_get(S, i, j);
  225         sjj = gretl_matrix_get(S, j, j);
  226         if (robust) {
  227             /* cumulate \hat{gamma}_{ij}^2 */
  228             sys->diag_test += (sij * sij) / rdenom[k++];
  229         } else {
  230             /* cumulate \hat{rho}_{ij}^2 */
  231             sys->diag_test += (sij * sij) / (sii * sjj);
  232         }
  233         }
  234     }
  235     if (robust) {
  236         free(rdenom);
  237     } else {
  238         /* supply the missing factor of T */
  239         sys->diag_test *= sys->T;
  240     }
  241     }
  242 
  243     /* scale the S matrix */
  244     if (geomean) {
  245     for (j=0; j<S->cols; j++) {
  246         for (i=j; i<S->rows; i++) {
  247         sij = gretl_matrix_get(S, i, j);
  248         sij /= system_vcv_denom(sys, i, j);
  249         gretl_matrix_set(S, i, j, sij);
  250         if (i != j) {
  251             gretl_matrix_set(S, j, i, sij);
  252         }
  253         }
  254     }
  255     } else {
  256     gretl_matrix_divide_by_scalar(S, sys->T);
  257     }
  258 
  259     return 0;
  260 }
  261 
  262 static void maybe_correct_model_dfd (equation_system *sys,
  263                      MODEL *pmod)
  264 {
  265     int nr = system_n_restrictions(sys);
  266 
  267     if (nr > 0) {
  268     /* adjust the df correction for the number of
  269        restrictions, averaged across the equations
  270     */
  271     double avr = nr / (double) sys->neqns;
  272 
  273     pmod->dfd = sys->T - pmod->ncoeff + floor(avr);
  274     }
  275 }
  276 
  277 /* Compute residuals and related quantities, for all cases
  278    other than FIML
  279 */
  280 
  281 static void
  282 sys_resids (equation_system *sys, int eq, const DATASET *dset)
  283 {
  284     MODEL *pmod = sys->models[eq];
  285     int yno = pmod->list[1];
  286     double yh;
  287     int i, t;
  288 
  289     pmod->ess = 0.0;
  290 
  291     for (t=pmod->t1; t<=pmod->t2; t++) {
  292     yh = 0.0;
  293     for (i=0; i<pmod->ncoeff; i++) {
  294         yh += pmod->coeff[i] * dset->Z[pmod->list[i+2]][t];
  295     }
  296     pmod->yhat[t] = yh;
  297     pmod->uhat[t] = dset->Z[yno][t] - yh;
  298     /* for cross-equation vcv */
  299     gretl_matrix_set(sys->E, t - pmod->t1, pmod->ID, pmod->uhat[t]);
  300     pmod->ess += pmod->uhat[t] * pmod->uhat[t];
  301     }
  302 
  303     /* df correction? */
  304     if (system_want_df_corr(sys)) {
  305     maybe_correct_model_dfd(sys, pmod);
  306     pmod->sigma = sqrt(pmod->ess / pmod->dfd);
  307     } else {
  308     pmod->sigma = sqrt(pmod->ess / pmod->nobs);
  309     }
  310 
  311     if (pmod->ifc && pmod->tss > 0) {
  312     /* R-squared based on actual-fitted correlation */
  313     double r;
  314 
  315     pmod->rsq = gretl_corr_rsq(pmod->t1, pmod->t2, dset->Z[yno],
  316                    pmod->yhat);
  317     r = 1 - pmod->rsq;
  318     pmod->adjrsq = 1.0 - (r * (pmod->nobs - 1) / pmod->dfd);
  319     } else {
  320     pmod->rsq = pmod->adjrsq = NADBL;
  321     }
  322 }
  323 
  324 /* use the per-equation residual standard deviations for scaling
  325    the covariance matrix, block by diagonal block
  326 */
  327 
  328 static void
  329 single_eqn_scale_vcv (MODEL *pmod, gretl_matrix *V, int offset)
  330 {
  331     double vij, s2 = pmod->sigma * pmod->sigma;
  332     int i, j, ii, jj;
  333 
  334     for (i=0; i<pmod->ncoeff; i++) {
  335     ii = i + offset;
  336     for (j=i; j<pmod->ncoeff; j++) {
  337         jj = j + offset;
  338         vij = gretl_matrix_get(V, ii, jj);
  339         vij *= s2;
  340         gretl_matrix_set(V, ii, jj, vij);
  341         gretl_matrix_set(V, jj, ii, vij);
  342     }
  343     }
  344 }
  345 
  346 /* write results from system estimation into the individual
  347    model structs
  348 */
  349 
  350 static void transcribe_sys_results (equation_system *sys,
  351                     const DATASET *dset,
  352                     int do_iters)
  353 {
  354     int do_bdiff = (sys->method == SYS_METHOD_3SLS && do_iters);
  355     double bij, oldb, bnum = 0, bden = 0;
  356     int offset = 0;
  357     int i, j, k;
  358 
  359     for (i=0; i<sys->neqns; i++) {
  360     MODEL *pmod = sys->models[i];
  361 
  362     /* update the coefficients */
  363     for (j=0; j<pmod->ncoeff; j++) {
  364         k = j + offset;
  365         bij = gretl_vector_get(sys->b, k);
  366         if (do_bdiff) {
  367         oldb = pmod->coeff[j];
  368         bnum += (bij - oldb) * (bij - oldb);
  369         bden += oldb * oldb;
  370         }
  371         pmod->coeff[j] = bij;
  372     }
  373     /* update residuals (and sigma-hat) */
  374     sys_resids(sys, i, dset);
  375     /* for single-equation methods, vcv needs scaling by
  376        an estimate of the residual variance
  377     */
  378     if (sys->method == SYS_METHOD_OLS ||
  379         sys->method == SYS_METHOD_TSLS ||
  380         sys->method == SYS_METHOD_LIML) {
  381         single_eqn_scale_vcv(pmod, sys->vcv, offset);
  382     }
  383     /* update standard errors */
  384     for (j=0; j<pmod->ncoeff; j++) {
  385         k = j + offset;
  386         pmod->sderr[j] = sqrt(gretl_matrix_get(sys->vcv, k, k));
  387     }
  388     offset += pmod->ncoeff;
  389     }
  390 
  391     if (do_bdiff) {
  392     sys->bdiff = sqrt(bnum / bden);
  393     }
  394 }
  395 
  396 /* compute SUR, 3SLS or LIML parameter estimates (or restricted OLS,
  397    TSLS, WLS)
  398 */
  399 
  400 static int
  401 calculate_sys_coeffs (equation_system *sys, const DATASET *dset,
  402               gretl_matrix *X, gretl_matrix *y,
  403               int mk, int nr, int do_iters)
  404 {
  405     gretl_matrix *V;
  406     int posdef = 0;
  407     int err = 0;
  408 
  409     V = gretl_matrix_copy(X);
  410     if (V == NULL) {
  411     return E_ALLOC;
  412     }
  413 
  414 #if SDEBUG
  415     gretl_matrix_print(X, "sys X");
  416     gretl_matrix_print(y, "sys y");
  417 #endif
  418 
  419     /* calculate the coefficients */
  420 
  421     if (sys->R == NULL) {
  422     /* hopefully X may be positive definite */
  423     err = gretl_cholesky_decomp_solve(X, y);
  424     if (!err) {
  425         posdef = 1;
  426         err = gretl_cholesky_invert(X);
  427     } else {
  428         /* nope: fall back to the LU solver */
  429         gretl_matrix_copy_values(X, V);
  430         err = 0;
  431     }
  432     }
  433 
  434     if (!posdef) {
  435     err = gretl_LU_solve_invert(X, y);
  436     }
  437 
  438     if (!err) {
  439     /* save the coefficient vector and covariance matrix */
  440     gretl_matrix *b = gretl_matrix_copy(y);
  441 
  442     if (b == NULL) {
  443         err = E_ALLOC;
  444     } else {
  445         system_attach_coeffs(sys, b);
  446         gretl_matrix_copy_values(V, X);
  447         system_attach_vcv(sys, V);
  448         /* transcribe stuff to the included models */
  449         transcribe_sys_results(sys, dset, do_iters);
  450     }
  451     }
  452 
  453     if (err) {
  454     gretl_matrix_free(V);
  455     }
  456 
  457     return err;
  458 }
  459 
  460 static int *
  461 system_model_list (equation_system *sys, int i, int *freeit)
  462 {
  463     int *list = NULL;
  464 
  465     *freeit = 0;
  466 
  467     if (sys_ols_ok(sys)) {
  468     return system_get_list(sys, i);
  469     }
  470 
  471     if (sys->method != SYS_METHOD_FIML) {
  472     list = system_get_list(sys, i);
  473     }
  474 
  475     if (sys->method == SYS_METHOD_3SLS ||
  476     sys->method == SYS_METHOD_TSLS ||
  477     sys->method == SYS_METHOD_LIML) {
  478     /* Is the list already in ivreg form?
  479        If not, ignore it. */
  480     if (list != NULL && !gretl_list_has_separator(list)) {
  481         list = NULL;
  482     }
  483     }
  484 
  485     if ((sys->method == SYS_METHOD_FIML ||
  486      sys->method == SYS_METHOD_LIML ||
  487      sys->method == SYS_METHOD_3SLS ||
  488      sys->method == SYS_METHOD_TSLS)
  489     && list == NULL) {
  490     list = compose_ivreg_list(sys, i);
  491     *freeit = 1;
  492     }
  493 
  494     return list;
  495 }
  496 
  497 /* Hansen-Sargan overidentification test for the system as a whole, as
  498    in Davidson and MacKinnon, ETM: p. 511 and equation (12.25) for the
  499    case of SUR; p. 532 and equation (12.61) for the case of 3SLS.  See
  500    also D & M, Estimation and Inference in Econometrics, equation
  501    (18.60), for a more computation-friendly statement of the criterion
  502    function.
  503 */
  504 
  505 static int hansen_sargan_test (equation_system *sys,
  506                    const DATASET *dset)
  507 {
  508     const int *exlist = system_get_instr_vars(sys);
  509     const double *Wi, *Wj;
  510     int nx = exlist[0];
  511     int m = sys->neqns;
  512     int T = sys->T;
  513     int df = system_get_overid_df(sys);
  514     gretl_matrix_block *B;
  515     gretl_matrix *WTW, *eW, *tmp;
  516     double x, X2;
  517     int i, j, t;
  518     int err = 0;
  519 
  520     if (df <= 0) {
  521     return 1;
  522     }
  523 
  524     B = gretl_matrix_block_new(&WTW, nx, nx,
  525                    &eW, m, nx,
  526                    &tmp, m, nx,
  527                    NULL);
  528     if (B == NULL) {
  529     return E_ALLOC;
  530     }
  531 
  532     /* construct W-transpose W in WTW */
  533     for (i=0; i<nx; i++) {
  534     Wi = dset->Z[exlist[i+1]] + sys->t1;
  535     for (j=i; j<nx; j++) {
  536         Wj = dset->Z[exlist[j+1]] + sys->t1;
  537         x = 0.0;
  538         for (t=0; t<sys->T; t++) {
  539         x += Wi[t] * Wj[t];
  540         }
  541         gretl_matrix_set(WTW, i, j, x);
  542         if (i != j) {
  543         gretl_matrix_set(WTW, j, i, x);
  544         }
  545     }
  546     }
  547 
  548     err = gretl_invert_symmetric_matrix(WTW);
  549 #if SDEBUG
  550     fprintf(stderr, "hansen_sargan: on invert, err=%d\n", err);
  551 #endif
  552     if (err) {
  553     sys->X2 = NADBL;
  554     goto bailout;
  555     }
  556 
  557     /* set up vectors of SUR or 3SLS residuals, transposed,
  558        times W; these are stacked in the m * nx matrix eW
  559     */
  560     for (i=0; i<m; i++) {
  561     for (j=0; j<nx; j++) {
  562         Wj = dset->Z[exlist[j+1]] + sys->t1;
  563         x = 0.0;
  564         for (t=0; t<T; t++) {
  565         x += gretl_matrix_get(sys->E, t, i) * Wj[t];
  566         }
  567         gretl_matrix_set(eW, i, j, x);
  568     }
  569     }
  570 
  571     /* multiply these vectors into (WTW)^{-1} */
  572     for (i=0; i<m; i++) {
  573     for (j=0; j<nx; j++) {
  574         x = 0.0;
  575         for (t=0; t<nx; t++) {
  576         x += gretl_matrix_get(eW, i, t) *
  577             gretl_matrix_get(WTW, t, j);
  578         }
  579         gretl_matrix_set(tmp, i, j, x);
  580     }
  581     }
  582 
  583     /* cumulate the Chi-square value */
  584     X2 = 0.0;
  585     for (i=0; i<m; i++) {
  586     for (j=0; j<m; j++) {
  587         x = 0.0;
  588         for (t=0; t<nx; t++) {
  589         x += gretl_matrix_get(tmp, i, t) *
  590             gretl_matrix_get(eW, j, t); /* transposed */
  591         }
  592         X2 += gretl_matrix_get(sys->S, i, j) * x;
  593     }
  594     }
  595 
  596 #if SDEBUG
  597     fprintf(stderr, "Hansen-Sargan: Chi-square(%d) = %g (p-value %g)\n",
  598         df, X2, chisq_cdf_comp(df, X2));
  599 #endif
  600     sys->X2 = X2;
  601 
  602  bailout:
  603 
  604     gretl_matrix_block_destroy(B);
  605 
  606     return err;
  607 }
  608 
  609 static int basic_system_allocate (equation_system *sys,
  610                   int mk, int nr,
  611                   gretl_matrix **X,
  612                   gretl_matrix **y)
  613 {
  614     int m = sys->neqns;
  615     int T = sys->T;
  616     int ldx = mk + nr;
  617 
  618     /* allocate a model for each stochastic equation */
  619     sys->models = gretl_model_array_new(m);
  620     if (sys->models == NULL) {
  621     return E_ALLOC;
  622     }
  623 
  624     sys->E = gretl_matrix_alloc(T, m);
  625     if (sys->E == NULL) {
  626     return E_ALLOC;
  627     }
  628 
  629     sys->S = gretl_matrix_alloc(m, m);
  630     if (sys->S == NULL) {
  631     return E_ALLOC;
  632     }
  633 
  634     if (X != NULL) {
  635     *X = gretl_matrix_alloc(ldx, ldx);
  636     if (*X == NULL) {
  637         return E_ALLOC;
  638     }
  639     }
  640 
  641     if (y != NULL) {
  642     *y = gretl_column_vector_alloc(ldx);
  643     if (*y == NULL) {
  644         return E_ALLOC;
  645     }
  646     }
  647 
  648     return 0;
  649 }
  650 
  651 /* compute log-likelihood for iterated SUR estimator */
  652 
  653 double sur_loglik (equation_system *sys)
  654 {
  655     int m = sys->neqns;
  656     int T = sys->T;
  657     gretl_matrix *tmp;
  658     double ldet;
  659     int err = 0;
  660 
  661     tmp = gretl_matrix_alloc(m, m);
  662     if (tmp == NULL) {
  663     return NADBL;
  664     }
  665 
  666     gls_sigma_from_uhat(sys, tmp, 0);
  667     ldet = gretl_vcv_log_determinant(tmp, &err);
  668 
  669     if (na(ldet)) {
  670     sys->ll = NADBL;
  671     } else {
  672     sys->ll = -(m * T / 2.0) * (LN_2_PI + 1.0);
  673     sys->ll -= (T / 2.0) * ldet;
  674     }
  675 
  676     gretl_matrix_free(tmp);
  677 
  678     return sys->ll;
  679 }
  680 
  681 /* If we're estimating with a specified set of linear restrictions,
  682    Rb = q, augment the X matrix with R and R-transpose.
  683 */
  684 
  685 static void
  686 augment_X_with_restrictions (gretl_matrix *X, int mk,
  687                  equation_system *sys)
  688 {
  689     int i, j, nr = sys->R->rows;
  690 
  691     /* place the R matrix */
  692     insert_sys_X_block(X, sys->R, mk, 0, 1.0);
  693 
  694     /* zero the bottom right-hand block */
  695     for (i=mk; i<mk+nr; i++) {
  696     for (j=mk; j<mk+nr; j++) {
  697         gretl_matrix_set(X, i, j, 0.0);
  698     }
  699     }
  700 }
  701 
  702 /* When estimating with a specified set of linear restrictions,
  703    Rb = q, augment the y matrix with q.
  704 */
  705 
  706 static int
  707 augment_y_with_restrictions (gretl_matrix *y, int mk, int nr,
  708                  equation_system *sys)
  709 {
  710     int i;
  711 
  712     if (sys->q == NULL) {
  713     return 1;
  714     }
  715 
  716     for (i=0; i<nr; i++) {
  717     gretl_vector_set(y, mk + i, gretl_vector_get(sys->q, i));
  718     }
  719 
  720     return 0;
  721 }
  722 
  723 #define SYS_MAX_ITER 100
  724 #define SYS_LL_TOL 1.0e-12
  725 #define SYS_BDIFF_TOL 1.0e-9
  726 
  727 /* check for convergence of iteration: we use the change in the
  728    log-likelihood when iterating SUR to the ML solution, or
  729    a measure of the change in the coefficients when iterating
  730    3SLS
  731 */
  732 
  733 static int sys_converged (equation_system *sys, double *llbak,
  734               gretlopt opt, PRN *prn, int *err)
  735 {
  736     double crit, tol = 0.0;
  737     int met = 0;
  738 
  739     if (sys->method == SYS_METHOD_SUR ||
  740     sys->method == SYS_METHOD_WLS) {
  741     double ll = sur_loglik(sys);
  742 
  743     tol = SYS_LL_TOL;
  744     crit = ll - *llbak;
  745 
  746     if (opt & OPT_V) {
  747         pprintf(prn, "Iteration %3d, ll = %#.8g\n", sys->iters, ll);
  748     }
  749     if (crit <= tol) {
  750         met = 1;
  751     } else if (sys->iters < SYS_MAX_ITER) {
  752         *llbak = ll;
  753     }
  754     } else if (sys->method == SYS_METHOD_3SLS) {
  755     tol = SYS_BDIFF_TOL;
  756     crit = sys->bdiff;
  757     if (opt & OPT_V) {
  758         pprintf(prn, "Iteration %3d, criterion = %g\n", sys->iters, crit);
  759     }
  760     if (crit <= tol) {
  761         met = 1;
  762     }
  763     }
  764 
  765     if (met && tol > 0 && (opt & OPT_V)) {
  766     pprintf(prn, "Tolerance of %g is met\n", tol);
  767     }
  768 
  769     if (!met && sys->iters >= SYS_MAX_ITER) {
  770     pprintf(prn, "Reached %d iterations without meeting "
  771         "tolerance of %g\n", sys->iters, tol);
  772     *err = E_NOCONV;
  773     }
  774 
  775     return met;
  776 }
  777 
  778 static void clean_up_models (equation_system *sys)
  779 {
  780     int i;
  781 
  782     if (sys->models == NULL) {
  783     /* "can't happen" */
  784     return;
  785     }
  786 
  787     sys->ess = 0.0;
  788 
  789     for (i=0; i<sys->neqns; i++) {
  790     sys->ess += sys->models[i]->ess;
  791     if (sys->method == SYS_METHOD_3SLS ||
  792         sys->method == SYS_METHOD_FIML ||
  793         sys->method == SYS_METHOD_TSLS ||
  794         sys->method == SYS_METHOD_LIML) {
  795         tsls_free_data(sys->models[i]);
  796     }
  797     if (sys->neqns > 1) {
  798         gretl_model_free(sys->models[i]);
  799     }
  800     }
  801 
  802     if (sys->flags & SYSTEM_LIML1) {
  803     /* ivreg (single-equation LIML) */
  804     sys->models[0]->rho = sys->models[0]->dw = NADBL;
  805     } else {
  806     free(sys->models);
  807     sys->models = NULL;
  808     }
  809 }
  810 
  811 /* Given a record of the redundant instruments dropped at the tsls
  812    stage, namely @idroplist, purge these series IDs from both the
  813    global system instrument list (sys->ilist) and the per-equation
  814    regression specification list (if the latter is given in tsls
  815    form).
  816 */
  817 
  818 static int drop_redundant_instruments (equation_system *sys,
  819                        const int *idroplist,
  820                        int eqn)
  821 {
  822     int i, j, di, pos;
  823     int err = 0;
  824 
  825     for (i=1; i<=idroplist[0]; i++) {
  826     di = idroplist[i];
  827     pos = in_gretl_list(sys->ilist, di);
  828     if (pos > 0) {
  829         gretl_list_delete_at_pos(sys->ilist, pos);
  830     } else {
  831         err = 1;
  832     }
  833     }
  834 
  835     pos = gretl_list_separator_position(sys->lists[eqn]);
  836 
  837     if (pos > 0) {
  838     int *eqnlist = sys->lists[eqn];
  839 
  840     for (i=1; i<=idroplist[0]; i++) {
  841         di = idroplist[i];
  842         for (j=pos+1; j<=eqnlist[0]; j++) {
  843         if (eqnlist[j] == di) {
  844             gretl_list_delete_at_pos(eqnlist, j);
  845             break;
  846         }
  847         }
  848     }
  849     }
  850 
  851     return err;
  852 }
  853 
  854 /* options to be passed in running initial 2SLS */
  855 
  856 static gretlopt sys_tsls_opt (const equation_system *sys)
  857 {
  858     gretlopt opt = (OPT_E | OPT_A);
  859 
  860     if (!(sys->flags & SYSTEM_DFCORR)) {
  861     opt |= OPT_N; /* suppress df correction */
  862     }
  863 
  864     if (sys->flags & SYSTEM_LIML1) {
  865     opt |= OPT_H; /* add "hatlist" of instrumented vars */
  866     }
  867 
  868     return opt;
  869 }
  870 
  871 /* options to be passed in running initial OLS */
  872 
  873 static gretlopt sys_ols_opt (const equation_system *sys,
  874                  int nr)
  875 {
  876     gretlopt opt = OPT_NONE;
  877 
  878     if (sys->method == SYS_METHOD_OLS) {
  879     if (!(sys->flags & SYSTEM_DFCORR)) {
  880         opt |= OPT_N; /* suppress df correction */
  881     }
  882     }
  883 
  884     if (sys->method == SYS_METHOD_OLS && nr == 0) {
  885     /* initial OLS will supply the estimates */
  886     if (sys->flags & SYSTEM_ROBUST) {
  887         opt |= OPT_R;
  888     }
  889     } else {
  890     /* treat initial OLS as auxiliary */
  891     opt |= OPT_A;
  892     }
  893 
  894     return opt;
  895 }
  896 
  897 static int allocate_Xi_etc (gretl_matrix **Xi,
  898                 gretl_matrix **Xj,
  899                 gretl_matrix **M,
  900                 int T, int k)
  901 {
  902     *Xi = gretl_matrix_alloc(T, k);
  903     *Xj = gretl_matrix_alloc(T, k);
  904     *M = gretl_matrix_alloc(k, k);
  905 
  906     if (*Xi == NULL || *Xj == NULL || *M == NULL) {
  907     return E_ALLOC;
  908     } else {
  909     return 0;
  910     }
  911 }
  912 
  913 static int ols_data_to_sys (equation_system *sys, int mk)
  914 {
  915     gretl_matrix *B, *V;
  916     MODEL *pmod;
  917     double vij;
  918     int i, j, k;
  919     int mi, mj;
  920     int vi, vj;
  921     int err = 0;
  922 
  923     B = gretl_matrix_alloc(mk, 1);
  924     V = gretl_zero_matrix_new(mk, mk);
  925     if (B == NULL || V == NULL) {
  926     return E_ALLOC;
  927     }
  928 
  929     j = vi = vj = 0;
  930     for (i=0; i<sys->neqns; i++) {
  931     pmod = sys->models[i];
  932     if (pmod->vcv == NULL) {
  933         err = makevcv(pmod, pmod->sigma);
  934         if (err) {
  935         break;
  936         }
  937     }
  938     k = pmod->ncoeff;
  939     for (mi=0; mi<k; mi++) {
  940         B->val[j++] = pmod->coeff[mi];
  941         for (mj=0; mj<k; mj++) {
  942         vij = gretl_model_get_vcv_element(pmod, mi, mj, k);
  943         gretl_matrix_set(V, vi+mi, vj+mj, vij);
  944         }
  945     }
  946     vi += k;
  947     vj += k;
  948     }
  949 
  950     if (err) {
  951     gretl_matrix_free(B);
  952     gretl_matrix_free(V);
  953     } else {
  954     system_attach_coeffs(sys, B);
  955     system_attach_vcv(sys, V);
  956     }
  957 
  958     return err;
  959 }
  960 
  961 /* general function that forms the basis for all specific system
  962    estimators */
  963 
  964 int system_estimate (equation_system *sys, DATASET *dset,
  965              gretlopt opt, PRN *prn)
  966 {
  967     int i, j, k, T, t;
  968     int v, l, mk, krow, nr;
  969     int orig_t1 = dset->t1;
  970     int orig_t2 = dset->t2;
  971     gretl_matrix *X = NULL;
  972     gretl_matrix *y = NULL;
  973     gretl_matrix *Xi = NULL;
  974     gretl_matrix *Xj = NULL;
  975     gretl_matrix *M = NULL;
  976     gretl_matrix **pX = NULL;
  977     gretl_matrix **py = NULL;
  978     MODEL **models = NULL;
  979     int method = sys->method;
  980     double llbak = -1.0e9;
  981     int single_equation = 0;
  982     int do_iteration = 0;
  983     int plain_ols = 0;
  984     int rsingle = 0;
  985     int do_diag = 0;
  986     int err = 0;
  987 
  988     sys->iters = 0;
  989 
  990     if (sys->flags & SYSTEM_ITERATE) {
  991     do_iteration = 1;
  992     }
  993 
  994     nr = system_n_restrictions(sys);
  995 
  996     if (method == SYS_METHOD_OLS || method == SYS_METHOD_TSLS ||
  997     method == SYS_METHOD_LIML || method == SYS_METHOD_WLS) {
  998     single_equation = 1;
  999     }
 1000 
 1001     if (method == SYS_METHOD_OLS && nr == 0) {
 1002     plain_ols = 1;
 1003     } else {
 1004     pX = &X;
 1005     py = &y;
 1006     }
 1007 
 1008     if (nr > 0 && !(opt & OPT_U)) {
 1009     if (method == SYS_METHOD_3SLS) {
 1010         /* doing 3SLS with restrictions: we want to obtain
 1011            restricted TSLS estimates as a starting point
 1012         */
 1013         rsingle = 1;
 1014     } else if (method == SYS_METHOD_SUR) {
 1015         /* doing SUR with restrictions: we want to obtain
 1016            restricted OLS estimates as a starting point
 1017         */
 1018         rsingle = 1;
 1019     }
 1020     }
 1021 
 1022     /* get uniform sample starting and ending points and check for
 1023        missing data */
 1024     err = system_adjust_t1t2(sys, dset);
 1025     if (err) {
 1026     return err;
 1027     }
 1028 
 1029     /* set sample for auxiliary regressions */
 1030     dset->t1 = sys->t1;
 1031     dset->t2 = sys->t2;
 1032 
 1033     /* max indep vars per equation */
 1034     k = system_max_indep_vars(sys);
 1035 
 1036     /* total indep vars, all equations */
 1037     mk = system_n_indep_vars(sys);
 1038 
 1039     /* set sample for auxiliary regressions */
 1040     dset->t1 = sys->t1;
 1041     dset->t2 = sys->t2;
 1042 
 1043     /* number of observations per series */
 1044     T = sys->T;
 1045 
 1046     /* allocate models etc */
 1047     err = basic_system_allocate(sys, mk, nr, pX, py);
 1048     if (err) {
 1049     goto cleanup;
 1050     }
 1051 
 1052     /* convenience pointers */
 1053     models = sys->models;
 1054 
 1055     if ((method == SYS_METHOD_FIML ||
 1056      method == SYS_METHOD_LIML) && !(opt & OPT_Q)) {
 1057     print_equation_system_info(sys, dset, OPT_H, prn);
 1058     }
 1059 
 1060     /* First estimate the equations separately (either by OLS or
 1061        TSLS), and put the single-equation residuals into the uhat
 1062        matrix.  Note that at this stage we are not in a position to
 1063        impose any cross-equation restrictions, since we're doing
 1064        straight equation-by-equation estimation.
 1065     */
 1066 
 1067     for (i=0; i<sys->neqns; i++) {
 1068     int freeit = 0;
 1069     int *list = system_model_list(sys, i, &freeit);
 1070     const int *idroplist = NULL;
 1071     MODEL *pmod = models[i];
 1072 
 1073     if (list == NULL) {
 1074         err = 1;
 1075         break;
 1076     }
 1077 
 1078     if (sys_ols_ok(sys)) {
 1079         *pmod = lsq(list, dset, OLS, sys_ols_opt(sys, nr));
 1080     } else {
 1081         *pmod = tsls(list, dset, sys_tsls_opt(sys));
 1082     }
 1083 
 1084     if (freeit) {
 1085         free(list);
 1086     }
 1087 
 1088     if ((err = pmod->errcode)) {
 1089         fprintf(stderr, "system_estimate: failed to estimate equation %d: "
 1090             "err = %d\n", i+1, err);
 1091         break;
 1092     }
 1093 
 1094     idroplist = gretl_model_get_list(pmod, "inst_droplist");
 1095     if (idroplist != NULL) {
 1096         drop_redundant_instruments(sys, idroplist, i);
 1097     }
 1098 
 1099     pmod->ID = i;
 1100     pmod->aux = AUX_SYS;
 1101     gretl_model_set_int(pmod, "method", method);
 1102 
 1103     /* save the sigma-squared for an LR test for a diagonal
 1104        covariance matrix */
 1105     if (method == SYS_METHOD_SUR && do_iteration && nr == 0) {
 1106         gretl_model_set_double(pmod, "ols_sigma_squared",
 1107                    pmod->ess / pmod->nobs);
 1108     }
 1109 
 1110     for (t=0; t<T; t++) {
 1111         gretl_matrix_set(sys->E, t, i, pmod->uhat[t + sys->t1]);
 1112     }
 1113     }
 1114 
 1115     if (err) {
 1116     fprintf(stderr, "system_estimate: after single-equation "
 1117         "estimation, err = %d\n", err);
 1118     goto cleanup;
 1119     }
 1120 
 1121     if (method == SYS_METHOD_LIML) {
 1122     /* compute the minimum eigenvalues and generate the
 1123        k-class data matrices */
 1124     err = liml_driver(sys, dset, prn);
 1125     if (err) goto cleanup;
 1126     }
 1127 
 1128     if (plain_ols) {
 1129     gls_sigma_from_uhat(sys, sys->S, 1);
 1130     err = ols_data_to_sys(sys, mk);
 1131     goto save_etc;
 1132     }
 1133 
 1134     /* marker for iterated versions of SUR, WLS, or 3SLS; also for
 1135        loopback in case of restricted 3SLS, where we want to compute
 1136        restricted TSLS estimates first
 1137     */
 1138  iteration_start:
 1139 
 1140     gls_sigma_from_uhat(sys, sys->S, 0);
 1141 
 1142     if (method == SYS_METHOD_WLS) {
 1143     gretl_matrix_zero(X);
 1144     err = gretl_invert_diagonal_matrix(sys->S);
 1145     } else if (single_equation || rsingle) {
 1146     gretl_matrix_zero(X);
 1147     } else {
 1148     err = gretl_invert_symmetric_matrix(sys->S);
 1149     }
 1150 
 1151 #if SDEBUG
 1152     fprintf(stderr, "system_estimate: on invert, err=%d\n", err);
 1153 #endif
 1154 
 1155     if (!err && Xi == NULL) {
 1156     /* the test against NULL here allows for the possibility
 1157        that we're iterating
 1158     */
 1159     err = allocate_Xi_etc(&Xi, &Xj, &M, T, k);
 1160     }
 1161 
 1162     if (err) goto cleanup;
 1163 
 1164     /* form the big stacked X matrix: Xi = data matrix for equation i,
 1165        specified in lists[i]
 1166     */
 1167 
 1168     krow = 0;
 1169     for (i=0; i<sys->neqns && !err; i++) {
 1170     int kcol = 0;
 1171 
 1172     err = make_sys_X_block(Xi, models[i], dset, sys->t1, method);
 1173 
 1174     for (j=0; j<=i && !err; j++) {
 1175         const gretl_matrix *Xk;
 1176         double sij;
 1177 
 1178         if (i != j) {
 1179         if (single_equation || rsingle) {
 1180             kcol += models[j]->ncoeff;
 1181             continue;
 1182         }
 1183         err = make_sys_X_block(Xj, models[j], dset, sys->t1, method);
 1184         Xk = Xj;
 1185         } else if (method == SYS_METHOD_LIML) {
 1186         err = make_liml_X_block(Xj, models[i], dset, sys->t1);
 1187         Xk = Xj;
 1188         } else {
 1189         Xk = Xi;
 1190         }
 1191 
 1192         M->rows = Xi->cols;
 1193         M->cols = Xk->cols;
 1194 
 1195         err = gretl_matrix_multiply_mod(Xi, GRETL_MOD_TRANSPOSE,
 1196                         Xk, GRETL_MOD_NONE,
 1197                         M, GRETL_MOD_NONE);
 1198 
 1199         if (rsingle || (single_equation && method != SYS_METHOD_WLS)) {
 1200         sij = 1.0;
 1201         } else {
 1202         sij = gretl_matrix_get(sys->S, i, j);
 1203         }
 1204 
 1205         insert_sys_X_block(X, M, krow, kcol, sij);
 1206         kcol += models[j]->ncoeff;
 1207     }
 1208 
 1209     krow += models[i]->ncoeff;
 1210     }
 1211 
 1212     if (err) {
 1213     fprintf(stderr, "after trying to make X matrix: err = %d\n", err);
 1214     goto cleanup;
 1215     }
 1216 
 1217     if (nr > 0) {
 1218     /* there are restrictions to be imposed */
 1219     augment_X_with_restrictions(X, mk, sys);
 1220     }
 1221 
 1222     if (!do_iteration && !rsingle) {
 1223     /* we're not coming back this way, so free some storage */
 1224     gretl_matrix_free(Xj);
 1225     Xj = NULL;
 1226     gretl_matrix_free(M);
 1227     M = NULL;
 1228     }
 1229 
 1230     /* form stacked Y column vector (m x k) */
 1231 
 1232     v = 0;
 1233     for (i=0; i<sys->neqns; i++) {
 1234 
 1235     /* loop over the m vertically-arranged blocks */
 1236     make_sys_X_block(Xi, models[i], dset, sys->t1, method);
 1237 
 1238     for (j=0; j<models[i]->ncoeff; j++) {
 1239         /* loop over the rows within each of the m blocks */
 1240         double yv = 0.0;
 1241         int lmin = 0, lmax = sys->neqns;
 1242 
 1243         if (single_equation || rsingle) {
 1244         /* no cross terms wanted */
 1245         lmin = i;
 1246         lmax = i + 1;
 1247         }
 1248 
 1249         for (l=lmin; l<lmax; l++) {
 1250         /* loop over the components that must be
 1251            added to form each element */
 1252         const double *yl = NULL;
 1253         double sil, xx = 0.0;
 1254 
 1255         if (method == SYS_METHOD_LIML) {
 1256             yl = gretl_model_get_data(models[l], "liml_y");
 1257         } else {
 1258             yl = dset->Z[system_get_depvar(sys, l)];
 1259         }
 1260 
 1261         /* multiply X'[l] into y */
 1262         for (t=0; t<T; t++) {
 1263             xx += gretl_matrix_get(Xi, t, j) * yl[t + sys->t1];
 1264         }
 1265 
 1266         if (rsingle || (single_equation && method != SYS_METHOD_WLS)) {
 1267             sil = 1.0;
 1268         } else {
 1269             sil = gretl_matrix_get(sys->S, i, l);
 1270         }
 1271 
 1272         yv += xx * sil;
 1273         }
 1274         gretl_vector_set(y, v++, yv);
 1275     }
 1276     }
 1277 
 1278     if (nr > 0) {
 1279     /* there are restrictions */
 1280     augment_y_with_restrictions(y, mk, nr, sys);
 1281     }
 1282 
 1283     /* The estimates calculated below will be SUR, 3SLS or LIML,
 1284        depending on how the data matrices above were constructed --
 1285        unless, that is, we're just doing restricted OLS, WLS or TSLS
 1286        estimates.
 1287     */
 1288     err = calculate_sys_coeffs(sys, dset, X, y, mk, nr,
 1289                    do_iteration);
 1290 
 1291     if (!err && rsingle) {
 1292     /* take one more pass */
 1293     rsingle = 0;
 1294     goto iteration_start;
 1295     }
 1296 
 1297     if (!err && do_iteration) {
 1298     /* check for convergence */
 1299     if (!sys_converged(sys, &llbak, opt, prn, &err)) {
 1300         if (!err) {
 1301         sys->iters += 1;
 1302         goto iteration_start;
 1303         }
 1304     }
 1305     }
 1306 
 1307     if (err) goto cleanup;
 1308 
 1309     if (nr == 0 && (method == SYS_METHOD_3SLS || method == SYS_METHOD_SUR)) {
 1310     /* compute this test while we have sigma-inverse available */
 1311     hansen_sargan_test(sys, dset);
 1312     }
 1313 
 1314     do_diag = !(method == SYS_METHOD_SUR && do_iteration);
 1315 
 1316     /* refresh sigma (non-inverted) */
 1317     gls_sigma_from_uhat(sys, sys->S, do_diag);
 1318 
 1319     if (method == SYS_METHOD_FIML) {
 1320     /* compute FIML estimates */
 1321     err = fiml_driver(sys, dset, opt, prn);
 1322     }
 1323 
 1324  save_etc:
 1325 
 1326     if (!err && !(sys->flags & SYSTEM_LIML1)) {
 1327     err = system_save_and_print_results(sys, dset, opt, prn);
 1328     }
 1329 
 1330  cleanup:
 1331 
 1332     gretl_matrix_free(Xi);
 1333     gretl_matrix_free(Xj);
 1334     gretl_matrix_free(M);
 1335     gretl_matrix_free(X);
 1336     gretl_matrix_free(y);
 1337 
 1338     clean_up_models(sys);
 1339 
 1340     dset->t1 = orig_t1;
 1341     dset->t2 = orig_t2;
 1342 
 1343     return err;
 1344 }