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 /* mechanisms for defining and handling systems of equations */
21
22 #define FULL_XML_HEADERS
23
24 #include "libgretl.h"
25 #include "system.h"
26 #include "objstack.h"
27 #include "usermat.h"
28 #include "gretl_xml.h"
29 #include "gretl_func.h"
30 #include "tsls.h"
31 #include "uservar.h"
32
33 #define SYSDEBUG 0
34
35 enum {
36 OP_PLUS,
37 OP_MINUS
38 } identity_ops;
39
40 enum {
41 ENDOG_LIST,
42 INSTR_LIST
43 } aux_list_types;
44
45 enum {
46 SYS_TEST_NONE,
47 SYS_TEST_LR,
48 SYS_TEST_F,
49 SYS_TEST_NOTIMP
50 } system_test_types;
51
52 enum {
53 SYS_UHAT,
54 SYS_YHAT
55 };
56
57 struct id_atom_ {
58 int op; /* operator (plus or minus) */
59 int varnum; /* ID number of variable to right of operator */
60 };
61
62 struct identity_ {
63 int n_atoms; /* number of "atomic" elements in identity */
64 int depvar; /* LHS variable in indentity */
65 id_atom *atoms; /* pointer to RHS "atoms" */
66 };
67
68 struct predet_ {
69 int id; /* ID number of lag variable */
70 int src; /* ID number of "parent" */
71 int lag; /* lag order */
72 };
73
74 struct liml_data_ {
75 double *lmin; /* min. eigenvalues */
76 double *ll; /* per-equation log likelihood */
77 int *idf; /* overidentification df */
78 };
79
80 const char *gretl_system_method_strings[] = {
81 "sur",
82 "3sls",
83 "fiml",
84 "liml",
85 "ols",
86 "tsls",
87 "wls",
88 NULL
89 };
90
91 const char *gretl_system_short_strings[] = {
92 N_("SUR"),
93 N_("3SLS"),
94 N_("FIML"),
95 N_("LIML"),
96 N_("OLS"),
97 N_("TSLS"),
98 N_("WLS"),
99 NULL
100 };
101
102 const char *gretl_system_long_strings[] = {
103 N_("Seemingly Unrelated Regressions"),
104 N_("Three-Stage Least Squares"),
105 N_("Full Information Maximum Likelihood"),
106 N_("Limited Information Maximum Likelihood"),
107 N_("Ordinary Least Squares"),
108 N_("Two-Stage Least Squares"),
109 N_("Weighted Least Squares"),
110 NULL
111 };
112
113 const char *nosystem = N_("No system of equations has been defined");
114 const char *badsystem = N_("Unrecognized equation system type");
115 const char *toofew = N_("An equation system must have at least two equations");
116
117 static void destroy_ident (identity *ident);
118 static int
119 add_predet_to_sys (equation_system *sys, const DATASET *dset,
120 int id, int src, int lag);
121 static int sys_check_lists (equation_system *sys, const DATASET *dset);
122
123 #define sys_anonymous(s) (strcmp(s, "$system") == 0)
124
125 static GList *sysstack;
126
get_anon_system_at_depth(int fd)127 static equation_system *get_anon_system_at_depth (int fd)
128 {
129 equation_system *sys;
130 GList *tmp = sysstack;
131
132 while (tmp != NULL) {
133 sys = tmp->data;
134 if (sys->fd == fd) {
135 return sys;
136 }
137 tmp = tmp->next;
138 }
139
140 return NULL;
141 }
142
get_anonymous_equation_system(void)143 equation_system *get_anonymous_equation_system (void)
144 {
145 int fd = gretl_function_depth();
146
147 return get_anon_system_at_depth(fd);
148 }
149
150 /* The use-case for push_anon_system() is when a user wants to define
151 a system and then estimate it by more than one method (or set it as
152 a target for "restrict" before estimation) but does not care to
153 give it a specific name via the "name <- system" mechanism -- or
154 we're inside a user-defined function where this naming mechanism is
155 not available.
156
157 In this case we save the system "anonymously" instead (under a name
158 of "$system").
159 */
160
push_anon_system(equation_system * sys)161 static void push_anon_system (equation_system *sys)
162 {
163 equation_system *old = get_anon_system_at_depth(sys->fd);
164
165 if (old != NULL) {
166 sysstack = g_list_remove(sysstack, old);
167 gretl_object_unref(old, GRETL_OBJ_SYS);
168 }
169
170 gretl_object_ref(sys, GRETL_OBJ_SYS);
171 sysstack = g_list_append(sysstack, sys);
172 }
173
174 /* For use when terminating function execution: if an
175 anonymous equation system was set up at the given
176 level of function execution, trash it.
177 */
178
delete_anonymous_equation_system(int level)179 void delete_anonymous_equation_system (int level)
180 {
181 equation_system *sys = get_anon_system_at_depth(level);
182
183 if (sys != NULL) {
184 sysstack = g_list_remove(sysstack, sys);
185 gretl_object_unref(sys, GRETL_OBJ_SYS);
186 }
187 }
188
189 static void
print_system_equation(const int * list,const DATASET * dset,PRN * prn)190 print_system_equation (const int *list, const DATASET *dset,
191 PRN *prn)
192 {
193 int i, v;
194
195 pputs(prn, "equation");
196
197 for (i=1; i<=list[0]; i++) {
198 v = list[i];
199 if (v == LISTSEP) {
200 pputs(prn, " ;");
201 } else if (v > 0 && v < dset->v) {
202 pprintf(prn, " %s", dset->varname[v]);
203 } else {
204 pprintf(prn, " %d", v);
205 }
206 }
207
208 pputc(prn, '\n');
209 }
210
211 static void
print_system_identity(const identity * ident,const DATASET * dset,gretlopt opt,PRN * prn)212 print_system_identity (const identity *ident, const DATASET *dset,
213 gretlopt opt, PRN *prn)
214 {
215 int i;
216
217 if (opt & OPT_H) {
218 pprintf(prn, "Identity: %s = %s ",
219 dset->varname[ident->depvar],
220 dset->varname[ident->atoms[0].varnum]);
221 } else {
222 pprintf(prn, "identity %s = %s ",
223 dset->varname[ident->depvar],
224 dset->varname[ident->atoms[0].varnum]);
225 }
226
227 for (i=1; i<ident->n_atoms; i++) {
228 pprintf(prn, "%c %s ", (ident->atoms[i].op == OP_PLUS)? '+' : '-',
229 dset->varname[ident->atoms[i].varnum]);
230 }
231
232 pputc(prn, '\n');
233 }
234
sys_max_predet_lag(const equation_system * sys)235 static int sys_max_predet_lag (const equation_system *sys)
236 {
237 int i, m = 0;
238
239 for (i=0; i<sys->plist[0]; i++) {
240 if (sys->pre_vars[i].lag > m) {
241 m = sys->pre_vars[i].lag;
242 }
243 }
244
245 return m;
246 }
247
get_predet_parent(const equation_system * sys,int v,int * lag)248 static int get_predet_parent (const equation_system *sys, int v, int *lag)
249 {
250 int i;
251
252 for (i=0; i<sys->plist[0]; i++) {
253 if (v == sys->pre_vars[i].id) {
254 *lag = sys->pre_vars[i].lag;
255 return sys->pre_vars[i].src;
256 }
257 }
258
259 *lag = 0;
260
261 return -1;
262 }
263
264 /**
265 * print_equation_system_info:
266 * @sys: gretl equation system.
267 * @dset: dataset information.
268 * @opt: use OPT_H for printing in a form designed to appear
269 * in the header when results of estimation are printed.
270 * @prn: printing struct.
271 *
272 * Prints details of an equation system to @prn.
273 */
274
275 void
print_equation_system_info(const equation_system * sys,const DATASET * dset,gretlopt opt,PRN * prn)276 print_equation_system_info (const equation_system *sys,
277 const DATASET *dset,
278 gretlopt opt, PRN *prn)
279 {
280 int header = (opt & OPT_H);
281 int i, vi, lag;
282
283 if (header && sys->name != NULL && !sys_anonymous(sys->name)) {
284 pprintf(prn, "%s %s\n", _("Equation system"), sys->name);
285 }
286
287 if (!header) {
288 for (i=0; i<sys->neqns; i++) {
289 print_system_equation(sys->lists[i], dset, prn);
290 }
291 }
292
293 for (i=0; i<sys->nidents; i++) {
294 print_system_identity(sys->idents[i], dset, opt, prn);
295 }
296
297 if (sys->ylist != NULL) {
298 pputs(prn, (header)? _("Endogenous variables:") : "endog");
299 for (i=1; i<=sys->ylist[0]; i++) {
300 vi = sys->ylist[i];
301 pprintf(prn, " %s", dset->varname[vi]);
302 }
303 pputc(prn, '\n');
304 }
305
306 if (header) {
307 if (sys->pre_vars != NULL) {
308 pputs(prn, _("Predetermined variables:"));
309 for (i=0; i<sys->plist[0]; i++) {
310 vi = sys->pre_vars[i].src;
311 lag = sys->pre_vars[i].lag;
312 pprintf(prn, " %s(-%d)", dset->varname[vi], lag);
313 }
314 pputc(prn, '\n');
315 }
316 if (sys->ilist != NULL && sys->ilist[0] > sys->plist[0]) {
317 pputs(prn, _("Exogenous variables:"));
318 for (i=1; i<=sys->ilist[0]; i++) {
319 vi = sys->ilist[i];
320 if (!in_gretl_list(sys->plist, vi)) {
321 pprintf(prn, " %s", dset->varname[vi]);
322 }
323 }
324 pputc(prn, '\n');
325 }
326 } else if (sys->ilist != NULL) {
327 pputs(prn, "instr");
328 for (i=1; i<=sys->ilist[0]; i++) {
329 vi = sys->ilist[i];
330 pprintf(prn, " %s", dset->varname[vi]);
331 }
332 pputc(prn, '\n');
333 }
334 }
335
system_method_from_string(const char * s)336 int system_method_from_string (const char *s)
337 {
338 int i = 0;
339
340 while (gretl_system_method_strings[i] != NULL) {
341 if (!strcmp(s, gretl_system_method_strings[i])) {
342 return i;
343 }
344 i++;
345 }
346
347 return i;
348 }
349
system_method_full_string(int method)350 const char *system_method_full_string (int method)
351 {
352 if (method >= SYS_METHOD_SUR && method < SYS_METHOD_MAX) {
353 return gretl_system_long_strings[method];
354 } else {
355 return NULL;
356 }
357 }
358
system_method_short_string(int method)359 const char *system_method_short_string (int method)
360 {
361 if (method >= SYS_METHOD_SUR && method < SYS_METHOD_MAX) {
362 return gretl_system_method_strings[method];
363 } else {
364 return NULL;
365 }
366 }
367
368 static equation_system *
equation_system_new(int method,const char * name,int * err)369 equation_system_new (int method, const char *name, int *err)
370 {
371 equation_system *sys;
372
373 if (method < 0 && name == NULL) {
374 *err = E_DATA;
375 return NULL;
376 }
377
378 sys = malloc(sizeof *sys);
379 if (sys == NULL) {
380 *err = E_ALLOC;
381 return NULL;
382 }
383
384 sys->name = NULL;
385
386 if (name != NULL) {
387 equation_system_set_name(sys, name);
388 }
389
390 sys->refcount = 0;
391 sys->fd = gretl_function_depth();
392 sys->method = method;
393
394 sys->t1 = sys->t2 = 0;
395 sys->smpl_t1 = sys->smpl_t2 = 0;
396 sys->T = sys->df = 0;
397
398 sys->neqns = 0;
399 sys->nidents = 0;
400
401 sys->q = sys->R = NULL;
402
403 sys->order = sys->flags = sys->iters = 0;
404
405 sys->ll = sys->llu = sys->ldet = NADBL;
406 sys->X2 = NADBL;
407 sys->ess = NADBL;
408 sys->diag_test = 0.0;
409 sys->bdiff = 0.0;
410
411 sys->b = sys->vcv = NULL;
412 sys->E = sys->S = NULL;
413 sys->yhat = NULL;
414
415 sys->Gamma = NULL;
416 sys->A = sys->B = NULL;
417 sys->F = sys->Sr = NULL;
418
419 sys->tslists = sys->lists = NULL;
420 sys->ilist = sys->ylist = NULL;
421 sys->plist = sys->xlist = NULL;
422 sys->biglist = NULL;
423
424 sys->pre_vars = NULL;
425 sys->idents = NULL;
426
427 sys->models = NULL;
428 sys->ldata = NULL;
429
430 return sys;
431 }
432
system_clear_restrictions(equation_system * sys)433 static void system_clear_restrictions (equation_system *sys)
434 {
435 if (sys->R != NULL) {
436 free(sys->R);
437 sys->R = NULL;
438 }
439
440 if (sys->q != NULL) {
441 free(sys->q);
442 sys->q = NULL;
443 }
444
445 sys->flags &= ~SYSTEM_RESTRICT;
446 }
447
system_clear_results(equation_system * sys)448 static void system_clear_results (equation_system *sys)
449 {
450 sys->iters = 0;
451
452 sys->t1 = sys->t2 = 0;
453 sys->smpl_t1 = sys->smpl_t2 = 0;
454 sys->T = sys->df = 0;
455
456 sys->ll = NADBL;
457 sys->llu = NADBL;
458 sys->ldet = NADBL;
459 sys->X2 = NADBL;
460 sys->ess = NADBL;
461
462 sys->diag_test = 0.0;
463 sys->bdiff = 0.0;
464
465 if (sys->b != NULL) {
466 gretl_matrix_replace(&sys->b, NULL);
467 }
468
469 if (sys->vcv != NULL) {
470 gretl_matrix_replace(&sys->vcv, NULL);
471 }
472
473 if (sys->S != NULL) {
474 gretl_matrix_replace(&sys->S, NULL);
475 }
476
477 if (sys->E != NULL) {
478 gretl_matrix_replace(&sys->E, NULL);
479 }
480
481 if (sys->yhat != NULL) {
482 gretl_matrix_replace(&sys->yhat, NULL);
483 }
484
485 if (sys->Gamma != NULL) {
486 gretl_matrix_replace(&sys->Gamma, NULL);
487 }
488
489 if (sys->B != NULL) {
490 gretl_matrix_replace(&sys->B, NULL);
491 }
492
493 if (sys->A != NULL) {
494 gretl_matrix_replace(&sys->A, NULL);
495 }
496
497 if (sys->Sr != NULL) {
498 gretl_matrix_replace(&sys->Sr, NULL);
499 }
500
501 if (sys->F != NULL) {
502 gretl_matrix_replace(&sys->F, NULL);
503 }
504
505 if (sys->ldata != NULL) {
506 free(sys->ldata->lmin);
507 free(sys->ldata->ll);
508 free(sys->ldata->idf);
509 free(sys->ldata);
510 sys->ldata = NULL;
511 }
512 }
513
equation_system_destroy(equation_system * sys)514 void equation_system_destroy (equation_system *sys)
515 {
516 int i;
517
518 #if SYSDEBUG
519 fprintf(stderr, "equation_system_destroy: %p\n", (void *) sys);
520 #endif
521
522 if (sys == NULL || sys->lists == NULL) {
523 return;
524 }
525
526 sys->refcount -= 1;
527 if (sys->refcount > 0) {
528 return;
529 }
530
531 for (i=0; i<sys->neqns; i++) {
532 free(sys->lists[i]);
533 }
534
535 free(sys->lists);
536 sys->lists = NULL;
537
538 for (i=0; i<sys->nidents; i++) {
539 destroy_ident(sys->idents[i]);
540 }
541
542 free(sys->idents);
543
544 free(sys->ylist);
545 free(sys->ilist);
546 free(sys->xlist);
547 free(sys->plist);
548 free(sys->biglist);
549
550 free(sys->pre_vars);
551
552 free(sys->name);
553
554 gretl_matrix_free(sys->R);
555 gretl_matrix_free(sys->q);
556
557 system_clear_results(sys);
558
559 free(sys);
560 }
561
562 #define SYS_PRUNE_LISTS 1
563
564 #if SYS_PRUNE_LISTS
565
566 /* 2020-12-09: When a system is defined using TSLS-style
567 regression lists, with instruments specified per
568 equation, there were problems estimating the system
569 via estimators that don't take such a specification.
570 The apparatus below is designed to handle this: we
571 remove the instruments immediately prior to estimation,
572 then restore them immediately afterwards.
573 */
574
575 /* Make a copy of the incoming per-equation regression
576 lists, with any instruments removed.
577 */
578
lists_strip_insts(equation_system * sys)579 static int lists_strip_insts (equation_system *sys)
580 {
581 int i, *eqlist;
582
583 sys->tslists = malloc(sys->neqns * sizeof *sys->tslists);
584
585 for (i=0; i<sys->neqns; i++) {
586 sys->tslists[i] = sys->lists[i];
587 if (gretl_list_has_separator(sys->lists[i])) {
588 gretl_list_split_on_separator(sys->tslists[i], &eqlist, NULL);
589 sys->lists[i] = eqlist;
590 } else {
591 sys->lists[i] = gretl_list_copy(sys->tslists[i]);
592 }
593 }
594
595 return 0;
596 }
597
598 /* Restore the original per-equation regression lists,
599 in case instruments were removed for estimation by
600 OLS, WLS, SUR or FIML.
601 */
602
lists_restore_insts(equation_system * sys)603 static void lists_restore_insts (equation_system *sys)
604 {
605 int i;
606
607 for (i=0; i<sys->neqns; i++) {
608 free(sys->lists[i]);
609 sys->lists[i] = sys->tslists[i];
610 sys->tslists[i] = NULL;
611 }
612
613 free(sys->tslists);
614 sys->tslists = NULL;
615 }
616
617 #endif /* SYS_PRUNE_LISTS */
618
sys_rearrange_eqn_lists(equation_system * sys,const DATASET * dset)619 static int sys_rearrange_eqn_lists (equation_system *sys,
620 const DATASET *dset)
621 {
622 int i, err = 0;
623
624 for (i=0; i<sys->neqns; i++) {
625 reglist_check_for_const(sys->lists[i], dset);
626 }
627
628 #if SYS_PRUNE_LISTS
629 if (sys->method != SYS_METHOD_TSLS &&
630 sys->method != SYS_METHOD_3SLS &&
631 sys->method != SYS_METHOD_LIML) {
632 /* we can't have ';' in equation lists */
633 for (i=0; i<sys->neqns; i++) {
634 if (gretl_list_has_separator(sys->lists[i])) {
635 lists_strip_insts(sys);
636 break;
637 }
638 }
639 }
640 #else
641 /* the status quo prior to 2020-12-09 */
642 if (sys->method != SYS_METHOD_TSLS &&
643 sys->method != SYS_METHOD_3SLS) {
644 /* we can't have ';' in equation lists */
645 int j;
646
647 for (i=0; i<sys->neqns && !err; i++) {
648 for (j=0; j<=sys->lists[i][0] && !err; j++) {
649 if (sys->lists[i][j] == LISTSEP) {
650 gretl_errmsg_sprintf("%s: tsls-style lists not supported",
651 gretl_system_short_strings[sys->method]);
652 err = E_DATA;
653 }
654 }
655 }
656 }
657 #endif
658
659 return err;
660 }
661
662 /* Form a list from row @i of matrix @m, ignoring trailing
663 zero elements. We check that series ID numbers are in
664 bounds, and also that the list does not contain
665 duplicated elements.
666 */
667
matrix_row_to_list(const gretl_matrix * m,int i,const DATASET * dset,int * err)668 static int *matrix_row_to_list (const gretl_matrix *m, int i,
669 const DATASET *dset,
670 int *err)
671 {
672 int *list = NULL;
673 int j, k, n = m->cols;
674
675 /* figure how many elements to read */
676 for (j=m->cols-1; j>=0; j--) {
677 if (gretl_matrix_get(m, i, j) == 0) {
678 n--;
679 } else {
680 break;
681 }
682 }
683
684 if (n == 0) {
685 *err = E_DATA;
686 return NULL;
687 }
688
689 /* check for out-of-bounds values */
690 for (j=0; j<n; j++) {
691 k = gretl_matrix_get(m, i, j);
692 if (k < 0 || k >= dset->v) {
693 *err = E_UNKVAR;
694 return NULL;
695 }
696 }
697
698 /* construct the list */
699 list = gretl_list_new(n);
700 if (list == NULL) {
701 *err = E_ALLOC;
702 } else {
703 k = 1;
704 for (j=0; j<n; j++) {
705 list[k++] = (int) gretl_matrix_get(m, i, j);
706 }
707 }
708
709 if (!*err) {
710 k = gretl_list_duplicates(list, EQUATION);
711 if (k >= 0) {
712 gretl_errmsg_sprintf(_("variable %d duplicated in the "
713 "command list."), k);
714 *err = E_DATA;
715 free(list);
716 list = NULL;
717 }
718 }
719
720 return list;
721 }
722
sys_check_eqn_list(const int * list)723 static int sys_check_eqn_list (const int *list)
724 {
725 int dupv = gretl_list_duplicates(list, EQUATION);
726
727 if (dupv >= 0) {
728 gretl_errmsg_sprintf(_("variable %d duplicated in the "
729 "command list."), dupv);
730 return E_DATA;
731 } else {
732 return 0;
733 }
734 }
735
736 /* @LY is a simple list of g regressands; either @LX is
737 simple list of common regressors or @AX is an array of
738 lists. We get either @LX or @AX as non-NULL, never both.
739 */
740
add_equations_from_lists(equation_system * sys,const int * LY,const int * LX,gretl_array * AX,const DATASET * dset)741 static int add_equations_from_lists (equation_system *sys,
742 const int *LY,
743 const int *LX,
744 gretl_array *AX,
745 const DATASET *dset)
746 {
747 int nx, n = sys->neqns;
748 int g = LY[0];
749 int i, j, err = 0;
750
751 if (AX != NULL && gretl_array_get_length(AX) != g) {
752 err = E_ARGS;
753 } else {
754 int **L = realloc(sys->lists, (n + g) * sizeof *sys->lists);
755
756 if (L == NULL) {
757 err = E_ALLOC;
758 } else {
759 sys->lists = L;
760 for (i=0; i<g; i++) {
761 sys->lists[n+i] = NULL;
762 }
763 }
764 }
765
766 if (err) {
767 return err;
768 }
769
770 if (AX != NULL) {
771 /* handle an array of regressor lists */
772 int *rhs, *list;
773
774 for (i=0; i<g && !err; i++) {
775 list = NULL;
776 rhs = gretl_array_get_element(AX, i, NULL, &err);
777 if (!err) {
778 nx = rhs[0];
779 list = gretl_list_new(1 + nx);
780 if (list == NULL) {
781 err = E_ALLOC;
782 }
783 }
784 if (!err) {
785 list[1] = LY[i+1];
786 for (j=1; j<=nx; j++) {
787 list[j+1] = rhs[j];
788 }
789 err = sys_check_eqn_list(list);
790 }
791 if (err) {
792 free(list);
793 } else {
794 sys->lists[n++] = list;
795 }
796 }
797 } else {
798 /* handle a list of common regressors */
799 int nx = LX[0];
800 int *list, *l0 = NULL;
801
802 for (i=0; i<g && !err; i++) {
803 if (i == 0) {
804 l0 = list = gretl_list_new(1 + nx);
805 } else {
806 list = gretl_list_copy(l0);
807 }
808 if (list == NULL) {
809 err = E_ALLOC;
810 } else {
811 list[1] = LY[i+1];
812 if (i == 0) {
813 for (j=1; j<=nx; j++) {
814 list[j+1] = LX[j];
815 }
816 }
817 err = sys_check_eqn_list(list);
818 }
819 if (err) {
820 free(list);
821 } else {
822 sys->lists[n++] = list;
823 }
824 }
825 }
826
827 if (!err) {
828 sys->neqns += g;
829 }
830
831 return err;
832 }
833
add_equations_from_matrix(equation_system * sys,const gretl_matrix * m,const DATASET * dset)834 static int add_equations_from_matrix (equation_system *sys,
835 const gretl_matrix *m,
836 const DATASET *dset)
837 {
838 int *list, **L;
839 int n = sys->neqns;
840 int i, err = 0;
841
842 L = realloc(sys->lists, (n + m->rows) * sizeof *sys->lists);
843 if (L == NULL) {
844 return E_ALLOC;
845 }
846
847 sys->lists = L;
848
849 for (i=0; i<m->rows && !err; i++) {
850 list = matrix_row_to_list(m, i, dset, &err);
851 if (!err) {
852 err = sys_check_eqn_list(list);
853 }
854 if (!err) {
855 sys->lists[n++] = list;
856 }
857 }
858
859 if (!err) {
860 sys->neqns += m->rows;
861 }
862
863 return err;
864 }
865
866 /**
867 * equation_system_append_multi:
868 * @sys: initialized equation system.
869 * @parm1: the name of a pre-defined matrix or list.
870 * @parm2: the name of list or array of lists (or %NULL).
871 * @dset: pointer to dataset struct.
872 *
873 * Adds one or more equations to @sys as follows.
874 *
875 * If @parm2 is %NULL then @parm1 is taken to be the
876 * name of a matrix, and we interpret the rows of the
877 * specified matrix as lists. Lists of differing length
878 * can be accommodated by padding unused trailing elements
879 * of short rows with zeros.
880 *
881 * If @parm2 is non-%NULL then @parm1 is taken to be the name
882 * of a list of regressands, one per equation; @parm2,
883 * pertaining to regressors, should be the name of a list
884 * if the regressors are in common across the equations
885 * or an array of lists (one per equation) otherwise.
886 *
887 * Returns: 0 on success, non-zero on failure, in which case
888 * @sys is destroyed.
889 */
890
equation_system_append_multi(equation_system * sys,const char * parm1,const char * parm2,const DATASET * dset)891 int equation_system_append_multi (equation_system *sys,
892 const char *parm1,
893 const char *parm2,
894 const DATASET *dset)
895 {
896 int err = 0;
897
898 if (sys == NULL) {
899 gretl_errmsg_set(_(nosystem));
900 return E_DATA;
901 } else if (parm1 == NULL) {
902 return E_ARGS;
903 }
904
905 if (parm2 != NULL) {
906 /* look for two lists, or list plus array */
907 const int *LY = get_list_by_name(parm1);
908 const int *LX = get_list_by_name(parm2);
909 gretl_array *AX = NULL;
910
911 if (LY == NULL) {
912 gretl_errmsg_sprintf(_("'%s': no such list"), parm1);
913 err = E_UNKVAR;
914 } else if (LX == NULL) {
915 /* try for an array on the right */
916 AX = get_array_by_name(parm2);
917 if (gretl_array_get_type(AX) != GRETL_TYPE_LISTS) {
918 gretl_errmsg_sprintf(_("'%s': not a list or array of lists"),
919 parm2);
920 err = E_UNKVAR;
921 }
922 }
923 if (!err) {
924 err = add_equations_from_lists(sys, LY, LX, AX, dset);
925 }
926 } else {
927 /* look for one matrix */
928 const gretl_matrix *m = get_matrix_by_name(parm1);
929
930 if (m == NULL) {
931 gretl_errmsg_sprintf(_("'%s': no such matrix"), parm1);
932 err = E_UNKVAR;
933 } else if (m->rows == 0 || m->cols == 0) {
934 err = E_DATA;
935 } else {
936 err = add_equations_from_matrix(sys, m, dset);
937 }
938 }
939
940 if (err) {
941 equation_system_destroy(sys);
942 }
943
944 return err;
945 }
946
947 /**
948 * equation_system_append:
949 * @sys: initialized equation system.
950 * @list: list containing dependent variable and regressors.
951 *
952 * Adds an equation, as represented by @list, to @sys.
953 *
954 * Returns: 0 on success, non-zero on failure, in which case
955 * @sys is destroyed.
956 */
957
equation_system_append(equation_system * sys,const int * list)958 int equation_system_append (equation_system *sys, const int *list)
959 {
960 int err = 0;
961
962 if (sys == NULL) {
963 gretl_errmsg_set(_(nosystem));
964 err = E_DATA;
965 } else {
966 int n = sys->neqns;
967 int **lists;
968
969 lists = realloc(sys->lists, (n + 1) * sizeof *sys->lists);
970 if (lists == NULL) {
971 err = E_ALLOC;
972 } else {
973 sys->lists = lists;
974 sys->lists[n] = gretl_list_copy(list);
975 #if SYSDEBUG
976 fprintf(stderr, "equation_system_append: added list %d\n", n);
977 printlist(list, "newly added list");
978 #endif
979 if (sys->lists[n] == NULL) {
980 err = E_ALLOC;
981 } else {
982 sys->neqns += 1;
983 }
984 }
985
986 if (err) {
987 equation_system_destroy(sys);
988 }
989 }
990
991 return err;
992 }
993
994 /* retrieve the name -- possibly quoted with embedded spaces -- for
995 an equation system */
996
get_system_name_from_line(const char * s)997 char *get_system_name_from_line (const char *s)
998 {
999 const char *p = NULL;
1000 char *name = NULL;
1001 int pchars = 0;
1002
1003 if (!strncmp(s, "method", 6)) {
1004 /* skip "method = whatever", with possible
1005 spaces around '='
1006 */
1007 char c = *(s+6);
1008
1009 if (c == ' ' || c == '=') {
1010 p = s + 6;
1011 p += strspn(s, " ");
1012 if (*p == '=') p++;
1013 p += strspn(s, " ");
1014 p += strcspn(s, " ");
1015 p += strspn(s, " ");
1016 s = p;
1017 }
1018 }
1019
1020 #if SYSDEBUG
1021 fprintf(stderr, "get_system_name_from_line: s = '%s'\n", s);
1022 #endif
1023
1024 if (*s == '"') {
1025 if (*(s + 1) != '\0') s++;
1026 p = s;
1027 while (*p && *p != '"') {
1028 if (!isspace((unsigned char) *p)) pchars++;
1029 p++;
1030 }
1031 if (*p != '"') {
1032 /* no closing quote */
1033 pchars = 0;
1034 }
1035 } else {
1036 p = s;
1037 while (*p && !isspace((unsigned char) *p)) {
1038 pchars++;
1039 p++;
1040 }
1041 }
1042
1043 if (pchars > 0) {
1044 name = gretl_strndup(s, p - s);
1045 }
1046
1047 return name;
1048 }
1049
1050 /* parse a system estimation method out of parameter @s */
1051
sys_get_estimator(const char * s)1052 static int sys_get_estimator (const char *s)
1053 {
1054 const char *p = strstr(s, "method");
1055 int method = -1;
1056
1057 if (p != NULL) {
1058 p += 6;
1059 p += strspn(s, " ");
1060 if (*p == '=') {
1061 char mstr[9];
1062
1063 p++;
1064 p += strspn(s, " ");
1065 if (sscanf(p, "%8s", mstr) == 1) {
1066 gretl_lower(mstr);
1067 method = system_method_from_string(mstr);
1068 }
1069 }
1070 }
1071
1072 return method;
1073 }
1074
parse_sys_start_param(const char * s,int * method,char ** name)1075 static int parse_sys_start_param (const char *s,
1076 int *method,
1077 char **name)
1078 {
1079 int len = strlen(s);
1080 int err = 0;
1081
1082 if (len > 7 && !strncmp(s, "method=", 7)) {
1083 char mstr[9];
1084
1085 *mstr = '\0';
1086 strncat(mstr, s + 7, 8);
1087 gretl_lower(mstr);
1088 *method = system_method_from_string(mstr);
1089 if (*method == SYS_METHOD_MAX) {
1090 /* invalid method was given */
1091 gretl_errmsg_set(_(badsystem));
1092 err = E_DATA;
1093 }
1094 } else if (len > 5 && !strncmp(s, "name=", 5)) {
1095 *name = gretl_strdup(s + 5);
1096 gretl_unquote(*name, &err);
1097 } else {
1098 err = E_PARSE;
1099 }
1100
1101 return err;
1102 }
1103
1104 /**
1105 * equation_system_start:
1106 * @param: may incude an estimation method spec.
1107 * @name: the name to be given to system, if any, otherwise NULL
1108 * or an empty string.
1109 * @opt: may include OPT_I for iterative estimation (will be
1110 * ignored if the the estimation method does not support it).
1111 * @err: location to receive error code.
1112 *
1113 * Start compiling an equation system. Either @param must contain
1114 * an estimation method or the system must be given a name.
1115 * If a method is given, the system will be estimated as soon as its
1116 * definition is complete. If a name is given, the system definition
1117 * is saved on a stack and it can subsequently be estimated via various
1118 * methods. If both a name and an estimation method are given, the system
1119 * is both estimated and saved.
1120 *
1121 * Returns: pointer to a new equation system, or %NULL on error.
1122 */
1123
equation_system_start(const char * param,const char * name,gretlopt opt,int * err)1124 equation_system *equation_system_start (const char *param,
1125 const char *name,
1126 gretlopt opt,
1127 int *err)
1128 {
1129 equation_system *sys = NULL;
1130 char *oldname = NULL;
1131 int method = -1;
1132
1133 #if SYSDEBUG > 1
1134 fprintf(stderr, "equation_system_start\n");
1135 #endif
1136
1137 if (param != NULL) {
1138 *err = parse_sys_start_param(param, &method, &oldname);
1139 if (*err) {
1140 return NULL;
1141 }
1142 }
1143
1144 if (name != NULL && *name != '\0') {
1145 sys = equation_system_new(method, name, err);
1146 } else if (oldname != NULL) {
1147 sys = equation_system_new(method, oldname, err);
1148 free(oldname);
1149 } else {
1150 /* treat the system as "anonymous" */
1151 sys = equation_system_new(method, "$system", err);
1152 if (!*err) {
1153 push_anon_system(sys);
1154 }
1155 }
1156
1157 if (sys != NULL) {
1158 if (opt & OPT_I) {
1159 sys->flags |= SYSTEM_ITERATE;
1160 }
1161 if (opt & OPT_Q) {
1162 sys->flags |= SYSTEM_QUIET;
1163 }
1164 }
1165
1166 #if SYSDEBUG > 1
1167 fprintf(stderr, "new system '%s' at %p, flags = %d\n", sys->name,
1168 (void *) sys, sys->flags);
1169 #endif
1170
1171 return sys;
1172 }
1173
1174 /* determine the degrees of freedom for the unrestricted
1175 estimation of @sys
1176 */
1177
system_get_dfu(const equation_system * sys)1178 static int system_get_dfu (const equation_system *sys)
1179 {
1180 int dfu = sys->T * sys->neqns; /* total observations */
1181 int i, pos;
1182
1183 for (i=0; i<sys->neqns; i++) {
1184 /* subtract the number of parameters */
1185 pos = gretl_list_separator_position(sys->lists[i]);
1186 if (pos > 0) {
1187 dfu -= pos - 2;
1188 } else {
1189 dfu -= sys->lists[i][0] - 1;
1190 }
1191 }
1192
1193 return dfu;
1194 }
1195
get_eqn_ref(const equation_system * sys,int j)1196 static int get_eqn_ref (const equation_system *sys, int j)
1197 {
1198 int i, pos, nparm = 0;
1199
1200 for (i=0; i<sys->neqns; i++) {
1201 pos = gretl_list_separator_position(sys->lists[i]);
1202 if (pos > 0) {
1203 nparm += pos - 2;
1204 } else {
1205 nparm += sys->lists[i][0] - 1;
1206 }
1207 if (nparm > j) {
1208 return i;
1209 }
1210 }
1211
1212 return -1;
1213 }
1214
1215 /* Determine whether the restrictions on @sys embodied in @R
1216 reference just one equation, and if so return the degrees
1217 of freedom for that equation. Otherwise return 0.
1218 */
1219
maybe_get_single_equation_dfu(const equation_system * sys,const gretl_matrix * R)1220 static int maybe_get_single_equation_dfu (const equation_system *sys,
1221 const gretl_matrix *R)
1222 {
1223 int i, j, eq = -1, eq0 = -1;
1224 int df = 0;
1225
1226 for (j=0; j<R->cols; j++) {
1227 for (i=0; i<R->rows; i++) {
1228 if (gretl_matrix_get(R, i, j) != 0) {
1229 eq = get_eqn_ref(sys, j);
1230 if (eq0 < 0) {
1231 eq0 = eq;
1232 } else if (eq != eq0) {
1233 /* multiple equations are referenced */
1234 return 0;
1235 }
1236 }
1237 }
1238 }
1239
1240 if (eq >= 0) {
1241 df = sys->T - (sys->lists[eq][0] - 1);
1242 }
1243
1244 return df;
1245 }
1246
1247 /* Asymptotic F-test, as in William Greene:
1248
1249 (1/J) * (Rb-q)' * [R*Var(b)*R']^{-1} * (Rb-q)
1250
1251 or Chi-square version if dfu = 0
1252 */
1253
multi_eqn_wald_test(const gretl_matrix * b,const gretl_matrix * V,const gretl_matrix * R,const gretl_matrix * q,int dfu,gretlopt opt,PRN * prn)1254 int multi_eqn_wald_test (const gretl_matrix *b,
1255 const gretl_matrix *V,
1256 const gretl_matrix *R,
1257 const gretl_matrix *q,
1258 int dfu, gretlopt opt,
1259 PRN *prn)
1260 {
1261 gretl_matrix *Rbq, *RvR;
1262 int Rrows, dfn;
1263 double test = NADBL;
1264 int err = 0;
1265
1266 if (R == NULL || q == NULL || b == NULL || V == NULL) {
1267 pputs(prn, "Missing matrix in system Wald test!\n");
1268 return E_DATA;
1269 }
1270
1271 Rrows = gretl_matrix_rows(R);
1272 dfn = gretl_matrix_rows(R);
1273
1274 Rbq = gretl_matrix_alloc(Rrows, 1);
1275 RvR = gretl_matrix_alloc(Rrows, Rrows);
1276
1277 if (Rbq == NULL || RvR == NULL) {
1278 err = E_ALLOC;
1279 } else {
1280 gretl_matrix_multiply(R, b, Rbq);
1281 gretl_matrix_subtract_from(Rbq, q);
1282 gretl_matrix_qform(R, GRETL_MOD_NONE, V,
1283 RvR, GRETL_MOD_NONE);
1284 err = gretl_invert_symmetric_matrix(RvR);
1285 }
1286
1287 if (!err) {
1288 test = gretl_scalar_qform(Rbq, RvR, &err);
1289 }
1290
1291 if (!err) {
1292 double pval;
1293
1294 if (dfu == 0) {
1295 pval = chisq_cdf_comp(dfn, test);
1296 } else {
1297 test /= dfn;
1298 pval = snedecor_cdf_comp(dfn, dfu, test);
1299 }
1300
1301 record_test_result(test, pval);
1302
1303 if (!(opt & OPT_Q)) {
1304 if (dfu == 0) {
1305 pprintf(prn, "%s:\n", _("Wald test for the specified restrictions"));
1306 pprintf(prn, " %s(%d) = %g [%.4f]\n", _("Chi-square"),
1307 dfn, test, pval);
1308 } else {
1309 pprintf(prn, "%s:\n", _("F test for the specified restrictions"));
1310 pprintf(prn, " F(%d,%d) = %g [%.4f]\n", dfn, dfu, test, pval);
1311 }
1312 pputc(prn, '\n');
1313 }
1314 }
1315
1316 gretl_matrix_free(Rbq);
1317 gretl_matrix_free(RvR);
1318
1319 return err;
1320 }
1321
get_wald_dfu(const equation_system * sys,const gretl_matrix * R)1322 static int get_wald_dfu (const equation_system *sys,
1323 const gretl_matrix *R)
1324 {
1325 int dfu = 0;
1326
1327 if (sys->method == SYS_METHOD_OLS || sys->method == SYS_METHOD_TSLS) {
1328 dfu = maybe_get_single_equation_dfu(sys, R);
1329 }
1330 if (dfu == 0) {
1331 dfu = system_get_dfu(sys);
1332 }
1333
1334 return dfu;
1335 }
1336
system_wald_test(const equation_system * sys,const gretl_matrix * R,const gretl_matrix * q,gretlopt opt,PRN * prn)1337 int system_wald_test (const equation_system *sys,
1338 const gretl_matrix *R,
1339 const gretl_matrix *q,
1340 gretlopt opt,
1341 PRN *prn)
1342 {
1343 const gretl_matrix *b = sys->b;
1344 const gretl_matrix *V = sys->vcv;
1345 int dfu = 0;
1346
1347 if (b == NULL || V == NULL) {
1348 gretl_errmsg_set("restrict: no estimates are available");
1349 return E_DATA;
1350 }
1351
1352 if (sys->method == SYS_METHOD_OLS || sys->method == SYS_METHOD_TSLS) {
1353 /* use the F-form? */
1354 dfu = maybe_get_single_equation_dfu(sys, R);
1355 }
1356
1357 return multi_eqn_wald_test(b, V, R, q, dfu, opt, prn);
1358 }
1359
system_do_LR_test(const equation_system * sys,double llu,gretlopt opt,PRN * prn)1360 static int system_do_LR_test (const equation_system *sys,
1361 double llu, gretlopt opt,
1362 PRN *prn)
1363 {
1364 double X2, pval, llr = sys->ll;
1365 int df = gretl_matrix_rows(sys->R);
1366 int err = 0;
1367
1368 if (na(llr) || na(llu) || llr == 0.0 || llu == 0.0 || df <= 0) {
1369 fputs("bad or missing data in system LR test\n", stderr);
1370 return E_DATA;
1371 }
1372
1373 X2 = 2.0 * (llu - llr);
1374 pval = chisq_cdf_comp(df, X2);
1375
1376 if (!(opt & OPT_Q)) {
1377 pprintf(prn, "%s:\n", _("LR test for the specified restrictions"));
1378 pprintf(prn, " %s = %g\n", _("Restricted log-likelihood"), llr);
1379 pprintf(prn, " %s = %g\n", _("Unrestricted log-likelihood"), llu);
1380 pprintf(prn, " %s(%d) = %g [%.4f]\n", _("Chi-square"), df, X2, pval);
1381 pputc(prn, '\n');
1382 }
1383
1384 return err;
1385 }
1386
sys_test_type(equation_system * sys)1387 static int sys_test_type (equation_system *sys)
1388 {
1389 int ret = SYS_TEST_NONE;
1390
1391 if (system_n_restrictions(sys) > 0) {
1392 if (sys->method == SYS_METHOD_SUR ||
1393 sys->method == SYS_METHOD_WLS) {
1394 if (sys->flags & SYSTEM_ITERATE) {
1395 ret = SYS_TEST_LR;
1396 } else {
1397 ret = SYS_TEST_F;
1398 }
1399 } else if (sys->method == SYS_METHOD_OLS ||
1400 sys->method == SYS_METHOD_TSLS ||
1401 sys->method == SYS_METHOD_3SLS) {
1402 ret = SYS_TEST_F;
1403 } else if (sys->method == SYS_METHOD_LIML) {
1404 /* experimental */
1405 ret = SYS_TEST_F;
1406 } else if (sys->method == SYS_METHOD_FIML) {
1407 ret = SYS_TEST_LR;
1408 }
1409 }
1410
1411 return ret;
1412 }
1413
1414 /* Handle the case where sys->b and sys->vcv have been augmented in
1415 order to test some restrictions: we now have to trim these matrices
1416 down to their final size. The vector @b contains unrestricted
1417 coefficient estimates, so its length can be used to gauge the
1418 true number of coeffs.
1419 */
1420
shrink_b_and_vcv(const gretl_matrix * b,equation_system * sys)1421 static int shrink_b_and_vcv (const gretl_matrix *b,
1422 equation_system *sys)
1423 {
1424 int nc = gretl_vector_get_length(b);
1425 gretl_matrix *V;
1426 double x;
1427 int i, j;
1428
1429 if (sys->vcv == NULL) {
1430 return E_DATA;
1431 }
1432
1433 if (sys->vcv->rows == nc) {
1434 /* no-op (shouldn't happen) */
1435 return 0;
1436 }
1437
1438 V = gretl_matrix_alloc(nc, nc);
1439 if (V == NULL) {
1440 return E_ALLOC;
1441 }
1442
1443 gretl_matrix_reuse(sys->b, nc, 1);
1444
1445 for (i=0; i<nc; i++) {
1446 for (j=0; j<nc; j++) {
1447 x = gretl_matrix_get(sys->vcv, i, j);
1448 gretl_matrix_set(V, i, j, x);
1449 }
1450 }
1451
1452 gretl_matrix_replace(&sys->vcv, V);
1453
1454 return 0;
1455 }
1456
estimate_with_test(equation_system * sys,DATASET * dset,int stest,int (* system_est)(),gretlopt opt,PRN * prn)1457 static int estimate_with_test (equation_system *sys, DATASET *dset,
1458 int stest, int (*system_est)(),
1459 gretlopt opt, PRN *prn)
1460 {
1461 gretl_matrix *vcv = NULL;
1462 gretl_matrix *b = NULL;
1463 double llu = 0.0;
1464 int err = 0;
1465
1466 /* estimate the unrestricted system first */
1467 sys->flags &= ~SYSTEM_RESTRICT;
1468 err = (* system_est) (sys, dset, opt | OPT_Q, prn);
1469 sys->flags ^= SYSTEM_RESTRICT;
1470 if (err) {
1471 goto bailout;
1472 }
1473
1474 /* grab the data from unrestricted estimation */
1475 b = sys->b;
1476 sys->b = NULL;
1477 if (stest == SYS_TEST_LR) {
1478 llu = sys->ll;
1479 } else if (stest == SYS_TEST_F) {
1480 vcv = sys->vcv;
1481 sys->vcv = NULL;
1482 }
1483
1484 /* now estimate the restricted system */
1485 system_clear_results(sys);
1486 err = (* system_est) (sys, dset, opt, prn);
1487 if (!err) {
1488 if (stest == SYS_TEST_LR) {
1489 err = system_do_LR_test(sys, llu, opt, prn);
1490 } else if (stest == SYS_TEST_F) {
1491 int dfu = get_wald_dfu(sys, sys->R);
1492
1493 err = multi_eqn_wald_test(b, vcv, sys->R, sys->q,
1494 dfu, opt, prn);
1495 }
1496 shrink_b_and_vcv(b, sys);
1497 }
1498
1499 bailout:
1500
1501 if (b != NULL) gretl_matrix_free(b);
1502 if (vcv != NULL) gretl_matrix_free(vcv);
1503
1504 return err;
1505 }
1506
1507 static void
adjust_sys_flags_for_method(equation_system * sys,int method)1508 adjust_sys_flags_for_method (equation_system *sys, int method)
1509 {
1510 char oldflags = sys->flags;
1511
1512 sys->flags = 0;
1513
1514 /* the robust option is for now OLS-only */
1515 if (oldflags & SYSTEM_ROBUST) {
1516 if (sys->method == SYS_METHOD_OLS) {
1517 sys->flags |= SYSTEM_ROBUST;
1518 }
1519 }
1520
1521 /* the iterate option is only available for WLS, SUR or 3SLS */
1522 if (oldflags & SYSTEM_ITERATE) {
1523 if (sys->method == SYS_METHOD_WLS ||
1524 sys->method == SYS_METHOD_SUR ||
1525 sys->method == SYS_METHOD_3SLS) {
1526 sys->flags |= SYSTEM_ITERATE;
1527 }
1528 }
1529
1530 /* by default, apply df correction for single-equation methods */
1531 if (sys->method == SYS_METHOD_OLS ||
1532 sys->method == SYS_METHOD_WLS ||
1533 sys->method == SYS_METHOD_TSLS ||
1534 sys->method == SYS_METHOD_LIML) {
1535 if (oldflags & SYSTEM_DFCORR) {
1536 sys->flags |= SYSTEM_DFCORR;
1537 }
1538 }
1539
1540 /* carry forward the GEOMEAN flag */
1541 if (oldflags & SYSTEM_VCV_GEOMEAN) {
1542 sys->flags |= SYSTEM_VCV_GEOMEAN;
1543 }
1544 }
1545
1546 static int
set_sys_flags_from_opt(equation_system * sys,gretlopt opt)1547 set_sys_flags_from_opt (equation_system *sys, gretlopt opt)
1548 {
1549 sys->flags = 0;
1550
1551 /* the robust option is for now OLS-only */
1552 if (opt & OPT_R) {
1553 if (sys->method == SYS_METHOD_OLS) {
1554 sys->flags |= SYSTEM_ROBUST;
1555 } else {
1556 gretl_errmsg_sprintf(_("Invalid option '--%s'"), "robust");
1557 return E_BADOPT;
1558 }
1559 }
1560
1561 /* the iterate option is available for WLS, SUR or 3SLS */
1562 if (opt & OPT_I) {
1563 if (sys->method == SYS_METHOD_WLS ||
1564 sys->method == SYS_METHOD_SUR ||
1565 sys->method == SYS_METHOD_3SLS) {
1566 sys->flags |= SYSTEM_ITERATE;
1567 }
1568 }
1569
1570 /* by default, apply a df correction for single-equation methods */
1571 if (sys->method == SYS_METHOD_OLS ||
1572 sys->method == SYS_METHOD_TSLS) {
1573 if (!(opt & OPT_N)) {
1574 sys->flags |= SYSTEM_DFCORR;
1575 }
1576 }
1577
1578 if (opt & OPT_M) {
1579 sys->flags |= SYSTEM_VCV_GEOMEAN;
1580 }
1581 if (opt & OPT_Q) {
1582 sys->flags |= SYSTEM_QUIET;
1583 }
1584 if (opt & OPT_S) {
1585 /* estimating single equation via LIML */
1586 sys->flags |= SYSTEM_LIML1;
1587 }
1588
1589 return 0;
1590 }
1591
1592 /**
1593 * equation_system_estimate:
1594 * @sys: pre-defined equation system.
1595 * @dset: dataset struct.
1596 * @opt: may include OPT_V for more verbose operation.
1597 * @prn: printing struct.
1598 *
1599 * Estimate a pre-defined equation system and print the results
1600 * to @prn.
1601 *
1602 * Returns: 0 on success, non-zero on error.
1603 */
1604
1605 int
equation_system_estimate(equation_system * sys,DATASET * dset,gretlopt opt,PRN * prn)1606 equation_system_estimate (equation_system *sys, DATASET *dset,
1607 gretlopt opt, PRN *prn)
1608 {
1609 int (*system_est) (equation_system *, DATASET *,
1610 gretlopt, PRN *);
1611 int err = 0;
1612
1613 #if SYSDEBUG
1614 fprintf(stderr, "*** equation_system_estimate\n");
1615 #endif
1616
1617 gretl_error_clear();
1618
1619 if (sys->xlist == NULL || sys->biglist == NULL) {
1620 /* allow for the possibility that we're looking at a
1621 system restored from a session file
1622 */
1623 err = sys_check_lists(sys, dset);
1624 if (err) {
1625 return err;
1626 }
1627 }
1628
1629 if (opt == OPT_UNSET) {
1630 opt = OPT_NONE;
1631 adjust_sys_flags_for_method(sys, sys->method);
1632 } else {
1633 err = set_sys_flags_from_opt(sys, opt);
1634 if (err) {
1635 return err;
1636 }
1637 }
1638
1639 /* if restrictions are in place, reset the restricted flag */
1640 if (sys->R != NULL && sys->q != NULL) {
1641 sys->flags |= SYSTEM_RESTRICT;
1642 }
1643
1644 /* in case we're re-estimating */
1645 system_clear_results(sys);
1646 err = sys_rearrange_eqn_lists(sys, dset);
1647
1648 if (!err) {
1649 system_est = get_plugin_function("system_estimate");
1650 if (system_est == NULL) {
1651 err = 1;
1652 }
1653 }
1654
1655 if (!err) {
1656 int stest = sys_test_type(sys);
1657
1658 if (stest == SYS_TEST_NOTIMP) {
1659 pputs(prn, _("Sorry, command not available for this estimator"));
1660 pputc(prn, '\n');
1661 err = 1;
1662 } else if (stest != SYS_TEST_NONE) {
1663 err = estimate_with_test(sys, dset, stest,
1664 system_est, opt, prn);
1665 } else {
1666 err = (*system_est) (sys, dset, opt, prn);
1667 }
1668 }
1669
1670 if (sys->tslists != NULL) {
1671 lists_restore_insts(sys);
1672 }
1673
1674 if (!err) {
1675 sys->smpl_t1 = dset->t1;
1676 sys->smpl_t2 = dset->t2;
1677 if (!(sys->flags & SYSTEM_LIML1)) {
1678 set_as_last_model(sys, GRETL_OBJ_SYS);
1679 }
1680 }
1681
1682 return err;
1683 }
1684
system_adjust_t1t2(equation_system * sys,const DATASET * dset)1685 int system_adjust_t1t2 (equation_system *sys, const DATASET *dset)
1686 {
1687 int err;
1688
1689 if (sys->biglist == NULL) {
1690 fprintf(stderr, "system_adjust_t1t2: no 'biglist' present!\n");
1691 return E_DATA;
1692 }
1693
1694 sys->t1 = dset->t1;
1695 sys->t2 = dset->t2;
1696
1697 err = list_adjust_sample(sys->biglist, &sys->t1, &sys->t2, dset, NULL);
1698 if (!err) {
1699 sys->T = sys->t2 - sys->t1 + 1;
1700 }
1701
1702 return err;
1703 }
1704
1705 /* construct a list containing all the variables referenced in the
1706 system */
1707
sys_make_biglist(equation_system * sys)1708 static int sys_make_biglist (equation_system *sys)
1709 {
1710 int *biglist;
1711 const int *slist;
1712 const identity *ident;
1713 int bign = 0;
1714 int i, j, vj, k;
1715
1716 for (i=0; i<sys->neqns; i++) {
1717 bign += sys->lists[i][0] - gretl_list_has_separator(sys->lists[i]);
1718 }
1719
1720 for (i=0; i<sys->nidents; i++) {
1721 bign += 1 + sys->idents[i]->n_atoms;
1722 }
1723
1724 if (sys->ilist != NULL) {
1725 bign += sys->ilist[0];
1726 }
1727
1728 biglist = gretl_list_new(bign);
1729 if (biglist == NULL) {
1730 return E_ALLOC;
1731 }
1732
1733 for (i=1; i<=bign; i++) {
1734 /* invalidate all elements */
1735 biglist[i] = -1;
1736 }
1737
1738 k = 0;
1739
1740 /* system equations */
1741 for (i=0; i<sys->neqns; i++) {
1742 slist = sys->lists[i];
1743 for (j=1; j<=slist[0]; j++) {
1744 vj = slist[j];
1745 if (vj != LISTSEP && !in_gretl_list(biglist, vj)) {
1746 biglist[++k] = vj;
1747 }
1748 }
1749 }
1750
1751 /* system identities */
1752 for (i=0; i<sys->nidents; i++) {
1753 ident = sys->idents[i];
1754 for (j=0; j<=ident->n_atoms; j++) {
1755 if (j == 0) {
1756 vj = ident->depvar;
1757 } else {
1758 vj = ident->atoms[j-1].varnum;
1759 }
1760 if (!in_gretl_list(biglist, vj)) {
1761 biglist[++k] = vj;
1762 }
1763 }
1764 }
1765
1766 if (sys->ilist != NULL) {
1767 /* system instruments */
1768 for (j=1; j<=sys->ilist[0]; j++) {
1769 vj = sys->ilist[j];
1770 if (!in_gretl_list(biglist, vj)) {
1771 biglist[++k] = vj;
1772 }
1773 }
1774 }
1775
1776 biglist[0] = k;
1777
1778 #if SYSDEBUG
1779 printlist(biglist, "system biglist (all vars)");
1780 #endif
1781
1782 free(sys->biglist);
1783 sys->biglist = biglist;
1784
1785 return 0;
1786 }
1787
sys_get_lag_src(const char * vname,const DATASET * dset)1788 static int sys_get_lag_src (const char *vname, const DATASET *dset)
1789 {
1790 int fd, src = series_index(dset, vname);
1791
1792 if (src == dset->v && vname != NULL &&
1793 (fd = gretl_function_depth()) > 0) {
1794 int i;
1795
1796 for (i=1; i<dset->v; i++) {
1797 if (fd == series_get_stack_level(dset, i) &&
1798 series_is_listarg(dset, i, NULL) &&
1799 !strcmp(dset->varname[i], vname)) {
1800 src = i;
1801 break;
1802 }
1803 }
1804 }
1805
1806 return src;
1807 }
1808
is_tsls_style_instrument(equation_system * sys,int v)1809 static int is_tsls_style_instrument (equation_system *sys, int v)
1810 {
1811 int i, j, pos;
1812
1813 for (i=0; i<sys->neqns; i++) {
1814 pos = gretl_list_separator_position(sys->lists[i]);
1815 if (pos > 0) {
1816 for (j=pos+1; j<=sys->lists[i][0]; j++) {
1817 if (sys->lists[i][j] == v) {
1818 return 1;
1819 }
1820 }
1821 }
1822 }
1823
1824 return 0;
1825 }
1826
1827 /* Here we deal with the case where the user has specified the
1828 system equations "tsls-style": besides the left-hand side
1829 variables, we also need to identify as endogenous any
1830 regressors that do not appear as instruments.
1831 */
1832
tsls_style_augment_ylist(equation_system * sys,int ** pylist)1833 static int tsls_style_augment_ylist (equation_system *sys,
1834 int **pylist)
1835 {
1836 const int *list;
1837 int pos, ny = (*pylist)[0];
1838 int i, j, vj;
1839
1840 for (i=0; i<sys->neqns; i++) {
1841 list = sys->lists[i];
1842 pos = gretl_list_separator_position(list);
1843 if (pos == 0) {
1844 continue;
1845 }
1846 for (j=2; j<pos; j++) {
1847 vj = list[j];
1848 if (vj == 0 || in_gretl_list(*pylist, vj)) {
1849 continue;
1850 }
1851 if (!is_tsls_style_instrument(sys, vj)) {
1852 if (ny < sys->neqns) {
1853 (*pylist)[++ny] = vj;
1854 (*pylist)[0] = ny;
1855 } else {
1856 gretl_list_append_term(pylist, vj);
1857 if (*pylist == NULL) {
1858 return E_ALLOC;
1859 }
1860 ny++;
1861 }
1862 }
1863 }
1864 }
1865
1866 return ny >= sys->neqns ? 0 : E_DATA;
1867 }
1868
1869 /* It's possible to specify instruments per equation, as in tsls,
1870 with a two-part regression list. But this cannot be mixed with
1871 any other means of specifying endogeneity/exogeneity of
1872 variables.
1873 */
1874
check_for_tsls_style_lists(equation_system * sys,int * err)1875 static int check_for_tsls_style_lists (equation_system *sys,
1876 int *err)
1877 {
1878 int i, tsls_style = 0;
1879
1880 for (i=0; i<sys->neqns; i++) {
1881 if (gretl_list_has_separator(sys->lists[i])) {
1882 tsls_style = 1;
1883 break;
1884 }
1885 }
1886
1887 if (tsls_style) {
1888 if (sys->ylist != NULL || sys->ilist != NULL) {
1889 gretl_errmsg_set("You can't mix tsls-style lists with 'endog' or 'instr'");
1890 *err = E_DATA;
1891 } else if (sys->nidents > 0) {
1892 gretl_errmsg_set("You can't mix identities with tsls-style lists");
1893 *err = E_DATA;
1894 }
1895 }
1896
1897 return tsls_style;
1898 }
1899
1900 /* Given a "partial" system, with equations expressed in tsls
1901 style and with more endogenous variables than equations, shift
1902 the "excess" endogenous terms into @xplist. They will still
1903 be instrumented when the equations are estimated, but they
1904 will be treated "as if" exogenous when it comes to building
1905 the structural form of the system.
1906 */
1907
tsls_style_shift_vars(equation_system * sys,int * xplist)1908 static void tsls_style_shift_vars (equation_system *sys, int *xplist)
1909 {
1910 int i, vi, nxp = xplist[0];
1911
1912 for (i=sys->ylist[0]; i>sys->neqns; i--) {
1913 vi = sys->ylist[i];
1914 gretl_list_delete_at_pos(sys->ylist, i);
1915 xplist[++nxp] = vi;
1916 }
1917
1918 xplist[0] = nxp;
1919 }
1920
1921 /* For consistency with the way in which per-equation results
1922 are organized, we should ensure that if a constant is present
1923 among the exogenous terms in a system it appears as the
1924 first element of the system's "xlist".
1925 */
1926
sys_xlist_reshuffle_const(int * list)1927 static void sys_xlist_reshuffle_const (int *list)
1928 {
1929 int i, cpos = 0;
1930
1931 for (i=1; i<=list[0]; i++) {
1932 if (list[i] == 0) {
1933 cpos = i;
1934 break;
1935 }
1936 }
1937
1938 if (cpos > 1) {
1939 for (i=cpos; i>1; i--) {
1940 list[i] = list[i-1];
1941 }
1942 list[1] = 0;
1943 }
1944 }
1945
1946 /* prior to system estimation, get all the required lists of variables
1947 in order
1948 */
1949
sys_check_lists(equation_system * sys,const DATASET * dset)1950 static int sys_check_lists (equation_system *sys,
1951 const DATASET *dset)
1952 {
1953 const int *slist;
1954 const char *vname;
1955 const identity *ident;
1956 int user_ylist = (sys->ylist != NULL);
1957 int user_ilist = (sys->ilist != NULL);
1958 int *ylist = NULL;
1959 int *xplist = NULL;
1960 int src, lag, nlhs;
1961 int tsls_style;
1962 int i, j, k, vj;
1963 int err = 0;
1964
1965 #if SYSDEBUG
1966 fprintf(stderr, "*** sys_check_lists\n");
1967 #endif
1968
1969 tsls_style = check_for_tsls_style_lists(sys, &err);
1970 if (err) {
1971 return err;
1972 }
1973
1974 /* start an empty list for predetermined variables */
1975
1976 sys->plist = gretl_null_list();
1977 if (sys->plist == NULL) {
1978 return E_ALLOC;
1979 }
1980
1981 /* Compose a minimal (and possibly incomplete) list of endogenous
1982 variables, including all vars that appear on the left-hand side
1983 of structural equations or identities.
1984 */
1985
1986 nlhs = sys->neqns + sys->nidents;
1987 ylist = gretl_list_new(nlhs);
1988 if (ylist == NULL) {
1989 return E_ALLOC;
1990 }
1991
1992 k = 0;
1993
1994 for (i=0; i<sys->neqns; i++) {
1995 #if SYSDEBUG
1996 printlist(sys->lists[i], "incoming equation list");
1997 #endif
1998 slist = sys->lists[i];
1999 vj = slist[1];
2000 if (!in_gretl_list(ylist, vj)) {
2001 ylist[++k] = vj;
2002 }
2003 }
2004
2005 if (tsls_style) {
2006 ylist[0] = k;
2007 err = tsls_style_augment_ylist(sys, &ylist);
2008 if (err) {
2009 goto bailout;
2010 }
2011 } else {
2012 for (i=0; i<sys->nidents; i++) {
2013 ident = sys->idents[i];
2014 vj = ident->depvar;
2015 if (!in_gretl_list(ylist, vj)) {
2016 ylist[++k] = vj;
2017 }
2018 }
2019 ylist[0] = k;
2020 }
2021
2022 #if SYSDEBUG
2023 printlist(ylist, "system auto ylist");
2024 #endif
2025
2026 /* If the user gave a list of endogenous vars (in sys->ylist), check
2027 that it contains all the presumably endogenous variables we found
2028 above and recorded in ylist; otherwise use the ylist as computed
2029 above.
2030 */
2031 if (user_ylist) {
2032 for (j=1; j<=ylist[0]; j++) {
2033 vj = ylist[j];
2034 if (!in_gretl_list(sys->ylist, vj)) {
2035 gretl_errmsg_sprintf("%s appears on the left-hand side "
2036 "of an equation but is not marked as endogenous",
2037 dset->varname[vj]);
2038 err = E_DATA;
2039 goto bailout;
2040 }
2041 }
2042 } else {
2043 sys->ylist = ylist;
2044 ylist = NULL;
2045 }
2046
2047 /* Now compose a big list of all the variables in the system
2048 (referenced in structural equations or identities, or specified
2049 as instruments).
2050 */
2051 err = sys_make_biglist(sys);
2052 if (err) {
2053 goto bailout;
2054 }
2055
2056 /* If the user gave an endogenous list, check that it doesn't
2057 contain any variables that are not actually present in the
2058 system.
2059 */
2060 if (user_ylist) {
2061 for (j=1; j<=sys->ylist[0]; j++) {
2062 vj = sys->ylist[j];
2063 if (!in_gretl_list(sys->biglist, vj)) {
2064 gretl_errmsg_sprintf("%s is marked as endogenous but is "
2065 "not present in the system",
2066 dset->varname[vj]);
2067 err = E_DATA;
2068 goto bailout;
2069 }
2070 }
2071 }
2072
2073 /* By elimination from biglist, using ylist, compose a maximal
2074 list of exogenous and predetermined variables, "xplist".
2075 */
2076
2077 xplist = gretl_list_new(sys->biglist[0]);
2078 if (xplist == NULL) {
2079 err = E_ALLOC;
2080 goto bailout;
2081 }
2082
2083 k = 0;
2084 for (j=1; j<=sys->biglist[0]; j++) {
2085 vj = sys->biglist[j];
2086 if (!in_gretl_list(sys->ylist, vj)) {
2087 xplist[++k] = vj;
2088 }
2089 }
2090 xplist[0] = k;
2091
2092 #if SYSDEBUG
2093 printlist(xplist, "system auto xplist (exog + predet)");
2094 #endif
2095
2096 /* If the user gave a list of instruments (sys->ilist), check
2097 that it does not contain anything that is in fact clearly
2098 endogenous; otherwise use the xplist we computed above.
2099 */
2100 if (user_ilist) {
2101 for (j=1; j<=sys->ilist[0]; j++) {
2102 vj = sys->ilist[j];
2103 if (in_gretl_list(sys->ylist, vj)) {
2104 gretl_errmsg_sprintf("%s is marked as an instrument "
2105 "but is endogenous",
2106 dset->varname[vj]);
2107 err = E_DATA;
2108 goto bailout;
2109 }
2110 }
2111 } else {
2112 sys->ilist = gretl_list_copy(xplist);
2113 if (sys->ilist == NULL) {
2114 err = E_ALLOC;
2115 }
2116 }
2117
2118 if (!err && sys->ylist[0] != nlhs) {
2119 /* Note: check added 2009-08-10, modified 2012-04-05,
2120 moved in front of the following "xplist" block
2121 2015-10-29.
2122 */
2123 if (tsls_style && sys->ylist[0] > nlhs) {
2124 /* from the pov of the structural form, endogenous regressors
2125 without an equation should be treated "as if" exogenous?
2126 */
2127 tsls_style_shift_vars(sys, xplist);
2128 } else {
2129 gretl_errmsg_sprintf("Found %d endogenous variables but %d equations",
2130 sys->ylist[0], nlhs);
2131 err = E_DATA;
2132 }
2133 }
2134
2135 /* Now narrow down "xplist", removing variables that are in fact
2136 lags of endogenous variables (and recording them as such): this
2137 should leave a list of truly exogenous variables.
2138 */
2139 for (j=1; j<=xplist[0] && !err; j++) {
2140 vj = xplist[j];
2141 lag = series_get_lag(dset, vj);
2142 if (lag > 0) {
2143 vname = series_get_parent_name(dset, vj);
2144 if (vname != NULL) {
2145 src = sys_get_lag_src(vname, dset);
2146 if (in_gretl_list(sys->ylist, src)) {
2147 err = add_predet_to_sys(sys, dset, vj, src, lag);
2148 gretl_list_delete_at_pos(xplist, j--);
2149 }
2150 }
2151 }
2152 }
2153
2154 if (!err) {
2155 sys_xlist_reshuffle_const(xplist);
2156 }
2157
2158 #if SYSDEBUG
2159 printlist(xplist, "final system exog list");
2160 #endif
2161
2162 if (!err) {
2163 free(sys->xlist);
2164 sys->xlist = xplist;
2165 xplist = NULL;
2166 }
2167
2168 bailout:
2169
2170 free(ylist);
2171 free(xplist);
2172
2173 return err;
2174 }
2175
sys_has_user_name(equation_system * sys)2176 static int sys_has_user_name (equation_system *sys)
2177 {
2178 return (sys->name != NULL && *sys->name != '\0' &&
2179 !sys_anonymous(sys->name));
2180 }
2181
2182 #define ALLOW_SINGLE 1
2183
2184 /**
2185 * equation_system_finalize:
2186 * @sys: pre-defined equation system.
2187 * @dset: dataset struct.
2188 * @opt: may include OPT_V for verbose operation, OPT_S
2189 * to permit estimation of a single equation, OPT_R for
2190 * a robust variance matrix (OLS only), OPT_N for no
2191 * df correction (OLS and TSLS only).
2192 * @prn: printing struct.
2193 *
2194 * Finalize an equation system, e.g. in response to "end system".
2195 * If the system has a specified name, we save it on a stack of
2196 * defined systems. If it has a specified estimation method,
2197 * we go ahead and estimate it. If it has both, we do both.
2198 *
2199 * Returns: 0 on success, non-zero on error. If the system is
2200 * mal-formed, it is destroyed.
2201 */
2202
equation_system_finalize(equation_system * sys,DATASET * dset,gretlopt opt,PRN * prn)2203 int equation_system_finalize (equation_system *sys, DATASET *dset,
2204 gretlopt opt, PRN *prn)
2205 {
2206 #if ALLOW_SINGLE
2207 int mineq = 1;
2208 #else
2209 int mineq = (opt & OPT_S)? 1 : 2;
2210 #endif
2211 int err = 0;
2212
2213 #if SYSDEBUG
2214 fprintf(stderr, "*** equation_system_finalize\n");
2215 #endif
2216
2217 gretl_error_clear();
2218
2219 if (sys == NULL) {
2220 gretl_errmsg_set(_(nosystem));
2221 return 1;
2222 }
2223
2224 if (sys->neqns < mineq) {
2225 gretl_errmsg_set(_(toofew));
2226 equation_system_destroy(sys);
2227 return 1;
2228 }
2229
2230 if (sys->method >= SYS_METHOD_MAX) {
2231 gretl_errmsg_set(_(badsystem));
2232 equation_system_destroy(sys);
2233 return 1;
2234 }
2235
2236 err = sys_check_lists(sys, dset);
2237
2238 if (!err && !(opt & OPT_S) && sys_has_user_name(sys)) {
2239 /* save the system for subsequent estimation: but note that we
2240 should not do this if given OPT_S, for single-equation
2241 LIML
2242 */
2243 err = gretl_stack_object_as(sys, GRETL_OBJ_SYS, sys->name);
2244 }
2245
2246 if (!err && sys->method >= 0) {
2247 if (sys->flags & SYSTEM_QUIET) {
2248 opt |= OPT_Q;
2249 }
2250 err = equation_system_estimate(sys, dset, opt, prn);
2251 }
2252
2253 return err;
2254 }
2255
2256 /* Implement the "estimate" command. The full form of the
2257 command line is
2258
2259 estimate <sysname> method=<method>
2260
2261 where <sysname> is the name of a previously defined
2262 equation system and <method> is a recognized estimation
2263 method.
2264
2265 However, if the "last model" is an equation system, we
2266 can accept just
2267
2268 estimate method=<method>
2269
2270 where the system is not named but is implicitly the
2271 last model.
2272
2273 The method=<method> field is optional, provided that the
2274 (named or last-model) system already has a specified
2275 estimation method; so the minimal form of the command
2276 line is simply "estimate".
2277 */
2278
estimate_named_system(const char * sysname,const char * param,DATASET * dset,gretlopt opt,PRN * prn)2279 int estimate_named_system (const char *sysname, const char *param,
2280 DATASET *dset, gretlopt opt, PRN *prn)
2281 {
2282 equation_system *sys = NULL;
2283 int err = 0;
2284
2285 #if SYSDEBUG
2286 fprintf(stderr, "*** estimate_named_system: '%s'\n", sysname);
2287 #endif
2288
2289 if (sysname != NULL) {
2290 /* we got a name */
2291 if (sys_anonymous(sysname)) {
2292 sys = get_anonymous_equation_system();
2293 } else {
2294 sys = get_equation_system_by_name(sysname);
2295 }
2296 if (sys == NULL) {
2297 gretl_errmsg_sprintf(_("'%s': unrecognized name"), sysname);
2298 err = E_DATA;
2299 }
2300 } else {
2301 /* no name given: try "last model"? */
2302 GretlObjType type;
2303 void *ptr;
2304
2305 ptr = get_last_model(&type);
2306 if (ptr == NULL || type != GRETL_OBJ_SYS) {
2307 gretl_errmsg_sprintf(_("%s: no system was specified"), "estimate");
2308 err = E_DATA;
2309 } else {
2310 sys = ptr;
2311 }
2312 }
2313
2314 #if 1 /* do we really want this? */
2315 if (err) {
2316 /* we haven't found a system to estimate yet */
2317 sys = get_anonymous_equation_system();
2318 if (sys != NULL) {
2319 gretl_error_clear();
2320 err = 0;
2321 }
2322 }
2323 #endif
2324
2325 if (!err) {
2326 int method;
2327
2328 if (param != NULL) {
2329 method = sys_get_estimator(param);
2330 } else {
2331 method = sys->method;
2332 }
2333
2334 if (method < 0 || method >= SYS_METHOD_MAX) {
2335 gretl_errmsg_set("estimate: no valid method was specified");
2336 err = E_DATA;
2337 } else {
2338 sys->method = method;
2339 err = equation_system_estimate(sys, dset, opt, prn);
2340 }
2341 }
2342
2343 #if SYSDEBUG
2344 fprintf(stderr, "*** estimate_named_system: returning %d\n", err);
2345 #endif
2346
2347 return err;
2348 }
2349
2350 /* effective list length, allowance made for IVREG-style lists */
2351
get_real_list_length(const int * list)2352 static int get_real_list_length (const int *list)
2353 {
2354 int i, len = list[0];
2355
2356 for (i=1; i<=list[0]; i++) {
2357 if (list[i] == LISTSEP) {
2358 len = i - 1;
2359 break;
2360 }
2361 }
2362
2363 return len;
2364 }
2365
2366 /* below: used in formulating restriction matrices for
2367 an equation system */
2368
system_get_list_length(const equation_system * sys,int i)2369 int system_get_list_length (const equation_system *sys, int i)
2370 {
2371 if (i >= 0 && i < sys->neqns) {
2372 return get_real_list_length(sys->lists[i]);
2373 } else {
2374 return 0;
2375 }
2376 }
2377
system_max_indep_vars(const equation_system * sys)2378 int system_max_indep_vars (const equation_system *sys)
2379 {
2380 int i, nvi, nv = 0;
2381
2382 for (i=0; i<sys->neqns; i++) {
2383 nvi = get_real_list_length(sys->lists[i]) - 1;
2384 if (nvi > nv) nv = nvi;
2385 }
2386
2387 return nv;
2388 }
2389
system_n_indep_vars(const equation_system * sys)2390 int system_n_indep_vars (const equation_system *sys)
2391 {
2392 int i, nvi, nv = 0;
2393
2394 for (i=0; i<sys->neqns; i++) {
2395 nvi = get_real_list_length(sys->lists[i]) - 1;
2396 nv += nvi;
2397 }
2398
2399 return nv;
2400 }
2401
system_short_string(const MODEL * pmod)2402 const char *system_short_string (const MODEL *pmod)
2403 {
2404 int i = gretl_model_get_int(pmod, "method");
2405
2406 return gretl_system_short_strings[i];
2407 }
2408
equation_system_set_name(equation_system * sys,const char * name)2409 void equation_system_set_name (equation_system *sys, const char *name)
2410 {
2411 if (name == sys->name) {
2412 return;
2413 }
2414
2415 if (sys->name == NULL) {
2416 sys->name = malloc(MAXSAVENAME);
2417 }
2418
2419 if (sys->name != NULL) {
2420 *sys->name = '\0';
2421 strncat(sys->name, name, MAXSAVENAME - 1);
2422 }
2423 }
2424
compose_ivreg_list(const equation_system * sys,int i)2425 int *compose_ivreg_list (const equation_system *sys, int i)
2426 {
2427 int *list;
2428 int j, k1, k2;
2429
2430 if (i >= sys->neqns) {
2431 return NULL;
2432 }
2433
2434 k1 = sys->lists[i][0];
2435 k2 = sys->ilist[0];
2436
2437 list = gretl_list_new(k1 + k2 + 1);
2438 if (list == NULL) {
2439 return NULL;
2440 }
2441
2442 for (j=1; j<=list[0]; j++) {
2443 if (j <= k1) {
2444 list[j] = sys->lists[i][j];
2445 } else if (j == k1 + 1) {
2446 list[j] = LISTSEP;
2447 } else {
2448 list[j] = sys->ilist[j - (k1 + 1)];
2449 }
2450 }
2451
2452 return list;
2453 }
2454
system_normality_test(const equation_system * sys,gretlopt opt,PRN * prn)2455 int system_normality_test (const equation_system *sys,
2456 gretlopt opt, PRN *prn)
2457 {
2458 int err = 0;
2459
2460 if (sys->E == NULL || sys->S == NULL) {
2461 err = 1;
2462 } else {
2463 err = multivariate_normality_test(sys->E,
2464 sys->S,
2465 opt,
2466 prn);
2467 }
2468
2469 return err;
2470 }
2471
system_get_resid_series(equation_system * sys,int eqnum,DATASET * dset,int * err)2472 double *system_get_resid_series (equation_system *sys, int eqnum,
2473 DATASET *dset, int *err)
2474 {
2475 double *u = NULL;
2476 int t;
2477
2478 if (sys->E == NULL) {
2479 *err = E_DATA;
2480 return NULL;
2481 }
2482
2483 u = malloc(dset->n * sizeof *u);
2484 if (u == NULL) {
2485 *err = E_ALLOC;
2486 return NULL;
2487 }
2488
2489 for (t=0; t<dset->n; t++) {
2490 if (t < sys->t1 || t > sys->t2) {
2491 u[t] = NADBL;
2492 } else {
2493 u[t] = gretl_matrix_get(sys->E, t - sys->t1, eqnum);
2494 }
2495 }
2496
2497 return u;
2498 }
2499
system_get_full_string(const equation_system * sys)2500 static gchar *system_get_full_string (const equation_system *sys)
2501 {
2502 const char *lstr = gretl_system_long_strings[sys->method];
2503 gchar *ret;
2504
2505 if (sys->flags & SYSTEM_ITERATE) {
2506 ret = g_strdup_printf(_("iterated %s"), _(lstr));
2507 } else {
2508 ret = g_strdup(_(lstr));
2509 }
2510
2511 return ret;
2512 }
2513
2514 /* simple accessor functions */
2515
system_want_df_corr(const equation_system * sys)2516 int system_want_df_corr (const equation_system *sys)
2517 {
2518 return sys->flags & SYSTEM_DFCORR;
2519 }
2520
system_n_restrictions(const equation_system * sys)2521 int system_n_restrictions (const equation_system *sys)
2522 {
2523 int nr = 0;
2524
2525 if (sys->R != NULL && (sys->flags & SYSTEM_RESTRICT)) {
2526 nr = gretl_matrix_rows(sys->R);
2527 }
2528
2529 return nr;
2530 }
2531
system_get_list(const equation_system * sys,int i)2532 int *system_get_list (const equation_system *sys, int i)
2533 {
2534 if (i >= sys->neqns) return NULL;
2535
2536 return sys->lists[i];
2537 }
2538
system_get_depvar(const equation_system * sys,int i)2539 int system_get_depvar (const equation_system *sys, int i)
2540 {
2541 if (i >= sys->neqns) return 0;
2542
2543 return sys->lists[i][1];
2544 }
2545
system_get_method(const equation_system * sys)2546 int system_get_method (const equation_system *sys)
2547 {
2548 return sys->method;
2549 }
2550
system_get_endog_vars(const equation_system * sys)2551 int *system_get_endog_vars (const equation_system *sys)
2552 {
2553 return sys->ylist;
2554 }
2555
system_get_instr_vars(const equation_system * sys)2556 int *system_get_instr_vars (const equation_system *sys)
2557 {
2558 return sys->ilist;
2559 }
2560
system_attach_coeffs(equation_system * sys,gretl_matrix * b)2561 void system_attach_coeffs (equation_system *sys, gretl_matrix *b)
2562 {
2563 gretl_matrix_replace(&sys->b, b);
2564 }
2565
system_attach_vcv(equation_system * sys,gretl_matrix * vcv)2566 void system_attach_vcv (equation_system *sys, gretl_matrix *vcv)
2567 {
2568 gretl_matrix_replace(&sys->vcv, vcv);
2569 }
2570
system_attach_sigma(equation_system * sys,gretl_matrix * S)2571 void system_attach_sigma (equation_system *sys, gretl_matrix *S)
2572 {
2573 gretl_matrix_replace(&sys->S, S);
2574 }
2575
system_attach_uhat(equation_system * sys,gretl_matrix * E)2576 void system_attach_uhat (equation_system *sys, gretl_matrix *E)
2577 {
2578 gretl_matrix_replace(&sys->E, E);
2579 }
2580
system_get_model(const equation_system * sys,int i)2581 MODEL *system_get_model (const equation_system *sys, int i)
2582 {
2583 if (sys->models == NULL || i >= sys->neqns) {
2584 return NULL;
2585 }
2586
2587 return sys->models[i];
2588 }
2589
2590 /* for case of applying df correction to cross-equation
2591 residual covariance matrix */
2592
system_vcv_geomean(const equation_system * sys)2593 int system_vcv_geomean (const equation_system *sys)
2594 {
2595 return sys->flags & SYSTEM_VCV_GEOMEAN;
2596 }
2597
sys_eqn_indep_coeffs(const equation_system * sys,int eq)2598 static int sys_eqn_indep_coeffs (const equation_system *sys, int eq)
2599 {
2600 int nc;
2601
2602 if (eq >= sys->neqns) {
2603 nc = 0;
2604 } else if (sys->R == NULL) {
2605 nc = sys->lists[eq][0] - 1;
2606 } else {
2607 int nr = gretl_matrix_rows(sys->R);
2608 int cols = gretl_matrix_cols(sys->R);
2609 int cmax, cmin = 0;
2610 int single, i, j;
2611
2612 nc = sys->lists[eq][0] - 1;
2613
2614 for (i=0; i<eq; i++) {
2615 cmin += sys->lists[i][0] - 1;
2616 }
2617
2618 cmax = cmin + nc;
2619
2620 for (i=0; i<nr; i++) {
2621 single = 1;
2622 for (j=0; j<cols; j++) {
2623 if (j >= cmin && j < cmax) {
2624 if (gretl_matrix_get(sys->R, i, j) != 0.0) {
2625 single = 0;
2626 break;
2627 }
2628 }
2629 }
2630 if (single) {
2631 nc--;
2632 }
2633 }
2634 }
2635
2636 return nc;
2637 }
2638
2639 double
system_vcv_denom(const equation_system * sys,int i,int j)2640 system_vcv_denom (const equation_system *sys, int i, int j)
2641 {
2642 double den = sys->T;
2643
2644 if ((sys->flags & SYSTEM_VCV_GEOMEAN) &&
2645 i < sys->neqns && j < sys->neqns) {
2646 int ki = sys_eqn_indep_coeffs(sys, i);
2647
2648 if (j == i) {
2649 den = sys->T - ki;
2650 } else {
2651 int kj = sys_eqn_indep_coeffs(sys, j);
2652
2653 den = (sys->T - ki) * (sys->T - kj);
2654 den = sqrt(den);
2655 }
2656 }
2657
2658 return den;
2659 }
2660
2661 /* for system over-identification test */
2662
system_get_overid_df(const equation_system * sys)2663 int system_get_overid_df (const equation_system *sys)
2664 {
2665 int gl, i, k = 0;
2666
2667 gl = sys->neqns * sys->ilist[0];
2668
2669 for (i=0; i<sys->neqns; i++) {
2670 k += sys->lists[i][0] - 1;
2671 }
2672
2673 return gl - k;
2674 }
2675
2676 /* dealing with identities (FIML, LIML)
2677 (really also needed for LIML? not clear why...) */
2678
rhs_var_in_identity(const equation_system * sys,int lhsvar,int rhsvar)2679 int rhs_var_in_identity (const equation_system *sys, int lhsvar,
2680 int rhsvar)
2681 {
2682 const identity *ident;
2683 int i, j;
2684
2685 for (i=0; i<sys->nidents; i++) {
2686 ident = sys->idents[i];
2687 if (ident->depvar == lhsvar) {
2688 for (j=0; j<ident->n_atoms; j++) {
2689 if (ident->atoms[j].varnum == rhsvar) {
2690 return (ident->atoms[j].op == OP_PLUS)? 1 : -1;
2691 }
2692 }
2693 }
2694 }
2695
2696 return 0;
2697 }
2698
destroy_ident(identity * ident)2699 static void destroy_ident (identity *ident)
2700 {
2701 free(ident->atoms);
2702 free(ident);
2703 }
2704
ident_new(int nv)2705 static identity *ident_new (int nv)
2706 {
2707 identity *ident;
2708
2709 ident = malloc(sizeof *ident);
2710 if (ident == NULL) return NULL;
2711
2712 ident->n_atoms = nv;
2713 ident->atoms = malloc(nv * sizeof *ident->atoms);
2714 if (ident->atoms == NULL) {
2715 free(ident);
2716 ident = NULL;
2717 }
2718
2719 return ident;
2720 }
2721
2722 static identity *
parse_identity(const char * str,DATASET * dset,int * err)2723 parse_identity (const char *str, DATASET *dset, int *err)
2724 {
2725 identity *ident;
2726 char *test;
2727 const char *p = str;
2728 char vname[VNAMELEN];
2729 int lag, gotop = 0;
2730 int i, v, nv, len;
2731
2732 #if SYSDEBUG
2733 fprintf(stderr, "parse_identity: input = '%s'\n", str);
2734 #endif
2735
2736 /* First pass: count the elements and check everything for validity.
2737 nv is the number of variables appearing on the right-hand side
2738 of the identity (nv >= 1).
2739 */
2740
2741 i = 0;
2742 while (*p && !*err) {
2743 p += strspn(p, " ");
2744 if (i == 0) {
2745 /* left-hand side variable */
2746 *err = extract_varname(vname, p, &len);
2747 if (!*err) {
2748 v = current_series_index(dset, vname);
2749 if (v < 0) {
2750 *err = E_UNKVAR;
2751 } else {
2752 p += len;
2753 p += strspn(p, " ");
2754 if (*p != '=') {
2755 *err = E_PARSE;
2756 } else {
2757 p++;
2758 i++;
2759 gotop = 1;
2760 }
2761 }
2762 }
2763 } else if (*p == '+' || *p == '-') {
2764 gotop = 1;
2765 p++;
2766 } else if (i > 1 && !gotop) {
2767 /* rhs: no + or - given */
2768 *err = E_PARSE;
2769 } else {
2770 /* right-hand size variable (may be lag) */
2771 *err = extract_varname(vname, p, &len);
2772 if (!*err && gotop && len == 0) {
2773 /* dangling operator */
2774 *err = E_PARSE;
2775 }
2776 if (!*err) {
2777 v = current_series_index(dset, vname);
2778 if (v < 0) {
2779 *err = E_UNKVAR;
2780 } else {
2781 p += len;
2782 if (*p == '(') {
2783 lag = strtol(p + 1, &test, 10);
2784 if (*test != ')') {
2785 *err = E_PARSE;
2786 } else if (lag >= 0) {
2787 *err = E_DATA;
2788 } else {
2789 v = laggenr(v, -lag, dset);
2790 if (v < 0) {
2791 *err = E_DATA;
2792 } else {
2793 p = test + 1;
2794 }
2795 }
2796 }
2797 i++;
2798 gotop = 0;
2799 }
2800 }
2801 }
2802 }
2803
2804 if (gotop) {
2805 /* trailing operator */
2806 *err = E_PARSE;
2807 }
2808
2809 if (!*err) {
2810 nv = i - 1;
2811 if (nv < 1) {
2812 *err = E_PARSE;
2813 } else {
2814 ident = ident_new(nv);
2815 if (ident == NULL) {
2816 *err = E_ALLOC;
2817 return NULL;
2818 }
2819 }
2820 }
2821
2822 if (*err) {
2823 return NULL;
2824 }
2825
2826 ident->atoms[0].op = OP_PLUS;
2827
2828 /* Second pass: fill out the identity */
2829
2830 p = str;
2831 i = 0;
2832 while (*p) {
2833 p += strspn(p, " ");
2834 if (i == 0) {
2835 extract_varname(vname, p, &len);
2836 ident->depvar = series_index(dset, vname);
2837 p = strchr(p, '=') + 1;
2838 i++;
2839 } else if (*p == '+' || *p == '-') {
2840 ident->atoms[i-1].op = (*p == '+')? OP_PLUS : OP_MINUS;
2841 p++;
2842 } else {
2843 extract_varname(vname, p, &len);
2844 v = series_index(dset, vname);
2845 p += len;
2846 if (*p == '(') {
2847 lag = strtol(p + 1, &test, 10);
2848 v = laggenr(v, -lag, dset);
2849 p = test + 1;
2850 }
2851 ident->atoms[i-1].varnum = v;
2852 i++;
2853 }
2854 }
2855
2856 return ident;
2857 }
2858
2859 /* add information regarding a predetermined regressor */
2860
2861 static int
add_predet_to_sys(equation_system * sys,const DATASET * dset,int id,int src,int lag)2862 add_predet_to_sys (equation_system *sys, const DATASET *dset,
2863 int id, int src, int lag)
2864 {
2865 int n = sys->plist[0];
2866 int *test;
2867 predet *pre;
2868
2869 if (id < 0 || src < 0 || id >= dset->v || src >= dset->v) {
2870 /* something screwy */
2871 return E_DATA;
2872 }
2873
2874 if (in_gretl_list(sys->plist, id)) {
2875 /* already present */
2876 return 0;
2877 }
2878
2879 pre = realloc(sys->pre_vars, (n + 1) * sizeof *pre);
2880 if (pre == NULL) {
2881 return E_ALLOC;
2882 }
2883
2884 sys->pre_vars = pre;
2885 sys->pre_vars[n].id = id;
2886 sys->pre_vars[n].src = src;
2887 sys->pre_vars[n].lag = lag;
2888
2889 /* add to list of predetermined variable IDs */
2890 test = gretl_list_append_term(&sys->plist, id);
2891 if (test == NULL) {
2892 return E_ALLOC;
2893 }
2894
2895 return 0;
2896 }
2897
2898 static int
add_identity_to_sys(equation_system * sys,const char * line,DATASET * dset)2899 add_identity_to_sys (equation_system *sys, const char *line,
2900 DATASET *dset)
2901 {
2902 identity **pident;
2903 identity *ident;
2904 int ni = sys->nidents;
2905 int err = 0;
2906
2907 ident = parse_identity(line, dset, &err);
2908 if (ident == NULL) return err;
2909
2910 /* connect the identity to the equation system */
2911 pident = realloc(sys->idents, (ni + 1) * sizeof *sys->idents);
2912 if (pident == NULL) {
2913 destroy_ident(ident);
2914 return E_ALLOC;
2915 }
2916
2917 sys->idents = pident;
2918 sys->idents[ni] = ident;
2919 sys->nidents += 1;
2920
2921 return 0;
2922 }
2923
2924 static int
add_aux_list_to_sys(equation_system * sys,const char * line,DATASET * dset,int which)2925 add_aux_list_to_sys (equation_system *sys, const char *line,
2926 DATASET *dset, int which)
2927 {
2928 int *list;
2929 int err = 0;
2930
2931 if (which == ENDOG_LIST) {
2932 if (sys->ylist != NULL) {
2933 gretl_errmsg_set("Only one list of endogenous variables may be given");
2934 return 1;
2935 }
2936 } else if (which == INSTR_LIST) {
2937 if (sys->ilist != NULL) {
2938 gretl_errmsg_set("Only one list of instruments may be given");
2939 return 1;
2940 }
2941 } else {
2942 return E_DATA;
2943 }
2944
2945 line += strspn(line, " ");
2946
2947 list = generate_list(line, dset, &err);
2948
2949 if (!err) {
2950 if (which == ENDOG_LIST) {
2951 sys->ylist = list;
2952 } else {
2953 sys->ilist = list;
2954 }
2955 }
2956
2957 return err;
2958 }
2959
2960 /**
2961 * system_parse_line:
2962 * @sys: initialized equation system.
2963 * @line: command line.
2964 * @dset: dataset struct.
2965 *
2966 * Modifies @sys according to the command supplied in @line,
2967 * which must start with "identity" (and supply an identity
2968 * to be added to the system, or "endog" (and supply a list
2969 * of endogenous variables), or "instr" (and supply a list of
2970 * instrumental variables).
2971 *
2972 * Returns: 0 on success, non-zero on failure, in which case
2973 * @sys is destroyed.
2974 */
2975
2976 int
system_parse_line(equation_system * sys,const char * line,DATASET * dset)2977 system_parse_line (equation_system *sys, const char *line,
2978 DATASET *dset)
2979 {
2980 int err = 0;
2981
2982 gretl_error_clear();
2983
2984 #if SYSDEBUG > 1
2985 fprintf(stderr, "*** system_parse_line: '%s'\n", line);
2986 #endif
2987
2988 line += strspn(line, " ");
2989
2990 if (!strncmp(line, "identity ", 9)) {
2991 err = add_identity_to_sys(sys, line + 8, dset);
2992 } else if (!strncmp(line, "endog ", 6)) {
2993 err = add_aux_list_to_sys(sys, line + 5, dset, ENDOG_LIST);
2994 } else if (!strncmp(line, "instr ", 6)) {
2995 err = add_aux_list_to_sys(sys, line + 5, dset, INSTR_LIST);
2996 } else {
2997 err = E_PARSE;
2998 }
2999
3000 if (err) {
3001 equation_system_destroy(sys);
3002 }
3003
3004 return err;
3005 }
3006
3007 void
system_set_restriction_matrices(equation_system * sys,gretl_matrix * R,gretl_matrix * q)3008 system_set_restriction_matrices (equation_system *sys,
3009 gretl_matrix *R, gretl_matrix *q)
3010 {
3011 system_clear_restrictions(sys);
3012
3013 sys->R = R;
3014 sys->q = q;
3015
3016 sys->flags |= SYSTEM_RESTRICT;
3017 }
3018
3019 double *
equation_system_get_series(const equation_system * sys,const DATASET * dset,int idx,const char * key,int * err)3020 equation_system_get_series (const equation_system *sys,
3021 const DATASET *dset,
3022 int idx, const char *key, int *err)
3023 {
3024 const gretl_matrix *M = NULL;
3025 double *x = NULL;
3026 const char *msel;
3027 int t, col = 0;
3028
3029 if (sys == NULL || (idx != M_UHAT && idx != M_YHAT)) {
3030 *err = E_BADSTAT;
3031 return NULL;
3032 }
3033
3034 msel = strchr(key, '[');
3035 if (msel == NULL || sscanf(msel, "[,%d]", &col) != 1) {
3036 *err = E_PARSE;
3037 } else if (col <= 0 || col > sys->neqns) {
3038 *err = E_DATA;
3039 }
3040
3041 if (!*err) {
3042 M = (idx == M_UHAT)? sys->E : sys->yhat;
3043 if (M == NULL) {
3044 *err = E_DATA;
3045 }
3046 }
3047
3048 if (!*err) {
3049 x = malloc(dset->n * sizeof *x);
3050 if (x == NULL) {
3051 *err = E_ALLOC;
3052 }
3053 }
3054
3055 if (!*err) {
3056 int s = 0;
3057
3058 col--; /* switch to 0-based */
3059 for (t=0; t<dset->n; t++) {
3060 if (t < sys->t1 || t > sys->t2) {
3061 x[t] = NADBL;
3062 } else {
3063 x[t] = gretl_matrix_get(M, s++, col);
3064 }
3065 }
3066 }
3067
3068 return x;
3069 }
3070
system_add_yhat_matrix(equation_system * sys)3071 static int system_add_yhat_matrix (equation_system *sys)
3072 {
3073 double x, avc;
3074 int k = 0;
3075 int i, s, t;
3076
3077 sys->yhat = gretl_matrix_alloc(sys->T, sys->neqns);
3078 if (sys->yhat == NULL) {
3079 return E_ALLOC;
3080 }
3081
3082 gretl_matrix_set_t1(sys->yhat, sys->t1);
3083 gretl_matrix_set_t2(sys->yhat, sys->t2);
3084
3085 for (i=0; i<sys->neqns; i++) {
3086 s = 0;
3087 for (t=sys->t1; t<=sys->t2; t++) {
3088 x = sys->models[i]->yhat[t];
3089 gretl_matrix_set(sys->yhat, s++, i, x);
3090 }
3091 k += sys->models[i]->ncoeff;
3092 }
3093
3094 k -= system_n_restrictions(sys);
3095
3096 avc = k / (double) sys->neqns;
3097 sys->df = sys->T - floor(avc);
3098
3099 return 0;
3100 }
3101
get_stderr_vec(const gretl_matrix * V)3102 static gretl_matrix *get_stderr_vec (const gretl_matrix *V)
3103 {
3104 gretl_matrix *se;
3105 int k = V->rows;
3106
3107 se = gretl_column_vector_alloc(k);
3108
3109 if (se != NULL) {
3110 double x;
3111 int i;
3112
3113 for (i=0; i<k; i++) {
3114 x = gretl_matrix_get(V, i, i);
3115 se->val[i] = sqrt(x);
3116 }
3117 }
3118
3119 return se;
3120 }
3121
3122 /* retrieve a copy of a specified matrix from an equation
3123 system */
3124
3125 gretl_matrix *
equation_system_get_matrix(const equation_system * sys,int idx,int * err)3126 equation_system_get_matrix (const equation_system *sys, int idx,
3127 int *err)
3128 {
3129 gretl_matrix *M = NULL;
3130
3131 if (sys == NULL) {
3132 *err = E_BADSTAT;
3133 return NULL;
3134 }
3135
3136 switch (idx) {
3137 case M_COEFF:
3138 if (sys->b == NULL) {
3139 *err = E_BADSTAT;
3140 } else {
3141 M = gretl_matrix_copy(sys->b);
3142 }
3143 break;
3144 case M_UHAT:
3145 M = gretl_matrix_copy(sys->E);
3146 break;
3147 case M_YHAT:
3148 M = gretl_matrix_copy(sys->yhat);
3149 break;
3150 case M_VCV:
3151 case M_SE:
3152 if (sys->vcv == NULL) {
3153 *err = E_BADSTAT;
3154 } else if (idx == M_SE) {
3155 M = get_stderr_vec(sys->vcv);
3156 } else {
3157 M = gretl_matrix_copy(sys->vcv);
3158 }
3159 break;
3160 case M_SIGMA:
3161 M = gretl_matrix_copy(sys->S);
3162 break;
3163 case M_SYSGAM:
3164 if (sys->Gamma == NULL) {
3165 *err = E_BADSTAT;
3166 } else {
3167 M = gretl_matrix_copy(sys->Gamma);
3168 }
3169 break;
3170 case M_SYSA:
3171 if (sys->A == NULL) {
3172 *err = E_BADSTAT;
3173 } else {
3174 M = gretl_matrix_copy(sys->A);
3175 }
3176 break;
3177 case M_SYSB:
3178 if (sys->B == NULL) {
3179 *err = E_BADSTAT;
3180 } else {
3181 M = gretl_matrix_copy(sys->B);
3182 }
3183 break;
3184 default:
3185 *err = E_BADSTAT;
3186 break;
3187 }
3188
3189 if (M == NULL && !*err) {
3190 *err = E_ALLOC;
3191 }
3192
3193 return M;
3194 }
3195
highest_numbered_var_in_system(const equation_system * sys,const DATASET * dset)3196 int highest_numbered_var_in_system (const equation_system *sys,
3197 const DATASET *dset)
3198 {
3199 int i, j, v, vmax = 0;
3200
3201 if (sys->biglist != NULL) {
3202 for (j=1; j<=sys->biglist[0]; j++) {
3203 v = sys->biglist[j];
3204 if (v > vmax) {
3205 vmax = v;
3206 }
3207 }
3208 } else {
3209 /* should not happen */
3210 for (i=0; i<sys->neqns; i++) {
3211 for (j=1; j<=sys->lists[i][0]; j++) {
3212 v = sys->lists[i][j];
3213 if (v == LISTSEP || v >= dset->v) {
3214 /* temporary variables, already gone? */
3215 continue;
3216 }
3217 if (v > vmax) {
3218 vmax = v;
3219 }
3220 }
3221 }
3222 }
3223
3224 return vmax;
3225 }
3226
sys_retrieve_identity(xmlNodePtr node,int * err)3227 static identity *sys_retrieve_identity (xmlNodePtr node, int *err)
3228 {
3229 identity *ident;
3230 xmlNodePtr cur;
3231 int n_atoms, depvar;
3232 int i, got = 0;
3233
3234 got += gretl_xml_get_prop_as_int(node, "n_atoms", &n_atoms);
3235 got += gretl_xml_get_prop_as_int(node, "depvar", &depvar);
3236 if (got < 2) {
3237 *err = E_DATA;
3238 return NULL;
3239 }
3240
3241 ident = ident_new(n_atoms);
3242 if (ident == NULL) {
3243 *err = E_ALLOC;
3244 return NULL;
3245 }
3246
3247 ident->depvar = depvar;
3248
3249 cur = node->xmlChildrenNode;
3250
3251 i = 0;
3252 while (cur != NULL && !*err) {
3253 if (!xmlStrcmp(cur->name, (XUC) "id_atom")) {
3254 got = gretl_xml_get_prop_as_int(cur, "op", &ident->atoms[i].op);
3255 got += gretl_xml_get_prop_as_int(cur, "varnum", &ident->atoms[i].varnum);
3256 if (got < 2) {
3257 *err = E_DATA;
3258 } else {
3259 i++;
3260 }
3261 }
3262 cur = cur->next;
3263 }
3264
3265 if (!*err && i != n_atoms) {
3266 *err = E_DATA;
3267 }
3268
3269 if (*err) {
3270 destroy_ident(ident);
3271 ident = NULL;
3272 }
3273
3274 return ident;
3275 }
3276
3277 equation_system *
equation_system_from_XML(xmlNodePtr node,xmlDocPtr doc,int * err)3278 equation_system_from_XML (xmlNodePtr node, xmlDocPtr doc, int *err)
3279 {
3280 equation_system *sys;
3281 xmlNodePtr cur;
3282 char *sname;
3283 int method = 0;
3284 int i, j, got = 0;
3285
3286 got += gretl_xml_get_prop_as_string(node, "name", &sname);
3287 got += gretl_xml_get_prop_as_int(node, "method", &method);
3288
3289 if (got < 2) {
3290 *err = E_DATA;
3291 return NULL;
3292 }
3293
3294 sys = equation_system_new(method, sname, err);
3295 if (*err) {
3296 return NULL;
3297 }
3298
3299 got = 0;
3300 got += gretl_xml_get_prop_as_int(node, "n_equations", &sys->neqns);
3301 got += gretl_xml_get_prop_as_int(node, "nidents", &sys->nidents);
3302 got += gretl_xml_get_prop_as_int(node, "flags", &sys->flags);
3303
3304 if (got < 3) {
3305 *err = E_DATA;
3306 goto bailout;
3307 }
3308
3309 gretl_xml_get_prop_as_int(node, "order", &sys->order);
3310
3311 sys->lists = malloc(sys->neqns * sizeof sys->lists);
3312 if (sys->lists == NULL) {
3313 *err = E_ALLOC;
3314 goto bailout;
3315 }
3316
3317 if (sys->nidents > 0) {
3318 sys->idents = malloc(sys->nidents * sizeof sys->idents);
3319 if (sys->idents == NULL) {
3320 *err = E_ALLOC;
3321 goto bailout;
3322 }
3323 }
3324
3325 cur = node->xmlChildrenNode;
3326
3327 i = j = 0;
3328 while (cur != NULL && !*err) {
3329 if (!xmlStrcmp(cur->name, (XUC) "eqnlist")) {
3330 sys->lists[i++] = gretl_xml_get_list(cur, doc, err);
3331 } else if (!xmlStrcmp(cur->name, (XUC) "endog_vars")) {
3332 sys->ylist = gretl_xml_get_list(cur, doc, err);
3333 } else if (!xmlStrcmp(cur->name, (XUC) "instr_vars")) {
3334 sys->ilist = gretl_xml_get_list(cur, doc, err);
3335 } else if (!xmlStrcmp(cur->name, (XUC) "identity")) {
3336 sys->idents[j++] = sys_retrieve_identity(cur, err);
3337 } else if (!xmlStrcmp(cur->name, (XUC) "R")) {
3338 sys->R = gretl_xml_get_matrix(cur, doc, err);
3339 } else if (!xmlStrcmp(cur->name, (XUC) "q")) {
3340 sys->q = gretl_xml_get_matrix(cur, doc, err);
3341 }
3342 cur = cur->next;
3343 }
3344
3345 if (!*err && (i != sys->neqns || j != sys->nidents)) {
3346 *err = E_DATA;
3347 }
3348
3349 if (*err) {
3350 equation_system_destroy(sys);
3351 sys = NULL;
3352 }
3353
3354 bailout:
3355
3356 return sys;
3357 }
3358
xml_print_identity(identity * ident,PRN * prn)3359 static void xml_print_identity (identity *ident, PRN *prn)
3360 {
3361 int i;
3362
3363 pprintf(prn, "<identity n_atoms=\"%d\" depvar=\"%d\">\n",
3364 ident->n_atoms, ident->depvar);
3365 for (i=0; i<ident->n_atoms; i++) {
3366 pprintf(prn, " <id_atom op=\"%d\" varnum=\"%d\"/>\n",
3367 ident->atoms[i].op, ident->atoms[i].varnum);
3368 }
3369 pputs(prn, "</identity>\n");
3370 }
3371
equation_system_serialize(equation_system * sys,SavedObjectFlags flags,PRN * prn)3372 int equation_system_serialize (equation_system *sys,
3373 SavedObjectFlags flags,
3374 PRN *prn)
3375 {
3376 char *xmlname = NULL;
3377 int tsls_style = 0;
3378 int i, err = 0;
3379
3380 if (sys->name == NULL || *sys->name == '\0') {
3381 xmlname = gretl_strdup("none");
3382 } else {
3383 xmlname = gretl_xml_encode(sys->name);
3384 }
3385
3386 pprintf(prn, "<gretl-equation-system name=\"%s\" saveflags=\"%d\" method=\"%d\" ",
3387 xmlname, flags, sys->method);
3388 free(xmlname);
3389
3390 pprintf(prn, "n_equations=\"%d\" nidents=\"%d\" flags=\"%d\" order=\"%d\">\n",
3391 sys->neqns, sys->nidents, sys->flags, sys->order);
3392
3393 for (i=0; i<sys->neqns; i++) {
3394 gretl_xml_put_tagged_list("eqnlist", sys->lists[i], prn);
3395 }
3396
3397 for (i=0; i<sys->neqns; i++) {
3398 if (gretl_list_has_separator(sys->lists[i])) {
3399 tsls_style = 1;
3400 break;
3401 }
3402 }
3403
3404 if (!tsls_style) {
3405 gretl_xml_put_tagged_list("endog_vars", sys->ylist, prn);
3406 gretl_xml_put_tagged_list("instr_vars", sys->ilist, prn);
3407 }
3408
3409 for (i=0; i<sys->nidents; i++) {
3410 xml_print_identity(sys->idents[i], prn);
3411 }
3412
3413 gretl_matrix_serialize(sys->R, "R", prn);
3414 gretl_matrix_serialize(sys->q, "q", prn);
3415
3416 pputs(prn, "</gretl-equation-system>\n");
3417
3418 return err;
3419 }
3420
bundlize_diag_test(const equation_system * sys)3421 static gretl_bundle *bundlize_diag_test (const equation_system *sys)
3422 {
3423 gretl_bundle *b = NULL;
3424 double test = NADBL;
3425 double pval = NADBL;
3426 int k = sys->S->rows;
3427 int df = k * (k - 1) / 2;
3428
3429 if (sys->method == SYS_METHOD_SUR && sys->iters > 0) {
3430 if (!na(sys->ldet) && sys->diag_test > 0.0) {
3431 test = sys->T * (sys->diag_test - sys->ldet);
3432 pval = chisq_cdf_comp(df, test);
3433
3434 }
3435 } else {
3436 test = sys->diag_test;
3437 pval = chisq_cdf_comp(df, test);
3438 }
3439
3440 if (!na(test) && !na(pval)) {
3441 b = gretl_bundle_new();
3442 if (b != NULL) {
3443 gretl_bundle_set_scalar(b, "test", test);
3444 gretl_bundle_set_scalar(b, "pvalue", pval);
3445 gretl_bundle_set_int(b, "df", df);
3446 }
3447 }
3448
3449 return b;
3450 }
3451
equation_system_bundlize(equation_system * sys,gretl_bundle * b)3452 int equation_system_bundlize (equation_system *sys,
3453 gretl_bundle *b)
3454 {
3455 const char *s;
3456 char lname[16];
3457 int tsls_style = 0;
3458 int i, err = 0;
3459
3460 gretl_bundle_set_string(b, "command", "system");
3461
3462 gretl_bundle_set_int(b, "neqns", sys->neqns);
3463 s = system_method_short_string(sys->method);
3464 if (s != NULL) {
3465 gretl_bundle_set_string(b, "method", s);
3466 }
3467 if (sys->flags & SYSTEM_ROBUST) {
3468 gretl_bundle_set_int(b, "robust", 1);
3469 }
3470
3471 /* convert to 1-based for user-space */
3472 gretl_bundle_set_int(b, "t1", sys->t1 + 1);
3473 gretl_bundle_set_int(b, "t2", sys->t2 + 1);
3474
3475 gretl_bundle_set_matrix(b, "coeff", sys->b);
3476 gretl_bundle_set_matrix(b, "vcv", sys->vcv);
3477 gretl_bundle_set_matrix(b, "sigma", sys->S);
3478 gretl_bundle_set_matrix(b, "uhat", sys->E);
3479 gretl_bundle_set_matrix(b, "yhat", sys->yhat);
3480
3481 if (sys->Gamma != NULL) {
3482 gretl_bundle_set_matrix(b, "Gamma", sys->Gamma);
3483 }
3484 if (sys->A != NULL) {
3485 gretl_bundle_set_matrix(b, "A", sys->A);
3486 }
3487 if (sys->B != NULL) {
3488 gretl_bundle_set_matrix(b, "B", sys->B);
3489 }
3490
3491 /* per-equation model structs? */
3492
3493 for (i=0; i<sys->neqns; i++) {
3494 sprintf(lname, "list%d", i+1);
3495 gretl_bundle_set_list(b, lname, sys->lists[i]);
3496 }
3497
3498 for (i=0; i<sys->neqns; i++) {
3499 if (gretl_list_has_separator(sys->lists[i])) {
3500 tsls_style = 1;
3501 break;
3502 }
3503 }
3504
3505 if (!tsls_style) {
3506 gretl_bundle_set_list(b, "endog_vars", sys->ylist);
3507 gretl_bundle_set_list(b, "instr_vars", sys->ilist);
3508 }
3509
3510 if (sys->xlist != NULL) {
3511 gretl_bundle_set_list(b, "xlist", sys->xlist);
3512 }
3513
3514 if (sys->plist != NULL) {
3515 gretl_bundle_set_list(b, "predet_vars", sys->plist);
3516 }
3517
3518 #if 0 /* FIXME */
3519 for (i=0; i<sys->nidents; i++) {
3520 xml_print_identity(sys->idents[i], fp);
3521 }
3522 #endif
3523
3524 /* restrictions? */
3525 if (sys->R != NULL) {
3526 gretl_bundle_set_matrix(b, "R", sys->R);
3527 }
3528 if (sys->q != NULL) {
3529 gretl_bundle_set_matrix(b, "q", sys->q);
3530 }
3531
3532 if (sys->diag_test > 0) {
3533 gretl_bundle *dt = bundlize_diag_test(sys);
3534
3535 if (dt != NULL) {
3536 gretl_bundle_donate_data(b, "diag_test", dt,
3537 GRETL_TYPE_BUNDLE, 0);
3538 }
3539 }
3540
3541 return err;
3542 }
3543
sur_ols_diag(equation_system * sys)3544 static int sur_ols_diag (equation_system *sys)
3545 {
3546 double s2, ls2sum = 0.0;
3547 int i, err = 0;
3548
3549 for (i=0; i<sys->neqns; i++) {
3550 s2 = gretl_model_get_double(sys->models[i], "ols_sigma_squared");
3551 if (na(s2)) {
3552 err = 1;
3553 break;
3554 }
3555 ls2sum += log(s2);
3556 }
3557
3558 if (!err) {
3559 sys->diag_test = ls2sum;
3560 }
3561
3562 return err;
3563 }
3564
3565 static void
print_system_overid_test(const equation_system * sys,PRN * prn)3566 print_system_overid_test (const equation_system *sys, PRN *prn)
3567 {
3568 int tex = tex_format(prn);
3569 int df = system_get_overid_df(sys);
3570 double pv;
3571
3572 if (sys->method == SYS_METHOD_FIML && df > 0) {
3573 double X2;
3574
3575 if (na(sys->ll) || na(sys->llu) ||
3576 sys->ll == 0.0 || sys->llu == 0.0) {
3577 return;
3578 }
3579
3580 X2 = 2.0 * (sys->llu - sys->ll);
3581 pv = chisq_cdf_comp(df, X2);
3582
3583 if (tex) {
3584 pprintf(prn, "%s:\\\\\n", _("LR over-identification test"));
3585 if (sys->ll < 0) {
3586 pprintf(prn, " %s = $-$%g", _("Restricted log-likelihood"), -sys->ll);
3587 } else {
3588 pprintf(prn, " %s = %g", _("Restricted log-likelihood"), sys->ll);
3589 }
3590 gretl_prn_newline(prn);
3591 if (sys->llu < 0) {
3592 pprintf(prn, " %s = $-$%g", _("Unrestricted log-likelihood"), -sys->llu);
3593 } else {
3594 pprintf(prn, " %s = %g", _("Unrestricted log-likelihood"), sys->llu);
3595 }
3596 gretl_prn_newline(prn);
3597 pprintf(prn, " $\\chi^2(%d)$ = %g [%.4f]\n", df, X2, pv);
3598 } else {
3599 pprintf(prn, "%s:\n", _("LR over-identification test"));
3600 pprintf(prn, " %s = %g\n", _("Restricted log-likelihood"), sys->ll);
3601 pprintf(prn, " %s = %g\n", _("Unrestricted log-likelihood"), sys->llu);
3602 pprintf(prn, " %s(%d) = %g [%.4f]\n\n", _("Chi-square"), df, X2, pv);
3603 }
3604 } else if ((sys->method == SYS_METHOD_3SLS ||
3605 sys->method == SYS_METHOD_SUR) && df > 0) {
3606 if (na(sys->X2) || sys->X2 <= 0.0) {
3607 if (!tex) {
3608 pputs(prn, _("Warning: the Hansen-Sargan over-identification test "
3609 "failed.\nThis probably indicates that the estimation "
3610 "problem is ill-conditioned.\n"));
3611 pputc(prn, '\n');
3612 }
3613 return;
3614 }
3615
3616 pv = chisq_cdf_comp(df, sys->X2);
3617
3618 if (tex) {
3619 pprintf(prn, "\\noindent %s:\\\\\n",
3620 _("Hansen--Sargan over-identification test"));
3621 pprintf(prn, " $\\chi^2(%d)$ = %g [%.4f]\\\\\n", df, sys->X2, pv);
3622 } else {
3623 pprintf(prn, "%s:\n", _("Hansen-Sargan over-identification test"));
3624 pprintf(prn, " %s(%d) = %g [%.4f]\n\n", _("Chi-square"),
3625 df, sys->X2, pv);
3626 }
3627 }
3628 }
3629
system_diag_test(const equation_system * sys,double * test,double * pval)3630 int system_diag_test (const equation_system *sys, double *test,
3631 double *pval)
3632 {
3633 int k, df, err = 0;
3634
3635 if (sys->S == NULL) {
3636 return E_BADSTAT;
3637 }
3638
3639 k = sys->S->rows;
3640 df = k * (k - 1) / 2;
3641
3642 if (sys->method == SYS_METHOD_SUR && sys->iters > 0) {
3643 /* iterated SUR */
3644 if (!na(sys->ldet) && sys->diag_test != 0.0) {
3645 double lr = sys->T * (sys->diag_test - sys->ldet);
3646
3647 if (test != NULL) {
3648 *test = lr;
3649 }
3650 if (pval != NULL) {
3651 *pval = chisq_cdf_comp(df, lr);
3652 }
3653 } else {
3654 err = E_BADSTAT;
3655 }
3656 } else if (sys->diag_test > 0) {
3657 /* other estimators */
3658 if (test != NULL) {
3659 *test = sys->diag_test;
3660 }
3661 if (pval != NULL) {
3662 *pval = chisq_cdf_comp(df, sys->diag_test);
3663 }
3664 } else {
3665 err = E_BADSTAT;
3666 }
3667
3668 return err;
3669 }
3670
system_print_sigma(const equation_system * sys,PRN * prn)3671 int system_print_sigma (const equation_system *sys, PRN *prn)
3672 {
3673 double test, pval;
3674 int tex = tex_format(prn);
3675 int k, df;
3676
3677 if (sys->S == NULL) {
3678 return E_DATA;
3679 }
3680
3681 k = sys->S->rows;
3682 df = k * (k - 1) / 2;
3683
3684 print_contemp_covariance_matrix(sys->S, sys->ldet, prn);
3685
3686 if (sys->method == SYS_METHOD_SUR && sys->iters > 0) {
3687 if (!na(sys->ldet) && sys->diag_test != 0.0) {
3688 test = sys->T * (sys->diag_test - sys->ldet);
3689 pval = chisq_cdf_comp(df, test);
3690 if (tex) {
3691 pprintf(prn, "%s:\\\\\n", _("LR test for diagonal covariance matrix"));
3692 pprintf(prn, " $\\chi^2(%d)$ = %g [%.4f]", df, test, pval);
3693 gretl_prn_newline(prn);
3694 } else {
3695 pprintf(prn, "%s:\n", _("LR test for diagonal covariance matrix"));
3696 pprintf(prn, " %s(%d) = %g [%.4f]\n", _("Chi-square"),
3697 df, test, pval);
3698 }
3699 }
3700 } else {
3701 const char *labels[] = {
3702 N_("Breusch--Pagan test for diagonal covariance matrix"),
3703 N_("Robust Breusch--Pagan test for diagonal covariance matrix"),
3704 N_("Breusch-Pagan test for diagonal covariance matrix"),
3705 N_("Robust Breusch-Pagan test for diagonal covariance matrix"),
3706 };
3707 const char *label;
3708
3709 test = sys->diag_test;
3710 if (sys->flags & SYSTEM_ROBUST) {
3711 label = tex ? labels[1] : labels[3];
3712 } else {
3713 label = tex ? labels[0] : labels[2];
3714 }
3715 if (test > 0) {
3716 pval = chisq_cdf_comp(df, test);
3717 if (tex) {
3718 pprintf(prn, "%s:", _(label));
3719 gretl_prn_newline(prn);
3720 pprintf(prn, " $\\chi^2(%d)$ = %g [%.4f]", df, test, pval);
3721 gretl_prn_newline(prn);
3722 } else {
3723 pprintf(prn, "%s:\n", _(label));
3724 pprintf(prn, " %s(%d) = %g [%.4f]\n", _("Chi-square"),
3725 df, test, pval);
3726 }
3727 }
3728 }
3729
3730 pputc(prn, '\n');
3731
3732 return 0;
3733 }
3734
3735 enum {
3736 ENDOG,
3737 EXOG,
3738 PREDET
3739 };
3740
categorize_variable(int vnum,const equation_system * sys,int * col,int * lag)3741 static int categorize_variable (int vnum, const equation_system *sys,
3742 int *col, int *lag)
3743 {
3744 int pos, ret = -1;
3745
3746 *lag = 0;
3747
3748 pos = in_gretl_list(sys->ylist, vnum);
3749 if (pos > 0) {
3750 ret = ENDOG;
3751 } else {
3752 pos = in_gretl_list(sys->xlist, vnum);
3753 if (pos > 0) {
3754 ret = EXOG;
3755 } else {
3756 int v = get_predet_parent(sys, vnum, lag);
3757
3758 pos = in_gretl_list(sys->ylist, v);
3759 if (pos > 0) {
3760 ret = PREDET;
3761 }
3762 }
3763 }
3764
3765 if (pos > 0) {
3766 *col = pos - 1;
3767 }
3768
3769 return ret;
3770 }
3771
get_col_and_lag(int vnum,const equation_system * sys,int * col,int * lag)3772 static int get_col_and_lag (int vnum, const equation_system *sys,
3773 int *col, int *lag)
3774 {
3775 int vp, pos, err = 0;
3776
3777 vp = get_predet_parent(sys, vnum, lag);
3778 pos = in_gretl_list(sys->ylist, vp);
3779 if (pos > 0) {
3780 *col = pos - 1;
3781 } else {
3782 *col = 0;
3783 err = 1;
3784 }
3785
3786 return err;
3787 }
3788
3789 /* reduced-form error covariance matrix */
3790
3791 static int
sys_add_RF_covariance_matrix(equation_system * sys,int n)3792 sys_add_RF_covariance_matrix (equation_system *sys, int n)
3793 {
3794 gretl_matrix *G = NULL, *S = NULL;
3795 int err = 0;
3796
3797 if (!gretl_is_identity_matrix(sys->Gamma)) {
3798 G = gretl_matrix_copy(sys->Gamma);
3799 if (G == NULL) {
3800 return E_ALLOC;
3801 }
3802 }
3803
3804 if (n > sys->S->rows) {
3805 S = gretl_zero_matrix_new(n, n);
3806 if (S == NULL) {
3807 gretl_matrix_free(G);
3808 return E_ALLOC;
3809 }
3810 gretl_matrix_inscribe_matrix(S, sys->S, 0, 0, GRETL_MOD_NONE);
3811 }
3812
3813 sys->Sr = gretl_matrix_alloc(n, n);
3814 if (sys->Sr == NULL) {
3815 gretl_matrix_free(G);
3816 gretl_matrix_free(S);
3817 return E_ALLOC;
3818 }
3819
3820 if (G != NULL) {
3821 err = gretl_SVD_invert_matrix(G);
3822 if (!err) {
3823 err = gretl_matrix_qform(G, GRETL_MOD_NONE,
3824 (S != NULL)? S : sys->S,
3825 sys->Sr, GRETL_MOD_NONE);
3826 }
3827 gretl_matrix_free(G);
3828 } else {
3829 gretl_matrix_copy_values(sys->Sr, (S != NULL)? S : sys->S);
3830 }
3831
3832 #if SYSDEBUG
3833 gretl_matrix_print(sys->Sr, "G^{-1} S G^{-1}'");
3834 #endif
3835
3836 if (S != NULL) {
3837 gretl_matrix_free(S);
3838 }
3839
3840 return err;
3841 }
3842
sys_add_structural_form(equation_system * sys,const DATASET * dset)3843 static int sys_add_structural_form (equation_system *sys,
3844 const DATASET *dset)
3845 {
3846 const int *ylist = sys->ylist;
3847 const int *xlist = sys->xlist;
3848 int ne = sys->neqns;
3849 int ni = sys->nidents;
3850 int n = ne + ni;
3851 double x = 0.0;
3852 int type, col = 0;
3853 int i, j, vj, lag;
3854 int err = 0;
3855
3856 #if SYSDEBUG
3857 printlist(ylist, "endogenous vars");
3858 printlist(xlist, "exogenous vars");
3859 #endif
3860
3861 if (ylist[0] > n) {
3862 gretl_errmsg_set("system: can't add structural form");
3863 return E_DATA;
3864 }
3865
3866 sys->order = sys_max_predet_lag(sys);
3867
3868 /* allocate coefficient matrices */
3869
3870 sys->Gamma = gretl_zero_matrix_new(n, n);
3871 if (sys->Gamma == NULL) {
3872 return E_ALLOC;
3873 }
3874
3875 if (sys->order > 0) {
3876 sys->A = gretl_zero_matrix_new(n, sys->order * n);
3877 if (sys->A == NULL) {
3878 return E_ALLOC;
3879 }
3880 }
3881
3882 if (xlist[0] > 0) {
3883 sys->B = gretl_zero_matrix_new(n, xlist[0]);
3884 if (sys->B == NULL) {
3885 return E_ALLOC;
3886 }
3887 }
3888
3889 /* process stochastic equations */
3890
3891 for (i=0; i<sys->neqns && !err; i++) {
3892 const MODEL *pmod = sys->models[i];
3893 const int *mlist = pmod->list;
3894
3895 for (j=1; j<=mlist[0] && mlist[j]!=LISTSEP; j++) {
3896 vj = mlist[j];
3897 type = categorize_variable(vj, sys, &col, &lag);
3898 x = (j > 1)? pmod->coeff[j-2] : 1.0;
3899 if (type == ENDOG) {
3900 if (j == 1) {
3901 /* left-hand side variable */
3902 gretl_matrix_set(sys->Gamma, i, col, 1.0);
3903 } else {
3904 gretl_matrix_set(sys->Gamma, i, col, -x);
3905 }
3906 } else if (type == EXOG) {
3907 gretl_matrix_set(sys->B, i, col, x);
3908 } else if (type == PREDET) {
3909 col += n * (lag - 1);
3910 gretl_matrix_set(sys->A, i, col, x);
3911 } else {
3912 gretl_errmsg_sprintf("structural form:\n couldn't categorize series "
3913 "%s in eqn %d\n", dset->varname[vj], i+1);
3914 err = E_DATA;
3915 }
3916 }
3917 }
3918
3919 /* process identities */
3920
3921 for (i=0; i<sys->nidents && !err; i++) {
3922 const identity *ident = sys->idents[i];
3923
3924 for (j=0; j<=ident->n_atoms; j++) {
3925 if (j == 0) {
3926 vj = ident->depvar;
3927 x = 1.0;
3928 } else {
3929 vj = ident->atoms[j-1].varnum;
3930 x = (ident->atoms[j-1].op)? -1.0 : 1.0;
3931 }
3932 type = categorize_variable(vj, sys, &col, &lag);
3933 if (type == ENDOG) {
3934 if (j == 0) {
3935 gretl_matrix_set(sys->Gamma, ne+i, col, 1.0);
3936 } else {
3937 gretl_matrix_set(sys->Gamma, ne+i, col, -x);
3938 }
3939 } else if (type == EXOG) {
3940 gretl_matrix_set(sys->B, ne+i, col, x);
3941 } else if (type == PREDET) {
3942 col += n * (lag - 1);
3943 gretl_matrix_set(sys->A, ne+i, col, x);
3944 } else {
3945 gretl_errmsg_sprintf("structural form:\n couldn't categorize series "
3946 "%s in identity %d\n", dset->varname[vj], i+1);
3947 err = E_DATA;
3948 }
3949 }
3950 }
3951
3952 if (!err) {
3953 err = sys_add_RF_covariance_matrix(sys, n);
3954 if (err) {
3955 fprintf(stderr, "error %d in sys_add_RF_covariance_matrix\n", err);
3956 }
3957 }
3958
3959 #if SYSDEBUG
3960 gretl_matrix_print(sys->Gamma, "sys->Gamma");
3961
3962 if (sys->A != NULL) {
3963 gretl_matrix_print(sys->A, "sys->A");
3964 } else {
3965 fputs("sys->A: no lagged endog variables used as instruments\n\n", stderr);
3966 }
3967
3968 if (sys->B != NULL) {
3969 gretl_matrix_print(sys->B, "sys->B");
3970 } else {
3971 fputs("sys->B: no truly exogenous variables present\n", stderr);
3972 }
3973 #endif
3974
3975 return err;
3976 }
3977
sys_companion_matrix(equation_system * sys,int * err)3978 static gretl_matrix *sys_companion_matrix (equation_system *sys,
3979 int *err)
3980 {
3981 int m = sys->A->rows;
3982 int n = sys->A->cols;
3983 gretl_matrix *C;
3984
3985 if (m == n) {
3986 C = sys->A;
3987 } else {
3988 C = gretl_zero_matrix_new(n, n);
3989 if (C == NULL) {
3990 *err = E_ALLOC;
3991 } else {
3992 gretl_matrix_inscribe_matrix(C, sys->A, 0, 0,
3993 GRETL_MOD_NONE);
3994 gretl_matrix_inscribe_I(C, m, 0, n - m);
3995 }
3996 }
3997
3998 #if SYSDEBUG
3999 gretl_matrix_print(C, "system companion matrix");
4000 #endif
4001
4002 return C;
4003 }
4004
sys_get_fcast_se(equation_system * sys,int periods,int * err)4005 static gretl_matrix *sys_get_fcast_se (equation_system *sys,
4006 int periods, int *err)
4007 {
4008 int n = sys->neqns + sys->nidents;
4009 gretl_matrix *Tmp = NULL;
4010 gretl_matrix *V0 = NULL, *Vt = NULL;
4011 gretl_matrix *C = NULL, *se = NULL;
4012 double vti;
4013 int i, k, t;
4014
4015 if (periods <= 0) {
4016 fprintf(stderr, "Invalid number of periods\n");
4017 *err = E_DATA;
4018 return NULL;
4019 }
4020
4021 if (sys->A == NULL) {
4022 fprintf(stderr, "sys->A is NULL\n");
4023 *err = E_DATA;
4024 return NULL;
4025 }
4026
4027 C = sys_companion_matrix(sys, err);
4028 if (*err) {
4029 return NULL;
4030 }
4031
4032 se = gretl_zero_matrix_new(periods, n);
4033 if (se == NULL) {
4034 *err = E_ALLOC;
4035 if (C != sys->A) {
4036 gretl_matrix_free(C);
4037 }
4038 return NULL;
4039 }
4040
4041 k = sys->A->cols;
4042
4043 Vt = gretl_matrix_alloc(k, k);
4044 V0 = gretl_zero_matrix_new(k, k);
4045 Tmp = gretl_matrix_alloc(k, k);
4046
4047 if (Vt == NULL || V0 == NULL || Tmp == NULL) {
4048 *err = E_ALLOC;
4049 goto bailout;
4050 }
4051
4052 for (t=0; t<periods; t++) {
4053 if (t == 0) {
4054 /* initial variance */
4055 gretl_matrix_inscribe_matrix(V0, sys->Sr, 0, 0, GRETL_MOD_NONE);
4056 gretl_matrix_copy_values(Vt, V0);
4057 } else {
4058 /* calculate further variances */
4059 gretl_matrix_copy_values(Tmp, Vt);
4060 gretl_matrix_qform(C, GRETL_MOD_NONE,
4061 Tmp, Vt, GRETL_MOD_NONE);
4062 gretl_matrix_add_to(Vt, V0);
4063 }
4064
4065 for (i=0; i<n; i++) {
4066 vti = gretl_matrix_get(Vt, i, i);
4067 gretl_matrix_set(se, t, i, sqrt(vti));
4068 }
4069 }
4070
4071 bailout:
4072
4073 gretl_matrix_free(V0);
4074 gretl_matrix_free(Vt);
4075 gretl_matrix_free(Tmp);
4076
4077 if (C != sys->A) {
4078 gretl_matrix_free(C);
4079 }
4080
4081 if (*err) {
4082 gretl_matrix_free(se);
4083 se = NULL;
4084 }
4085
4086 return se;
4087 }
4088
sys_add_fcast_variance(equation_system * sys,gretl_matrix * F,int n_static)4089 static int sys_add_fcast_variance (equation_system *sys, gretl_matrix *F,
4090 int n_static)
4091 {
4092 gretl_matrix *se = NULL;
4093 double ftj, vti;
4094 int n = sys->neqns + sys->nidents;
4095 int k = F->rows - n_static;
4096 int i, j, s;
4097 int err = 0;
4098
4099 if (k > 0) {
4100 se = sys_get_fcast_se(sys, k, &err);
4101 if (err) {
4102 goto bailout;
4103 }
4104 }
4105
4106 for (j=0; j<n; j++) {
4107 for (s=0; s<F->rows; s++) {
4108 ftj = gretl_matrix_get(F, s, j);
4109 if (na(ftj)) {
4110 gretl_matrix_set(F, s, n + j, NADBL);
4111 } else {
4112 i = s - n_static;
4113 if (i < 0) {
4114 if (j < sys->neqns) {
4115 vti = sqrt(gretl_matrix_get(sys->Sr, j, j));
4116 } else {
4117 /* LHS of identity */
4118 vti = 0.0;
4119 }
4120 } else {
4121 vti = gretl_matrix_get(se, i, j);
4122 }
4123 gretl_matrix_set(F, s, n + j, vti);
4124 }
4125 }
4126 }
4127
4128 if (se != NULL) {
4129 gretl_matrix_free(se);
4130 }
4131
4132 bailout:
4133
4134 if (err) {
4135 for (i=0; i<n; i++) {
4136 for (s=0; s<F->rows; s++) {
4137 gretl_matrix_set(F, s, n + i, NADBL);
4138 }
4139 }
4140 }
4141
4142 return err;
4143 }
4144
4145 /* Experimental: generate "fitted values " (in-sample, prior to the
4146 forecast period) for variables that do not appear on the left-hand
4147 side of any stochastic equation. But I'm not quite sure how these
4148 should be defined.
4149 */
4150
sys_get_fitted_values(equation_system * sys,int v,int t1,int t2,const DATASET * dset,int * err)4151 gretl_matrix *sys_get_fitted_values (equation_system *sys,
4152 int v, int t1, int t2,
4153 const DATASET *dset,
4154 int *err)
4155 {
4156 gretl_matrix *F = NULL;
4157 gretl_matrix *G = NULL;
4158 gretl_matrix *yl = NULL, *x = NULL;
4159 gretl_matrix *y = NULL, *yh = NULL;
4160 const int *ylist = sys->ylist;
4161 const int *xlist = sys->xlist;
4162 const int *plist = sys->plist;
4163 int n = sys->neqns + sys->nidents;
4164 double xit;
4165 int col, T;
4166 int i, vi, s, t, lag;
4167
4168 if (sys->Gamma == NULL) {
4169 *err = E_DATA;
4170 return NULL;
4171 }
4172
4173 #if SYSDEBUG
4174 printlist(ylist, "ylist");
4175 printlist(xlist, "xlist");
4176 printlist(plist, "plist");
4177 fprintf(stderr, "sys->order = %d\n", sys->order);
4178 #endif
4179
4180 if (!gretl_is_identity_matrix(sys->Gamma)) {
4181 G = gretl_matrix_copy(sys->Gamma);
4182 if (G == NULL) {
4183 *err = E_ALLOC;
4184 } else {
4185 *err = gretl_SVD_invert_matrix(G);
4186 }
4187 if (*err) {
4188 gretl_matrix_free(G);
4189 return NULL;
4190 }
4191 }
4192
4193 T = t2 - t1 + 1;
4194
4195 y = gretl_matrix_alloc(n, 1);
4196 yh = gretl_matrix_alloc(n, 1);
4197 F = gretl_zero_matrix_new(T, 1);
4198
4199 if (y == NULL || yh == NULL || F == NULL) {
4200 *err = E_ALLOC;
4201 goto bailout;
4202 }
4203
4204 if (sys->order > 0) {
4205 yl = gretl_matrix_alloc(sys->order * ylist[0], 1);
4206 if (yl == NULL) {
4207 *err = E_ALLOC;
4208 goto bailout;
4209 }
4210 }
4211
4212 if (xlist[0] > 0) {
4213 x = gretl_matrix_alloc(xlist[0], 1);
4214 if (x == NULL) {
4215 *err = E_ALLOC;
4216 goto bailout;
4217 }
4218 }
4219
4220 for (t=t1, s=0; t<=t2; t++, s++) {
4221 int miss = 0;
4222
4223 /* lags of endogenous vars */
4224 if (sys->order > 0) {
4225 gretl_matrix_zero(yl);
4226 for (i=1; i<=plist[0] && !miss; i++) {
4227 vi = plist[i];
4228 get_col_and_lag(vi, sys, &col, &lag);
4229 xit = dset->Z[vi][t];
4230 col += n * (lag - 1);
4231 if (na(xit)) {
4232 miss = 1;
4233 } else {
4234 yl->val[col] = xit;
4235 }
4236 }
4237 if (!miss) {
4238 gretl_matrix_multiply(sys->A, yl, y);
4239 }
4240 } else {
4241 gretl_matrix_zero(y);
4242 }
4243
4244 #if SYSDEBUG > 1
4245 gretl_matrix_print(yl, "yl");
4246 #endif
4247
4248 /* exogenous vars */
4249 if (xlist[0] > 0 && !miss) {
4250 for (i=1; i<=xlist[0] && !miss; i++) {
4251 xit = dset->Z[xlist[i]][t];
4252 if (na(xit)) {
4253 miss = 1;
4254 } else {
4255 x->val[i-1] = xit;
4256 }
4257 }
4258 if (!miss) {
4259 gretl_matrix_multiply_mod(sys->B, GRETL_MOD_NONE,
4260 x, GRETL_MOD_NONE,
4261 y, GRETL_MOD_CUMULATE);
4262 }
4263 }
4264
4265 if (miss) {
4266 gretl_vector_set(sys->F, s, NADBL);
4267 } else {
4268 /* multiply by Gamma^{-1} */
4269 if (G != NULL) {
4270 gretl_matrix_multiply(G, y, yh);
4271 } else {
4272 gretl_matrix_multiply(sys->Gamma, y, yh);
4273 }
4274 gretl_vector_set(F, s, yh->val[v]);
4275 }
4276 }
4277
4278 bailout:
4279
4280 gretl_matrix_free(y);
4281 gretl_matrix_free(yh);
4282 gretl_matrix_free(yl);
4283 gretl_matrix_free(x);
4284 gretl_matrix_free(G);
4285
4286 if (*err) {
4287 gretl_matrix_free(F);
4288 F = NULL;
4289 } else {
4290 gretl_matrix_set_t1(F, t1);
4291 gretl_matrix_set_t2(F, t2);
4292 }
4293
4294 #if SYSDEBUG
4295 gretl_matrix_print(F, "F: calculated fitted values");
4296 #endif
4297
4298 return F;
4299 }
4300
4301 /* trim the first sys->order rows off sys->F */
4302
sys_F_trim(equation_system * sys)4303 static int sys_F_trim (equation_system *sys)
4304 {
4305 int err, p = sys->order;
4306 int r = sys->F->rows - p;
4307 gretl_matrix *tmp = gretl_matrix_alloc(r, sys->F->cols);
4308
4309 if (tmp == NULL) {
4310 err = E_ALLOC;
4311 } else {
4312 err = gretl_matrix_extract_matrix(tmp, sys->F, p, 0,
4313 GRETL_MOD_NONE);
4314 }
4315
4316 if (!err) {
4317 gretl_matrix_free(sys->F);
4318 sys->F = tmp;
4319 }
4320
4321 return err;
4322 }
4323
sys_add_forecast(equation_system * sys,int t1,int t2,const DATASET * dset,gretlopt opt)4324 static int sys_add_forecast (equation_system *sys,
4325 int t1, int t2,
4326 const DATASET *dset,
4327 gretlopt opt)
4328 {
4329 gretl_matrix *G = NULL;
4330 gretl_matrix *yl = NULL, *x = NULL;
4331 gretl_matrix *y = NULL, *yh = NULL;
4332 const int *ylist = sys->ylist;
4333 const int *xlist = sys->xlist;
4334 const int *plist = sys->plist;
4335 int n = sys->neqns + sys->nidents;
4336 double xit, xitd;
4337 int tdyn, col, T, ncols;
4338 int i, vi, s, t, lag, s0 = 0;
4339 int err = 0;
4340
4341 if (sys->Gamma == NULL || sys->Gamma->rows != n) {
4342 fprintf(stderr, "sys_add_forecast: Gamma is broken!\n");
4343 return E_DATA;
4344 }
4345
4346 #if SYSDEBUG
4347 fprintf(stderr, "*** sys_add_forecast\n");
4348 printlist(ylist, "ylist");
4349 printlist(xlist, "xlist");
4350 printlist(plist, "plist");
4351 fprintf(stderr, "sys->order = %d\n", sys->order);
4352 fprintf(stderr, "t1=%d, t2=%d, out-of-sample=%d\n",
4353 t1, t2, (opt & OPT_O)? 1 : 0);
4354 #endif
4355
4356 if (!gretl_is_identity_matrix(sys->Gamma)) {
4357 G = gretl_matrix_copy(sys->Gamma);
4358 if (G == NULL) {
4359 err = E_ALLOC;
4360 } else {
4361 err = gretl_SVD_invert_matrix(G);
4362 }
4363 if (err) {
4364 gretl_matrix_free(G);
4365 return err;
4366 }
4367 }
4368
4369 /* At which observation (tdyn) does the forecast turn
4370 dynamic, if at all? */
4371
4372 if ((opt & OPT_S) || sys->A == NULL) {
4373 /* got the --static option (or no dynamics): never */
4374 tdyn = t2 + 1;
4375 } else if (opt & OPT_D) {
4376 /* got the --dynamic option: from the start */
4377 tdyn = t1;
4378 } else {
4379 /* by default, for a model with dynamics: just after
4380 the estimation sample ends */
4381 tdyn = sys->t2 + 1;
4382 }
4383
4384 T = t2 - t1 + 1;
4385 ncols = 2 * n;
4386
4387 if (sys->order > 0 && (opt & OPT_O)) {
4388 /* out-of-sample option: we'll need some extra rows
4389 for pre-sample values
4390 */
4391 T += sys->order;
4392 }
4393
4394 y = gretl_matrix_alloc(n, 1);
4395 yh = gretl_matrix_alloc(n, 1);
4396 sys->F = gretl_zero_matrix_new(T, ncols);
4397
4398 if (y == NULL || yh == NULL || sys->F == NULL) {
4399 err = E_ALLOC;
4400 goto bailout;
4401 }
4402
4403 if (sys->order > 0) {
4404 yl = gretl_matrix_alloc(sys->order * ylist[0], 1);
4405 if (yl == NULL) {
4406 err = E_ALLOC;
4407 goto bailout;
4408 }
4409 if (opt & OPT_O) {
4410 /* out-of-sample */
4411 int p = sys->order;
4412 double yti;
4413
4414 for (t=t1-p, s=0; t<t1; t++, s++) {
4415 for (i=0; i<n; i++) {
4416 yti = dset->Z[ylist[i+1]][t];
4417 gretl_matrix_set(sys->F, s, i, yti);
4418 }
4419 }
4420 s0 = p;
4421 }
4422 }
4423
4424 if (xlist[0] > 0) {
4425 x = gretl_matrix_alloc(xlist[0], 1);
4426 if (x == NULL) {
4427 err = E_ALLOC;
4428 goto bailout;
4429 }
4430 }
4431
4432 for (t=t1, s=s0; t<=t2; t++, s++) {
4433 int miss = 0;
4434
4435 /* lags of endogenous vars */
4436 if (sys->order > 0) {
4437 gretl_matrix_zero(yl);
4438 for (i=1; i<=plist[0] && !miss; i++) {
4439 vi = plist[i];
4440 get_col_and_lag(vi, sys, &col, &lag);
4441 xitd = NADBL;
4442 if (t < tdyn || s - lag < 0) {
4443 /* pre-forecast value */
4444 xit = dset->Z[vi][t];
4445 } else {
4446 /* prior forecast value preferred */
4447 if (s - lag >= 0) {
4448 xitd = xit = dset->Z[vi][t];
4449 }
4450 xit = gretl_matrix_get(sys->F, s - lag, col);
4451 }
4452 col += n * (lag - 1);
4453 if (!na(xit)) {
4454 yl->val[col] = xit;
4455 } else if (!na(xitd)) {
4456 yl->val[col] = xitd;
4457 } else {
4458 miss = 1;
4459 }
4460 }
4461 if (!miss) {
4462 gretl_matrix_multiply(sys->A, yl, y);
4463 }
4464 } else {
4465 gretl_matrix_zero(y);
4466 }
4467
4468 #if SYSDEBUG > 1
4469 gretl_matrix_print(yl, "yl");
4470 #endif
4471
4472 /* exogenous vars */
4473 if (xlist[0] > 0 && !miss) {
4474 for (i=1; i<=xlist[0] && !miss; i++) {
4475 xit = dset->Z[xlist[i]][t];
4476 if (na(xit)) {
4477 miss = 1;
4478 } else {
4479 x->val[i-1] = xit;
4480 }
4481 }
4482 if (!miss) {
4483 gretl_matrix_multiply_mod(sys->B, GRETL_MOD_NONE,
4484 x, GRETL_MOD_NONE,
4485 y, GRETL_MOD_CUMULATE);
4486 }
4487 }
4488
4489 if (miss) {
4490 for (i=0; i<n; i++) {
4491 gretl_matrix_set(sys->F, s, i, NADBL);
4492 }
4493 } else {
4494 /* multiply by Gamma^{-1} */
4495 if (G != NULL) {
4496 gretl_matrix_multiply(G, y, yh);
4497 } else {
4498 gretl_matrix_multiply(sys->Gamma, y, yh);
4499 }
4500 for (i=0; i<n; i++) {
4501 gretl_matrix_set(sys->F, s, i, yh->val[i]);
4502 }
4503 }
4504 }
4505
4506 bailout:
4507
4508 gretl_matrix_free(y);
4509 gretl_matrix_free(yh);
4510 gretl_matrix_free(yl);
4511 gretl_matrix_free(x);
4512 gretl_matrix_free(G);
4513
4514 if (!err && s0 > 0) {
4515 err = sys_F_trim(sys);
4516 }
4517
4518 if (err) {
4519 gretl_matrix_free(sys->F);
4520 sys->F = NULL;
4521 } else {
4522 gretl_matrix_set_t1(sys->F, t1);
4523 gretl_matrix_set_t2(sys->F, t2);
4524 }
4525
4526 #if SYSDEBUG
4527 gretl_matrix_print(sys->F, "sys->F, forecasts only");
4528 #endif
4529
4530 if (!err) {
4531 sys_add_fcast_variance(sys, sys->F, tdyn - t1);
4532 }
4533
4534 #if SYSDEBUG
4535 gretl_matrix_print(sys->F, "sys->F, with variance");
4536 #endif
4537
4538 return err;
4539 }
4540
4541 const gretl_matrix *
system_get_forecast_matrix(equation_system * sys,int t1,int t2,DATASET * dset,gretlopt opt,int * err)4542 system_get_forecast_matrix (equation_system *sys, int t1, int t2,
4543 DATASET *dset, gretlopt opt,
4544 int *err)
4545 {
4546 if (sys->F != NULL) {
4547 gretl_matrix_replace(&sys->F, NULL);
4548 }
4549
4550 *err = sys_add_forecast(sys, t1, t2, dset, opt);
4551
4552 return sys->F;
4553 }
4554
4555 /* attach to the system some extra statistics generated in the course
4556 of estimation via LIML */
4557
sys_attach_ldata(equation_system * sys)4558 static int sys_attach_ldata (equation_system *sys)
4559 {
4560 int i, n = sys->neqns;
4561 int err = 0;
4562
4563 sys->ldata = malloc(sizeof *sys->ldata);
4564
4565 if (sys->ldata == NULL) {
4566 err = E_ALLOC;
4567 } else {
4568 sys->ldata->lmin = NULL;
4569 sys->ldata->ll = NULL;
4570 sys->ldata->idf = NULL;
4571
4572 sys->ldata->lmin = malloc(n * sizeof *sys->ldata->lmin);
4573 sys->ldata->ll = malloc(n * sizeof *sys->ldata->ll);
4574 sys->ldata->idf = malloc(n * sizeof *sys->ldata->idf);
4575
4576 if (sys->ldata->lmin == NULL ||
4577 sys->ldata->ll == NULL ||
4578 sys->ldata->idf == NULL) {
4579 free(sys->ldata->lmin);
4580 free(sys->ldata->ll);
4581 free(sys->ldata->idf);
4582 free(sys->ldata);
4583 sys->ldata = NULL;
4584 err = E_ALLOC;
4585 }
4586 }
4587
4588 if (!err) {
4589 const MODEL *pmod;
4590
4591 for (i=0; i<n; i++) {
4592 pmod = sys->models[i];
4593 sys->ldata->lmin[i] = gretl_model_get_double(pmod, "lmin");
4594 sys->ldata->ll[i] = pmod->lnL;
4595 sys->ldata->idf[i] = gretl_model_get_int(pmod, "idf");
4596 }
4597 }
4598
4599 return err;
4600 }
4601
sys_print_reconstituted_models(const equation_system * sys,const DATASET * dset,PRN * prn)4602 static int sys_print_reconstituted_models (const equation_system *sys,
4603 const DATASET *dset,
4604 PRN *prn)
4605 {
4606 MODEL mod;
4607 const double *y;
4608 double x;
4609 int print_insts = 0;
4610 int ifc = 0, nc = 0, ncmax = 0;
4611 int i, j, t, k = 0;
4612 int err = 0;
4613
4614 gretl_model_init(&mod, dset);
4615
4616 mod.t1 = sys->t1;
4617 mod.t2 = sys->t2;
4618 mod.nobs = sys->T;
4619 mod.aux = AUX_SYS;
4620 gretl_model_set_int(&mod, "method", sys->method);
4621
4622 mod.lnL = NADBL;
4623
4624 if (sys->method == SYS_METHOD_TSLS ||
4625 sys->method == SYS_METHOD_3SLS ||
4626 sys->method == SYS_METHOD_FIML ||
4627 sys->method == SYS_METHOD_LIML) {
4628 mod.ci = IVREG;
4629 } else {
4630 mod.ci = OLS;
4631 }
4632
4633 print_insts = sys->method == SYS_METHOD_TSLS ||
4634 sys->method == SYS_METHOD_3SLS;
4635
4636 for (i=0; i<sys->neqns && !err; i++) {
4637 int *mlist = sys->lists[i];
4638 int freelist = 0;
4639
4640 mod.ID = i;
4641
4642 if (print_insts && !gretl_list_has_separator(mlist)) {
4643 mod.list = compose_ivreg_list(sys, i);
4644 if (mod.list == NULL) {
4645 err = E_ALLOC;
4646 } else {
4647 freelist = 1;
4648 }
4649 } else {
4650 mod.list = mlist;
4651 }
4652
4653 if (!err) {
4654 ifc = nc = 0;
4655 for (j=2; j<=mod.list[0]; j++) {
4656 if (mod.list[j] == 0) {
4657 ifc = 1;
4658 } else if (mod.list[j] == LISTSEP) {
4659 break;
4660 }
4661 nc++;
4662 }
4663 }
4664
4665 if (!err && nc > ncmax) {
4666 mod.coeff = realloc(mod.coeff, nc * sizeof *mod.coeff);
4667 mod.sderr = realloc(mod.sderr, nc * sizeof *mod.sderr);
4668 ncmax = nc;
4669 }
4670
4671 if (mod.coeff == NULL || mod.sderr == NULL) {
4672 err = E_ALLOC;
4673 break;
4674 }
4675
4676 mod.ncoeff = nc;
4677 mod.dfn = nc - ifc;
4678
4679 if (sys->flags & SYSTEM_DFCORR) {
4680 gretl_model_set_int(&mod, "dfcorr", 1);
4681 mod.dfd = mod.nobs - nc;
4682 } else {
4683 mod.dfd = mod.nobs;
4684 }
4685
4686 y = dset->Z[mod.list[1]];
4687 mod.ybar = gretl_mean(mod.t1, mod.t2, y);
4688 mod.sdy = gretl_stddev(mod.t1, mod.t2, y);
4689
4690 mod.ess = 0.0;
4691 for (t=0; t<sys->T; t++) {
4692 x = gretl_matrix_get(sys->E, t, i);
4693 mod.ess += x * x;
4694 }
4695
4696 if (sys->method == SYS_METHOD_OLS ||
4697 sys->method == SYS_METHOD_TSLS ||
4698 sys->method == SYS_METHOD_LIML) {
4699 /* single-equation methods */
4700 mod.sigma = sqrt(mod.ess / mod.dfd);
4701 } else {
4702 mod.sigma = sqrt(gretl_matrix_get(sys->S, i, i));
4703 }
4704
4705 for (j=0; j<mod.ncoeff; j++) {
4706 mod.coeff[j] = sys->b->val[k];
4707 mod.sderr[j] = sqrt(gretl_matrix_get(sys->vcv, k, k));
4708 k++;
4709 }
4710
4711 if (sys->method == SYS_METHOD_LIML && sys->ldata != NULL) {
4712 gretl_model_set_double(&mod, "lmin", sys->ldata->lmin[i]);
4713 mod.lnL = sys->ldata->ll[i];
4714 gretl_model_set_int(&mod, "idf", sys->ldata->idf[i]);
4715 }
4716
4717 printmodel(&mod, dset, OPT_NONE, prn);
4718
4719 if (freelist) {
4720 free(mod.list);
4721 }
4722
4723 mod.list = NULL;
4724 }
4725
4726 clear_model(&mod);
4727
4728 return err;
4729 }
4730
gretl_system_print(equation_system * sys,const DATASET * dset,gretlopt opt,PRN * prn)4731 int gretl_system_print (equation_system *sys, const DATASET *dset,
4732 gretlopt opt, PRN *prn)
4733 {
4734 const char *name = sys->name;
4735 gchar *sysstr = NULL;
4736 int tex = tex_format(prn);
4737 int nr = system_n_restrictions(sys);
4738 int i;
4739
4740 if (sys->models != NULL &&
4741 sys->method == SYS_METHOD_LIML &&
4742 sys->ldata == NULL) {
4743 sys_attach_ldata(sys);
4744 }
4745
4746 if (name != NULL && sys_anonymous(name)) {
4747 /* don't print internal reference name */
4748 name = NULL;
4749 }
4750
4751 if (tex) {
4752 pputs(prn, "\\begin{center}\n");
4753 } else {
4754 pputc(prn, '\n');
4755 }
4756
4757 sysstr = system_get_full_string(sys);
4758 if (name != NULL) {
4759 if (tex) {
4760 pprintf(prn, "%s, %s\\\\\n", _("Equation system"), name);
4761 } else {
4762 pprintf(prn, "%s, %s\n", _("Equation system"), name);
4763 }
4764 pprintf(prn, "%s: %s", _("Estimator"), sysstr);
4765 } else {
4766 pprintf(prn, "%s, %s", _("Equation system"), sysstr);
4767 }
4768 g_free(sysstr);
4769
4770 if (sys->iters > 0) {
4771 gretl_prn_newline(prn);
4772 pprintf(prn, _("Convergence achieved after %d iterations\n"), sys->iters);
4773 if (sys->method == SYS_METHOD_SUR ||
4774 sys->method == SYS_METHOD_FIML) {
4775 if (tex) {
4776 gretl_prn_newline(prn);
4777 pprintf(prn, "%s = ", _("Log-likelihood"));
4778 if (sys->ll < 0) {
4779 pprintf(prn, "$-$%g", -sys->ll);
4780 } else {
4781 pprintf(prn, "%g", sys->ll);
4782 }
4783 } else {
4784 pprintf(prn, "%s = %g\n", _("Log-likelihood"), sys->ll);
4785 }
4786 }
4787 }
4788
4789 if (tex) {
4790 pputs(prn, "\n\\end{center}\n\n");
4791 } else {
4792 pputc(prn, '\n');
4793 }
4794
4795 if (sys->models != NULL) {
4796 for (i=0; i<sys->neqns; i++) {
4797 if (sys->flags & SYSTEM_DFCORR) {
4798 gretl_model_set_int(sys->models[i], "dfcorr", 1);
4799 }
4800 printmodel(sys->models[i], dset, OPT_NONE, prn);
4801 }
4802 } else {
4803 sys_print_reconstituted_models(sys, dset, prn);
4804 }
4805
4806 system_print_sigma(sys, prn);
4807
4808 if (nr == 0 && (sys->method == SYS_METHOD_FIML ||
4809 sys->method == SYS_METHOD_3SLS ||
4810 sys->method == SYS_METHOD_SUR)) {
4811 print_system_overid_test(sys, prn);
4812 }
4813
4814 return 0;
4815 }
4816
ensure_asy_printout(equation_system * sys)4817 static void ensure_asy_printout (equation_system *sys)
4818 {
4819 int i;
4820
4821 if (sys->models != NULL) {
4822 for (i=0; i<sys->neqns; i++) {
4823 sys->models[i]->ci = IVREG;
4824 }
4825 }
4826 }
4827
4828 int
system_save_and_print_results(equation_system * sys,DATASET * dset,gretlopt opt,PRN * prn)4829 system_save_and_print_results (equation_system *sys, DATASET *dset,
4830 gretlopt opt, PRN *prn)
4831 {
4832 int nr = system_n_restrictions(sys);
4833 int err = 0;
4834
4835 if (sys->E != NULL) {
4836 gretl_matrix_set_t1(sys->E, sys->t1);
4837 gretl_matrix_set_t2(sys->E, sys->t2);
4838 }
4839
4840 if (sys->iters > 0 && sys->method == SYS_METHOD_SUR && nr == 0) {
4841 sur_ols_diag(sys);
4842 }
4843
4844 sys->ldet = gretl_vcv_log_determinant(sys->S, &err);
4845
4846 if (!err) {
4847 err = system_add_yhat_matrix(sys);
4848 }
4849
4850 if (!err) {
4851 err = sys_add_structural_form(sys, dset);
4852 }
4853
4854 if (!(opt & OPT_Q)) {
4855 if (sys->method == SYS_METHOD_FIML) {
4856 ensure_asy_printout(sys);
4857 }
4858 gretl_system_print(sys, dset, opt, prn);
4859 }
4860
4861 return err;
4862 }
4863
system_autocorrelation_test(equation_system * sys,int order,gretlopt opt,PRN * prn)4864 int system_autocorrelation_test (equation_system *sys, int order,
4865 gretlopt opt, PRN *prn)
4866 {
4867 double *u, lb;
4868 int i, err = 0;
4869
4870 for (i=0; i<sys->neqns && !err; i++) {
4871 if (!(opt & OPT_Q)) {
4872 pprintf(prn, "%s %d:\n", _("Equation"), i + 1);
4873 }
4874 u = sys->E->val + (i * sys->T);
4875 lb = ljung_box(order, 0, sys->T - 1, u, &err);
4876 if (!err && !((opt & OPT_Q))) {
4877 pprintf(prn, "%s: %s(%d) = %g [%.4f]\n\n",
4878 _("Ljung-Box Q'"), _("Chi-square"), order,
4879 lb, chisq_cdf_comp(order, lb));
4880 }
4881 }
4882
4883 return err;
4884 }
4885
system_arch_test(equation_system * sys,int order,gretlopt opt,PRN * prn)4886 int system_arch_test (equation_system *sys, int order,
4887 gretlopt opt, PRN *prn)
4888 {
4889 const double *u;
4890 int i, err = 0;
4891
4892 if (!(opt & OPT_I)) {
4893 pputc(prn, '\n');
4894 }
4895
4896 for (i=0; i<sys->neqns && !err; i++) {
4897 if (!(opt & OPT_I)) {
4898 pprintf(prn, "%s %d:\n", _("Equation"), i + 1);
4899 }
4900 u = sys->E->val + (i * sys->T);
4901 err = array_arch_test(u, sys->T, order, opt | OPT_M, prn);
4902 }
4903
4904 return err;
4905 }
4906
system_supports_method(equation_system * sys,int method)4907 int system_supports_method (equation_system *sys, int method)
4908 {
4909 if (sys != NULL) {
4910 int i;
4911
4912 for (i=0; i<sys->neqns; i++) {
4913 if (gretl_list_has_separator(sys->lists[i])) {
4914 return method == SYS_METHOD_TSLS ||
4915 method == SYS_METHOD_3SLS;
4916 }
4917 }
4918 }
4919
4920 return 1;
4921 }
4922
finalize_liml_model(MODEL * pmod,equation_system * sys)4923 static void finalize_liml_model (MODEL *pmod, equation_system *sys)
4924 {
4925 *pmod = *sys->models[0];
4926
4927 #if SYSDEBUG
4928 display_model_data_items(pmod);
4929 #endif
4930
4931 gretl_model_destroy_data_item(pmod, "tslsX");
4932 gretl_model_destroy_data_item(pmod, "endog");
4933 gretl_model_destroy_data_item(pmod, "method");
4934 gretl_model_destroy_data_item(pmod, "liml_y");
4935
4936 #if SYSDEBUG
4937 display_model_data_items(pmod);
4938 #endif
4939
4940 free(sys->models[0]);
4941 free(sys->models);
4942 sys->models = NULL;
4943
4944 pmod->aux = AUX_NONE;
4945 pmod->opt |= OPT_L;
4946 pmod->rsq = pmod->adjrsq = NADBL;
4947 pmod->fstt = NADBL;
4948 set_model_id(pmod, OPT_NONE);
4949 }
4950
4951 /* implement the --liml option to "tsls" */
4952
single_equation_liml(const int * list,DATASET * dset,gretlopt opt)4953 MODEL single_equation_liml (const int *list, DATASET *dset,
4954 gretlopt opt)
4955 {
4956 int *mlist = NULL, *ilist = NULL;
4957 equation_system *sys = NULL;
4958 MODEL model;
4959 int err = 0;
4960
4961 gretl_model_init(&model, dset);
4962
4963 err = ivreg_process_lists(list, &mlist, &ilist);
4964
4965 if (!err) {
4966 sys = equation_system_new(SYS_METHOD_LIML, NULL, &err);
4967 }
4968
4969 if (!err) {
4970 err = equation_system_append(sys, mlist);
4971 }
4972
4973 if (!err) {
4974 sys->ilist = ilist;
4975 err = equation_system_finalize(sys, dset, OPT_S, NULL);
4976 }
4977
4978 if (err) {
4979 model.errcode = err;
4980 } else {
4981 finalize_liml_model(&model, sys);
4982 }
4983
4984 equation_system_destroy(sys);
4985 free(mlist);
4986
4987 return model;
4988 }
4989
gretl_system_get_sample(const equation_system * sys,int * t1,int * t2)4990 int gretl_system_get_sample (const equation_system *sys,
4991 int *t1, int *t2)
4992 {
4993 if (sys != NULL) {
4994 *t1 = sys->smpl_t1;
4995 *t2 = sys->smpl_t2;
4996 return 0;
4997 }
4998
4999 return E_DATA;
5000 }
5001