1 /*
2 * gretl -- Gnu Regression, Econometrics and Time-series Library
3 * Copyright (C) 2001 Allin Cottrell and Riccardo "Jack" Lucchetti
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, either version 3 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program. If not, see <http://www.gnu.org/licenses/>.
17 *
18 */
19
20 #include "libgretl.h"
21 #include "system.h"
22 #include "var.h"
23 #include "objstack.h"
24 #include "uservar.h"
25 #include "matrix_extra.h"
26 #include "gretl_restrict.h"
27 #include "bootstrap.h"
28 #include "gretl_func.h"
29 #include "gretl_bfgs.h"
30 #include "cmd_private.h"
31
32 #define RDEBUG 0
33 #define R_UNSPEC -9
34
35 enum {
36 VECM_NONE = 0,
37 VECM_B = 1 << 0,
38 VECM_A = 1 << 1
39 };
40
41 typedef struct rrow_ rrow;
42
43 struct rrow_ {
44 int nterms; /* number of terms in restriction */
45 double *mult; /* array of numerical multipliers on coeffs */
46 int *eq; /* array of equation numbers (for multi-equation case) */
47 int *bnum; /* array of coeff numbers */
48 char *letter; /* array of coeff letters */
49 double rhs; /* numerical value on right-hand side */
50 };
51
52 struct gretl_restriction_ {
53 int g; /* number of restrictions (rows of R) */
54 int gmax; /* max. possible restrictions */
55 int bmulti; /* flag for case where coeffs need two indices */
56 int amulti; /* VECM only */
57 int gb, ga; /* restrictions on beta, alpha (VECM only) */
58 int bcols, acols; /* VECM only */
59 int vecm; /* pertains to VECM beta or alpha? */
60 gretl_matrix *R; /* LHS restriction matrix */
61 gretl_matrix *q; /* RHS restriction matrix */
62 gretl_matrix *Ra; /* second LHS restriction matrix */
63 gretl_matrix *qa; /* second RHS restriction matrix */
64 char *mask; /* selection mask for coeffs */
65 rrow **rows; /* "atomic" restrictions */
66 void *obj; /* pointer to model, system */
67 GretlObjType otype; /* type of model, system */
68 gretlopt opt; /* OPT_C is used for "coeffsum" */
69 char *rfunc; /* name of nonlinear restriction function */
70 double test; /* test statistic */
71 double pval; /* p-value of test statistic */
72 double lnl; /* loglikelihood */
73 double bsum; /* sum of coefficients */
74 double bse; /* standard error of bsum */
75 int code; /* identifier for statistic type */
76 };
77
78 #define targ_specified(r,i,j,c) (((c=='b' && r->bmulti) || (c=='a' && r->amulti)) \
79 && r->rows[i]->eq != NULL \
80 && r->rows[i]->eq[j] != R_UNSPEC)
81
82 static int
83 real_restriction_set_parse_line (gretl_restriction *rset,
84 const char *line,
85 const DATASET *dset,
86 int first);
87
88 /* @R is putatively the matrix in R \beta -q = 0. Check that it makes
89 sense in that role, i.e., that RR' is non-singular. */
90
check_R_matrix(const gretl_matrix * R)91 static int check_R_matrix (const gretl_matrix *R)
92 {
93 gretl_matrix *m;
94 int err = 0;
95
96 m = gretl_matrix_alloc(R->rows, R->rows);
97
98 if (m == NULL) {
99 err = E_ALLOC;
100 } else {
101 gretl_matrix_multiply_mod(R, GRETL_MOD_NONE,
102 R, GRETL_MOD_TRANSPOSE,
103 m, GRETL_MOD_NONE);
104
105 err = gretl_invert_general_matrix(m);
106
107 if (err == E_SINGULAR) {
108 gretl_errmsg_set(_("Matrix inversion failed: restrictions may be "
109 "inconsistent or redundant"));
110 }
111
112 gretl_matrix_free(m);
113 }
114
115 return err;
116 }
117
118 /* Here we're adding to an existing restriction on a VECM, either in
119 respect of beta (@letter = 'b') or in respect of alpha
120 (@letter = 'a').
121 */
122
augment_vecm_restriction(gretl_restriction * rset,char letter)123 static int augment_vecm_restriction (gretl_restriction *rset,
124 char letter)
125 {
126 GRETL_VAR *vecm = rset->obj;
127 const gretl_matrix *R0, *q0;
128 gretl_matrix *R2 = NULL, *q2 = NULL;
129 int err = 0;
130
131 if (letter == 'b') {
132 R0 = gretl_VECM_R_matrix(vecm);
133 q0 = gretl_VECM_q_matrix(vecm);
134 } else {
135 R0 = gretl_VECM_Ra_matrix(vecm);
136 q0 = gretl_VECM_qa_matrix(vecm);
137 }
138
139 R2 = gretl_matrix_row_concat(R0, rset->R, &err);
140 if (err) {
141 return err;
142 }
143
144 err = check_R_matrix(R2);
145 if (err) {
146 gretl_matrix_free(R2);
147 return err;
148 }
149
150 if (q0 == NULL) {
151 q2 = gretl_column_vector_alloc(R2->rows);
152 if (q2 == NULL) {
153 err = E_ALLOC;
154 } else {
155 int i, n = R0->rows;
156
157 for (i=0; i<R2->rows; i++) {
158 q2->val[i] = (i < n)? 0.0 : rset->q->val[i-n];
159 }
160 }
161 } else {
162 q2 = gretl_matrix_row_concat(q0, rset->q, &err);
163 }
164
165 if (err) {
166 gretl_matrix_free(R2);
167 } else if (letter == 'b') {
168 gretl_matrix_replace(&rset->R, R2);
169 gretl_matrix_replace(&rset->q, q2);
170 } else {
171 gretl_matrix_replace(&rset->Ra, R2);
172 gretl_matrix_replace(&rset->qa, q2);
173 }
174
175 return err;
176 }
177
178 /* Adding an existing restriction on a VECM, in respect of beta or
179 alpha, to a new restriction in respect of the other matrix.
180 */
181
add_old_vecm_restriction(gretl_restriction * rset,char letter)182 static int add_old_vecm_restriction (gretl_restriction *rset,
183 char letter)
184 {
185 GRETL_VAR *vecm = rset->obj;
186 int err = 0;
187
188 if (letter == 'b') {
189 const gretl_matrix *R = gretl_VECM_R_matrix(vecm);
190 const gretl_matrix *q = gretl_VECM_q_matrix(vecm);
191
192 rset->R = gretl_matrix_copy(R);
193 if (rset->R == NULL) {
194 err = E_ALLOC;
195 } else if (q != NULL) {
196 rset->q = gretl_matrix_copy(q);
197 if (rset->q == NULL) {
198 err = E_ALLOC;
199 }
200 }
201 } else {
202 const gretl_matrix *Ra = gretl_VECM_Ra_matrix(vecm);
203 const gretl_matrix *qa = gretl_VECM_qa_matrix(vecm);
204
205 rset->Ra = gretl_matrix_copy(Ra);
206 if (rset->Ra == NULL) {
207 err = E_ALLOC;
208 } else if (qa != NULL) {
209 rset->qa = gretl_matrix_copy(qa);
210 if (rset->qa == NULL) {
211 err = E_ALLOC;
212 }
213 }
214 }
215
216 return err;
217 }
218
219 /* We set things up here such that the restrictions are
220 R_b * vec(beta) = q; R_a * vec(alpha') = 0.
221 */
222
get_R_vecm_column(const gretl_restriction * rset,int i,int j,char letter)223 static int get_R_vecm_column (const gretl_restriction *rset,
224 int i, int j, char letter)
225 {
226 const rrow *r = rset->rows[i];
227 GRETL_VAR *var = rset->obj;
228 int col = 0;
229
230 if (letter == 'b') {
231 col = r->bnum[j];
232 if (rset->bmulti) {
233 col += r->eq[j] * gretl_VECM_n_beta(var);
234 }
235 } else if (letter == 'a') {
236 col = r->bnum[j];
237 if (rset->amulti) {
238 col += r->eq[j] * gretl_VECM_n_alpha(var);
239 }
240 }
241
242 return col;
243 }
244
get_R_var_column(const gretl_restriction * rset,int i,int j)245 static int get_R_var_column (const gretl_restriction *rset,
246 int i, int j)
247 {
248 const rrow *r = rset->rows[i];
249 GRETL_VAR *var = rset->obj;
250 int k, col = r->bnum[j];
251
252 for (k=0; k<r->eq[j]; k++) {
253 col += var->ncoeff;
254 }
255
256 return col;
257 }
258
get_R_sys_column(const gretl_restriction * rset,int i,int j)259 static int get_R_sys_column (const gretl_restriction *rset,
260 int i, int j)
261 {
262 const rrow *r = rset->rows[i];
263 equation_system *sys = rset->obj;
264 int col = r->bnum[j];
265 int k, n;
266
267 for (k=0; k<r->eq[j]; k++) {
268 n = system_get_list_length(sys, k);
269 if (n > 0) {
270 col += n - 1;
271 }
272 }
273
274 return col;
275 }
276
get_restriction_param(const rrow * r,int k)277 static double get_restriction_param (const rrow *r, int k)
278 {
279 int i;
280
281 for (i=0; i<r->nterms; i++) {
282 if (r->bnum[i] == k) {
283 return r->mult[i];
284 }
285 }
286
287 return 0.0;
288 }
289
vecm_cross_error(void)290 static void vecm_cross_error (void)
291 {
292 gretl_errmsg_set("VECM: beta/alpha cross restrictions are "
293 "not handled");
294 }
295
296 /* See if we have both beta and alpha terms; check that there
297 are no cross beta/alpha restrictions. */
298
vecm_ab_check(gretl_restriction * rset)299 static int vecm_ab_check (gretl_restriction *rset)
300 {
301 int atotal = 0, btotal = 0;
302 int i, j, err = 0;
303
304 for (i=0; i<rset->g && !err; i++) {
305 rrow *r = rset->rows[i];
306 int a = 0, b = 0;
307
308 for (j=0; j<r->nterms && !err; j++) {
309 if (r->letter[j] == 'a') {
310 a++;
311 atotal++;
312 } else {
313 b++;
314 btotal++;
315 }
316 if (a > 0 && b > 0) {
317 vecm_cross_error();
318 err = E_NOTIMP;
319 }
320 }
321 }
322
323 if (!err) {
324 /* incoming default is VECM_B */
325 if (atotal > 0) {
326 rset->vecm |= VECM_A;
327 }
328 if (btotal == 0) {
329 rset->vecm &= ~VECM_B;
330 }
331 }
332
333 return err;
334 }
335
336 /* check integrity of vecm restrictions, either beta or alpha */
337
vecm_x_check(gretl_restriction * rset,char letter)338 static int vecm_x_check (gretl_restriction *rset, char letter)
339 {
340 int *multi = (letter == 'a')? &rset->amulti : &rset->bmulti;
341 int *ki = (letter == 'a')? &rset->ga : &rset->gb;
342 int unspec = 0, anyspec = 0;
343 int i, j, err = 0;
344
345 for (i=0; i<rset->g && !err; i++) {
346 rrow *r = rset->rows[i];
347
348 if (r->letter[0] != letter) {
349 continue;
350 }
351
352 *ki += 1;
353
354 for (j=0; j<r->nterms && !err; j++) {
355 if (r->eq != NULL) {
356 if (r->eq[j] == R_UNSPEC) {
357 unspec = 1;
358 } else if (r->eq[j] >= 0) {
359 anyspec = 1;
360 }
361 }
362 if (anyspec && unspec) {
363 err = E_PARSE;
364 }
365 }
366 }
367
368 if (!err && unspec && *multi) {
369 for (i=0; i<rset->g; i++) {
370 if (rset->rows[i]->letter[0] == letter) {
371 free(rset->rows[i]->eq);
372 rset->rows[i]->eq = NULL;
373 }
374 }
375 *multi = 0;
376 }
377
378 return err;
379 }
380
381 /* Check the validity of a set of VECM restrictions */
382
vecm_restriction_check(gretl_restriction * rset)383 static int vecm_restriction_check (gretl_restriction *rset)
384 {
385 int r = gretl_VECM_rank(rset->obj);
386 int err;
387
388 err = vecm_ab_check(rset);
389
390 if (!err && (rset->vecm & VECM_B)) {
391 err = vecm_x_check(rset, 'b');
392 }
393
394 if (!err && (rset->vecm & VECM_A)) {
395 err = vecm_x_check(rset, 'a');
396 }
397
398 if (!err && (rset->vecm & VECM_B)) {
399 rset->bcols = gretl_VECM_n_beta(rset->obj);
400 if (rset->bmulti) {
401 rset->bcols *= r;
402 }
403 }
404
405 if (!err && (rset->vecm & VECM_A)) {
406 rset->acols = gretl_VECM_n_alpha(rset->obj);
407 if (rset->amulti) {
408 rset->acols *= r;
409 }
410 }
411
412 if (!err) {
413 if (rset->gb > rset->bcols) {
414 gretl_errmsg_set("Too many restrictions");
415 err = E_NONCONF;
416 } else if (rset->ga > rset->acols) {
417 fprintf(stderr, "rset->ga = %d, rset->acols = %d\n",
418 rset->ga, rset->acols);
419 gretl_errmsg_set("Too many restrictions");
420 err = E_NONCONF;
421 }
422 }
423
424 return err;
425 }
426
427 /* If we have been given restrictions in row-wise fashion,
428 allocate the required R and q matrices */
429
rset_allocate_R_q(gretl_restriction * rset,int nc)430 static int rset_allocate_R_q (gretl_restriction *rset, int nc)
431 {
432 rset->R = gretl_zero_matrix_new(rset->g, nc);
433 rset->q = gretl_zero_matrix_new(rset->g, 1);
434
435 #if RDEBUG
436 fprintf(stderr, "rset_allocate_R_q: on output R=%p (%dx%d), q=%p (%dx%d)\n",
437 (void *) rset->R, rset->g, nc, (void *) rset->q, rset->g, 1);
438 #endif
439
440 if (rset->R == NULL || rset->q == NULL) {
441 gretl_matrix_free(rset->R);
442 gretl_matrix_free(rset->q);
443 rset->R = rset->q = NULL;
444 return E_ALLOC;
445 }
446
447 return 0;
448 }
449
450 /* For the single-equation case, allocate and fill out the
451 R and q matrices. */
452
equation_form_matrices(gretl_restriction * rset)453 static int equation_form_matrices (gretl_restriction *rset)
454 {
455 MODEL *pmod = rset->obj;
456 int nc = 0;
457 int col, i, j;
458
459 if (rset->mask == NULL) {
460 nc = pmod->ncoeff;
461 } else {
462 for (i=0; i<pmod->ncoeff; i++) {
463 if (rset->mask[i]) {
464 nc++;
465 }
466 }
467 }
468
469 if (rset_allocate_R_q(rset, nc)) {
470 return E_ALLOC;
471 }
472
473 for (i=0; i<rset->g; i++) {
474 rrow *r = rset->rows[i];
475 double x;
476
477 col = 0;
478 for (j=0; j<pmod->ncoeff; j++) {
479 if (rset->mask[j]) {
480 x = get_restriction_param(r, j);
481 gretl_matrix_set(rset->R, i, col++, x);
482 }
483 }
484 gretl_vector_set(rset->q, i, r->rhs);
485 }
486
487 return 0;
488 }
489
490 /* For the equation-system case, allocate and fill out the
491 R and q matrices. */
492
sys_form_matrices(gretl_restriction * rset)493 static int sys_form_matrices (gretl_restriction *rset)
494 {
495 int nc = system_n_indep_vars(rset->obj);
496 int col, i, j;
497 int err = 0;
498
499 if (rset_allocate_R_q(rset, nc)) {
500 return E_ALLOC;
501 }
502
503 for (i=0; i<rset->g; i++) {
504 rrow *r = rset->rows[i];
505
506 for (j=0; j<r->nterms; j++) {
507 col = get_R_sys_column(rset, i, j);
508 gretl_matrix_set(rset->R, i, col, r->mult[j]);
509 }
510 gretl_vector_set(rset->q, i, r->rhs);
511 }
512
513 err = check_R_matrix(rset->R);
514
515 return err;
516 }
517
518 /* For the VAR case, allocate and fill out the restriction
519 matrices. */
520
var_form_matrices(gretl_restriction * rset)521 static int var_form_matrices (gretl_restriction *rset)
522 {
523 GRETL_VAR *var = rset->obj;
524 int nc = var->neqns * var->ncoeff;
525 int col, i, j;
526 int err = 0;
527
528 if (rset_allocate_R_q(rset, nc)) {
529 return E_ALLOC;
530 }
531
532 for (i=0; i<rset->g; i++) {
533 rrow *r = rset->rows[i];
534
535 for (j=0; j<r->nterms; j++) {
536 col = get_R_var_column(rset, i, j);
537 gretl_matrix_set(rset->R, i, col, r->mult[j]);
538 }
539 gretl_vector_set(rset->q, i, r->rhs);
540 }
541
542 err = check_R_matrix(rset->R);
543
544 return err;
545 }
546
547 /* For the VECM case, allocate and fill out the restriction
548 matrices. */
549
vecm_form_matrices(gretl_restriction * rset)550 static int vecm_form_matrices (gretl_restriction *rset)
551 {
552 int i, j, m, err;
553
554 err = vecm_restriction_check(rset);
555 if (err) {
556 return err;
557 }
558
559 if (rset->vecm & VECM_B) {
560 rset->R = gretl_zero_matrix_new(rset->gb, rset->bcols);
561 rset->q = gretl_zero_matrix_new(rset->gb, 1);
562 if (rset->R == NULL || rset->q == NULL) {
563 return E_ALLOC;
564 }
565 }
566
567 if (rset->vecm & VECM_A) {
568 rset->Ra = gretl_zero_matrix_new(rset->ga, rset->acols);
569 rset->qa = gretl_zero_matrix_new(rset->ga, 1);
570 if (rset->Ra == NULL || rset->qa == NULL) {
571 return E_ALLOC;
572 }
573 }
574
575 /* construct the restriction matrices, beta first then alpha
576 (if both are given) */
577
578 for (m=0; m<2; m++) {
579 gretl_matrix *R = (m == 0)? rset->R : rset->Ra;
580 gretl_matrix *q = (m == 0)? rset->q : rset->qa;
581 char letter = (m == 0)? 'b' : 'a';
582 int col, row = 0;
583
584 if ((letter == 'b' && rset->vecm == VECM_A) ||
585 (letter == 'a' && rset->vecm == VECM_B)) {
586 continue;
587 }
588
589 for (i=0; i<rset->g; i++) {
590 rrow *r = rset->rows[i];
591
592 if (r->letter[0] == letter) {
593 for (j=0; j<r->nterms; j++) {
594 col = get_R_vecm_column(rset, i, j, letter);
595 gretl_matrix_set(R, row, col, r->mult[j]);
596 }
597 gretl_vector_set(q, row, r->rhs);
598 row++;
599 }
600 }
601
602 #if RDEBUG
603 if (letter == 'a') {
604 gretl_matrix_print(R, "Ra, constructed");
605 gretl_matrix_print(q, "qa, constructed");
606 }
607 #endif
608 }
609
610 /* cumulate prior restrictions if needed */
611
612 if (!err && beta_restricted_VECM(rset->obj)) {
613 if (rset->vecm & VECM_B) {
614 err = augment_vecm_restriction(rset, 'b');
615 } else {
616 err = add_old_vecm_restriction(rset, 'b');
617 }
618 }
619
620 if (!err && alpha_restricted_VECM(rset->obj)) {
621 if (rset->vecm & VECM_A) {
622 err = augment_vecm_restriction(rset, 'a');
623 } else {
624 err = add_old_vecm_restriction(rset, 'a');
625 }
626 }
627
628 if (!err && rset->R != NULL) {
629 err = check_R_matrix(rset->R);
630 }
631
632 if (!err && rset->Ra != NULL) {
633 err = check_R_matrix(rset->Ra);
634 }
635
636 #if RDEBUG
637 fprintf(stderr, "vecm_form_matrices: returning %d\n", err);
638 #endif
639
640 return err;
641 }
642
643 /* We were given R and/or q directly by the user: check them
644 for sanity. */
645
test_user_matrices(gretl_restriction * rset)646 static int test_user_matrices (gretl_restriction *rset)
647 {
648 gretl_matrix *R = rset->R;
649 gretl_matrix *q = rset->q;
650 int err = 0;
651
652 if (R == NULL || q == NULL) {
653 /* we didn't get both parts */
654 return E_DATA;
655 }
656
657 if (R->rows != q->rows || q->cols != 1) {
658 /* R and q don't work */
659 return E_NONCONF;
660 }
661
662 if (rset->vecm) {
663 /* this is supposed to be a restriction on beta:
664 test it as such */
665 int nb = gretl_VECM_n_beta(rset->obj);
666 int nbr = nb * gretl_VECM_rank(rset->obj);
667
668 if (R->cols == nb && R->rows <= nb) {
669 /* common restriction on columns of beta */
670 rset->bmulti = 0;
671 } else if (R->cols == nbr && R->rows <= nbr) {
672 ; /* OK, full restriction on beta */
673 } else {
674 gretl_errmsg_set(_("Invalid restriction on beta"));
675 err = E_NONCONF;
676 }
677 if (!err) {
678 rset->bcols = R->cols;
679 }
680 } else {
681 /* the non-VECM case */
682 if (R->rows > rset->gmax) {
683 gretl_errmsg_sprintf(_("Too many restrictions: the maximum "
684 "number is %d"), rset->gmax);
685 err = E_NONCONF;
686 } else if (R->cols != rset->gmax) {
687 gretl_errmsg_sprintf(_("The R matrix must have %d columns"),
688 rset->gmax);
689 err = E_NONCONF;
690 }
691 }
692
693 return err;
694 }
695
696 /* Set up the matrices needed for testing a set of restrictions:
697 general driver function that calls more specific functions
698 depending on the case. */
699
restriction_set_form_matrices(gretl_restriction * rset)700 static int restriction_set_form_matrices (gretl_restriction *rset)
701 {
702 int err = 0;
703
704 if (rset->R != NULL || rset->q != NULL) {
705 err = test_user_matrices(rset);
706 } else if (rset->otype == GRETL_OBJ_EQN) {
707 err = equation_form_matrices(rset);
708 } else if (rset->otype == GRETL_OBJ_VAR) {
709 if (gretl_VECM_rank(rset->obj) > 0) {
710 err = vecm_form_matrices(rset);
711 } else {
712 err = var_form_matrices(rset);
713 }
714 } else if (rset->g > 0) {
715 /* mutivariate system */
716 err = sys_form_matrices(rset);
717 }
718
719 #if RDEBUG
720 gretl_matrix_print(rset->R, "R");
721 gretl_matrix_print(rset->q, "q");
722 #endif
723
724 return err;
725 }
726
727 /* Make a mask with 1s in positions in the array of coeffs where
728 a coeff is referenced in one or more restrictions, 0s otherwise.
729 We do this only for single-equation restriction sets, and
730 only when the restrictions are given in row-wise form.
731 */
732
restriction_set_make_mask(gretl_restriction * rset)733 static int restriction_set_make_mask (gretl_restriction *rset)
734 {
735 MODEL *pmod;
736 rrow *r;
737 int i, j;
738
739 if (rset->otype != GRETL_OBJ_EQN || rset->obj == NULL) {
740 return E_DATA;
741 }
742
743 pmod = rset->obj;
744 rset->mask = calloc(pmod->ncoeff, 1);
745
746 if (rset->mask == NULL) {
747 destroy_restriction_set(rset);
748 return E_ALLOC;
749 }
750
751 for (i=0; i<rset->g; i++) {
752 r = rset->rows[i];
753 for (j=0; j<r->nterms; j++) {
754 rset->mask[r->bnum[j]] = 1;
755 }
756 }
757
758 return 0;
759 }
760
count_ops(const char * p)761 static int count_ops (const char *p)
762 {
763 int n = 0 ;
764
765 while (*p) {
766 if (*p == '+' || *p == '-') n++;
767 if (*p == '=') break;
768 p++;
769 }
770
771 return n;
772 }
773
coeff_index_from_scalar(const char * s)774 static int coeff_index_from_scalar (const char *s)
775 {
776 int err = 0;
777 double x = gretl_scalar_get_value(s, &err);
778
779 return (err || na(x))? -1 : (int) x;
780 }
781
782 /* Try to retrieve the coefficient number in a model to be
783 restricted, based on the parameter name, or possibly the
784 name of a scalar variable.
785 */
786
787 static int
bnum_from_name(gretl_restriction * r,const DATASET * dset,const char * s)788 bnum_from_name (gretl_restriction *r, const DATASET *dset,
789 const char *s)
790 {
791 int k = -1;
792
793 if (dset == NULL || r->obj == NULL) {
794 ;
795 } else if (r->otype == GRETL_OBJ_VAR) {
796 GRETL_VAR *var = r->obj;
797 int i, v;
798
799 for (i=1; i<=var->ylist[0]; i++) {
800 v = var->ylist[i];
801 if (v > 0 && v < dset->v && !strcmp(dset->varname[v], s)) {
802 k = i;
803 break;
804 }
805 }
806 } else if (r->otype == GRETL_OBJ_EQN) {
807 const MODEL *pmod = r->obj;
808
809 k = gretl_model_get_param_number(pmod, dset, s);
810 if (k >= 0) {
811 /* convert to 1-based for compatibility with numbers read
812 directly: the index will be converted to 0-base below
813 */
814 k++;
815 }
816 }
817
818 if (k < 0) {
819 /* try for a scalar variable */
820 k = coeff_index_from_scalar(s);
821 }
822 if (k < 0) {
823 gretl_errmsg_sprintf(_("%s: couldn't interpret as coefficient number"), s);
824 }
825
826 #if RDEBUG
827 fprintf(stderr, "bnum_from_name: %s -> bnum = %d\n", s, k);
828 #endif
829
830 return k;
831 }
832
833 /* Pick apart strings of the form "b[X]" or "b[X,Y]". If the ",Y" is
834 present the "X" element must be an equation number, and the "Y" may
835 be a coefficient number or the name of a variable. If the ",Y" is
836 not present, "X" may be a coefficient number or the name of a
837 variable.
838
839 Such terms may be preceded by a multiplier, but they should not
840 be followed by anything, so we flag an error if the closing "]"
841 is not the end of the string.
842 */
843
pick_apart(gretl_restriction * r,const char * s,int * eq,int * bnum,const DATASET * dset)844 static int pick_apart (gretl_restriction *r, const char *s,
845 int *eq, int *bnum,
846 const DATASET *dset)
847 {
848 int parens_ok = 0;
849 const char *junk;
850 char s1[VNAMELEN] = {0};
851 char s2[VNAMELEN] = {0};
852 char *targ = s1;
853 char rdelim;
854 int i, j, k;
855 int err = 0;
856
857 #if RDEBUG
858 fprintf(stderr, "pick_apart: looking at '%s'\n", s);
859 #endif
860
861 /* register then skip past the opening delimiter */
862 rdelim = (*s == '[')? ']' : ')';
863 s++;
864
865 *eq = *bnum = -1;
866
867 k = gretl_charpos(rdelim, s);
868 if (k <= 0 || k > 2*(VNAMELEN-1)) {
869 return E_PARSE;
870 }
871
872 junk = s + k + 1;
873 while (isspace(*junk)) junk++;
874 if (*junk != '\0') {
875 fprintf(stderr, "Got trailing junk '%s'\n", junk);
876 return E_PARSE;
877 }
878
879 j = 0;
880 for (i=0; i<k; i++) {
881 if (s[i] == ',') {
882 targ = s2;
883 j = 0;
884 } else if (!isspace(s[i])) {
885 if (j == VNAMELEN-1) {
886 return E_PARSE;
887 }
888 targ[j++] = s[i];
889 }
890 }
891
892 #if RDEBUG
893 fprintf(stderr, " s1 = '%s', s2 = '%s'\n", s1, s2);
894 #endif
895
896 if (targ == s2) {
897 /* got a comma separator: [eqn,bnum] */
898 *eq = positive_int_from_string(s1);
899 if (*eq <= 0) {
900 err = E_PARSE;
901 }
902 if (isdigit(*s2)) {
903 *bnum = positive_int_from_string(s2);
904 } else {
905 *bnum = bnum_from_name(r, dset, s2);
906 }
907 } else {
908 /* only one field: [bnum] */
909 *eq = R_UNSPEC;
910 if (isdigit(*s1)) {
911 *bnum = positive_int_from_string(s1);
912 } else {
913 parens_ok = 1;
914 *bnum = bnum_from_name(r, dset, s1);
915 }
916 }
917
918 if (!err && rdelim == ')' && !parens_ok) {
919 err = E_PARSE;
920 }
921
922 return err;
923 }
924
bbit_trailing_garbage(const char * s)925 static int bbit_trailing_garbage (const char *s)
926 {
927 if (*s) {
928 while (isspace(*s)) s++;
929 if (*s) {
930 if (*s == '/') {
931 gretl_errmsg_sprintf(_("%s: division not allowed here"), s-1);
932 }
933 return 1;
934 }
935 }
936
937 return 0;
938 }
939
940 /* We got a leading letter ('a' or 'b') and now we parse
941 what follows. The simplest case is just a number, giving
942 e.g., "b1". Then we assess the validity of what we've
943 found.
944 */
945
parse_b_bit(gretl_restriction * r,const char * s,int * eq,int * bnum,const DATASET * dset)946 static int parse_b_bit (gretl_restriction *r, const char *s,
947 int *eq, int *bnum,
948 const DATASET *dset)
949 {
950 int err = E_PARSE;
951
952 if (isdigit((unsigned char) *s)) {
953 char *test;
954
955 *bnum = strtol(s, &test, 10);
956 if (bbit_trailing_garbage(test)) {
957 return err;
958 }
959 if (r->otype == GRETL_OBJ_VAR) {
960 *eq = R_UNSPEC;
961 }
962 err = 0;
963 } else if (*s == '[' || *s == '(') {
964 err = pick_apart(r, s, eq, bnum, dset);
965 }
966
967 if (*bnum < 1) {
968 if (!gretl_errmsg_is_set() && *bnum >= 0) {
969 gretl_errmsg_sprintf(_("Coefficient number (%d) is out of range"),
970 *bnum);
971 }
972 err = 1;
973 } else {
974 *bnum -= 1; /* convert to zero base */
975 }
976
977 if (*eq == R_UNSPEC) {
978 /* didn't get a specific equation number */
979 if (r->otype == GRETL_OBJ_EQN) {
980 *eq = 0;
981 } else if (r->otype != GRETL_OBJ_VAR) {
982 err = E_PARSE;
983 }
984 } else if (*eq == 0) {
985 gretl_errmsg_sprintf(_("Equation number (%d) is out of range"),
986 *eq);
987 err = 1;
988 } else {
989 *eq -= 1; /* convert to zero base */
990 }
991
992 return err;
993 }
994
r_get_scalar(const char * s,double * x)995 static int r_get_scalar (const char *s, double *x)
996 {
997 if (sscanf(s, "%lf", x)) {
998 return 1;
999 } else {
1000 int n;
1001
1002 s += strspn(s, " ");
1003 n = gretl_namechar_spn(s);
1004
1005 if (n > 0) {
1006 char tmp[VNAMELEN];
1007
1008 *tmp = '\0';
1009 strncat(tmp, s, n);
1010 if (gretl_is_scalar(tmp)) {
1011 *x = gretl_scalar_get_value(tmp, NULL);
1012 if (!na(*x)) {
1013 return 1;
1014 }
1015 }
1016 }
1017 }
1018
1019 return 0;
1020 }
1021
get_named_param(const char * s,MODEL * pmod,int * bnum)1022 static int get_named_param (const char *s, MODEL *pmod,
1023 int *bnum)
1024 {
1025 if (NONLIST_MODEL(pmod->ci) && pmod->params != NULL) {
1026 char test[VNAMELEN];
1027 int i;
1028
1029 *test = '\0';
1030 strncat(test, s, VNAMELEN - 1);
1031 tailstrip(test);
1032
1033 for (i=0; i<pmod->nparams; i++) {
1034 if (!strcmp(test, pmod->params[i])) {
1035 *bnum = i;
1036 return 1;
1037 }
1038 }
1039 }
1040
1041 return 0;
1042 }
1043
starts_b(const char * s,gretl_restriction * r)1044 static int starts_b (const char *s, gretl_restriction *r)
1045 {
1046 if (*s == 'b' || (r->vecm && *s == 'a')) {
1047 char c = *(s+1);
1048
1049 return isdigit(c) || c == '[' || c == '(';
1050 }
1051
1052 return 0;
1053 }
1054
1055 #define could_be_param(r,s) (r->otype == GRETL_OBJ_EQN && \
1056 !isdigit(*s))
1057
1058 static int
parse_coeff_chunk(gretl_restriction * r,const char * s,double * x,int * eq,int * bnum,char * letter,const DATASET * dset)1059 parse_coeff_chunk (gretl_restriction *r, const char *s, double *x,
1060 int *eq, int *bnum, char *letter,
1061 const DATASET *dset)
1062 {
1063 const char *s0 = s;
1064 int done = 0;
1065 int err = E_PARSE;
1066
1067 *eq = 1;
1068
1069 while (isspace((unsigned char) *s)) s++;
1070
1071 #if RDEBUG
1072 fprintf(stderr, "parse_coeff_chunk: s='%s'\n", s);
1073 #endif
1074
1075 if (starts_b(s, r)) {
1076 /* something like "b[" or "b2" */
1077 *letter = *s;
1078 s++;
1079 #if RDEBUG
1080 fprintf(stderr, " first branch: s='%s'\n", s);
1081 #endif
1082 err = parse_b_bit(r, s, eq, bnum, dset);
1083 *x = 1.0;
1084 done = 1;
1085 } else if (could_be_param(r, s)) {
1086 /* e.g. "beta" for MLE or GMM */
1087 done = get_named_param(s, r->obj, bnum);
1088 if (done) {
1089 *letter = 'b';
1090 *x = 1.0;
1091 *eq = 0;
1092 err = 0;
1093 }
1094 }
1095
1096 if (!done && r_get_scalar(s, x)) {
1097 s += strspn(s, " ");
1098 s += strcspn(s, " *");
1099 s += strspn(s, " *");
1100 if (starts_b(s, r)) {
1101 *letter = *s;
1102 s++;
1103 err = parse_b_bit(r, s, eq, bnum, dset);
1104 } else if (*s == '\0') {
1105 /* plain numeric value, saved in *x */
1106 *bnum = 0;
1107 *letter = '\0';
1108 err = 0;
1109 }
1110 }
1111
1112 if (err) {
1113 gretl_errmsg_sprintf(_("parse error in '%s'\n"), s0);
1114 }
1115
1116 #if RDEBUG
1117 fprintf(stderr, "parse_coeff_chunk: x=%g, eq=%d, bnum=%d, letter=%c\n",
1118 *x, *eq, *bnum, *letter);
1119 #endif
1120
1121 return err;
1122 }
1123
destroy_restriction(rrow * r)1124 static void destroy_restriction (rrow *r)
1125 {
1126 if (r != NULL) {
1127 free(r->mult);
1128 free(r->eq);
1129 free(r->bnum);
1130 free(r->letter);
1131 free(r);
1132 }
1133 }
1134
destroy_restriction_set(gretl_restriction * rset)1135 void destroy_restriction_set (gretl_restriction *rset)
1136 {
1137 if (rset->rows != NULL) {
1138 int i;
1139
1140 for (i=0; i<rset->g; i++) {
1141 destroy_restriction(rset->rows[i]);
1142 }
1143 free(rset->rows);
1144 }
1145
1146 free(rset->mask);
1147
1148 #if RDEBUG
1149 fprintf(stderr, "destroy_restriction_set: R at %p, q at %p\n",
1150 (void *) rset->R, (void *) rset->q);
1151 fprintf(stderr, " Ra at %p, qa at %p\n",
1152 (void *) rset->Ra, (void *) rset->qa);
1153 #endif
1154
1155 gretl_matrix_free(rset->R);
1156 gretl_matrix_free(rset->q);
1157 gretl_matrix_free(rset->Ra);
1158 gretl_matrix_free(rset->qa);
1159
1160 free(rset->rfunc);
1161
1162 free(rset);
1163 }
1164
restriction_new(int n,int multi,int vecm)1165 static rrow *restriction_new (int n, int multi, int vecm)
1166 {
1167 rrow *r;
1168 int i;
1169
1170 r = malloc(sizeof *r);
1171 if (r == NULL) {
1172 return NULL;
1173 }
1174
1175 r->mult = NULL;
1176 r->eq = NULL;
1177 r->bnum = NULL;
1178 r->letter = NULL;
1179
1180 r->mult = malloc(n * sizeof *r->mult);
1181 r->bnum = malloc(n * sizeof *r->bnum);
1182 if (r->mult == NULL || r->bnum == NULL) {
1183 destroy_restriction(r);
1184 return NULL;
1185 }
1186
1187 for (i=0; i<n; i++) {
1188 r->mult[i] = 0.0;
1189 r->bnum[i] = 0;
1190 }
1191
1192 if (multi) {
1193 r->eq = malloc(n * sizeof *r->eq);
1194 if (r->eq == NULL) {
1195 destroy_restriction(r);
1196 return NULL;
1197 }
1198 }
1199
1200 if (vecm) {
1201 r->letter = malloc(n * sizeof *r->letter);
1202 if (r->letter == NULL) {
1203 destroy_restriction(r);
1204 return NULL;
1205 }
1206 }
1207
1208 r->nterms = n;
1209 r->rhs = 0.0;
1210
1211 return r;
1212 }
1213
restriction_set_add_row(gretl_restriction * rset,int n_terms)1214 static rrow *restriction_set_add_row (gretl_restriction *rset,
1215 int n_terms)
1216 {
1217 rrow **rows = NULL;
1218 int n = rset->g;
1219
1220 rows = realloc(rset->rows, (n + 1) * sizeof *rows);
1221 if (rows == NULL) {
1222 return NULL;
1223 }
1224
1225 rset->rows = rows;
1226
1227 rset->rows[n] = restriction_new(n_terms, rset->bmulti, rset->vecm);
1228 if (rset->rows[n] == NULL) {
1229 return NULL;
1230 }
1231
1232 rset->g += 1;
1233
1234 return rset->rows[n];
1235 }
1236
print_mult(double mult,int first,PRN * prn)1237 static void print_mult (double mult, int first, PRN *prn)
1238 {
1239 if (mult == 1.0) {
1240 if (!first) pputs(prn, " + ");
1241 } else if (mult == -1.0) {
1242 if (first) pputs(prn, "-");
1243 else pputs(prn, " - ");
1244 } else if (mult > 0.0) {
1245 if (first) pprintf(prn, "%g*", mult);
1246 else pprintf(prn, " + %g*", mult);
1247 } else if (mult < 0.0) {
1248 if (first) pprintf(prn, "%g*", mult);
1249 else pprintf(prn, " - %g*", fabs(mult));
1250 }
1251 }
1252
print_restriction(const gretl_restriction * rset,int i,const DATASET * dset,PRN * prn)1253 static void print_restriction (const gretl_restriction *rset,
1254 int i, const DATASET *dset,
1255 PRN *prn)
1256 {
1257 const rrow *r = rset->rows[i];
1258 char letter, vname[VNAMELEN];
1259 int j, k;
1260
1261 for (j=0; j<r->nterms; j++) {
1262 letter = (r->letter != NULL)? r->letter[j] : 'b';
1263 k = r->bnum[j];
1264 print_mult(r->mult[j], j == 0, prn);
1265 if (targ_specified(rset, i, j, letter)) {
1266 pprintf(prn, "%c[%d,%d]", letter, r->eq[j] + 1, k + 1);
1267 } else if (rset->otype == GRETL_OBJ_VAR) {
1268 pprintf(prn, "%c[%d]", letter, k + 1);
1269 } else {
1270 MODEL *pmod = rset->obj;
1271
1272 gretl_model_get_param_name(pmod, dset, k, vname);
1273 if (NONLIST_MODEL(pmod->ci)) {
1274 pputs(prn, vname);
1275 } else {
1276 pprintf(prn, "b[%s]", vname);
1277 }
1278 }
1279 }
1280
1281 pprintf(prn, " = %g\n", r->rhs);
1282 }
1283
1284 /* Pretty-print the set of restrictions. We do this only if
1285 the restrictions were given row-wise by the user.
1286 */
1287
print_restriction_set(const gretl_restriction * rset,const DATASET * dset,gretlopt opt,PRN * prn)1288 static void print_restriction_set (const gretl_restriction *rset,
1289 const DATASET *dset,
1290 gretlopt opt, PRN *prn)
1291 {
1292 int i;
1293
1294 if (prn == NULL) {
1295 return;
1296 }
1297
1298 if (rset->g > 1) {
1299 pputs(prn, _("Restriction set"));
1300 } else {
1301 pprintf(prn, "%s:", _("Restriction"));
1302 }
1303 pputc(prn, '\n');
1304
1305 for (i=0; i<rset->g; i++) {
1306 if (rset->g > 1) {
1307 pprintf(prn, " %d: ", i + 1);
1308 } else {
1309 pputc(prn, ' ');
1310 }
1311 print_restriction(rset, i, dset, prn);
1312 }
1313
1314 if (opt & OPT_V) {
1315 if (rset->R != NULL || rset->q != NULL ||
1316 rset->Ra != NULL || rset->qa != NULL) {
1317 pputc(prn, '\n');
1318 }
1319 if (rset->R != NULL) {
1320 gretl_matrix_print_to_prn(rset->R, "R", prn);
1321 }
1322 if (rset->q != NULL) {
1323 gretl_matrix_print_to_prn(rset->q, "q", prn);
1324 }
1325 if (rset->Ra != NULL) {
1326 gretl_matrix_print_to_prn(rset->Ra, "R(alpha)", prn);
1327 }
1328 if (rset->qa != NULL) {
1329 gretl_matrix_print_to_prn(rset->qa, "q(alpha)", prn);
1330 }
1331 }
1332 }
1333
print_restriction_from_matrices(const gretl_matrix * R,const gretl_matrix * q,char letter,int npar,PRN * prn)1334 void print_restriction_from_matrices (const gretl_matrix *R,
1335 const gretl_matrix *q,
1336 char letter, int npar,
1337 PRN *prn)
1338 {
1339 double x;
1340 int eqn, coeff, started;
1341 int i, j;
1342
1343 for (i=0; i<R->rows; i++) {
1344 started = 0;
1345 coeff = 1;
1346 eqn = (R->cols > npar)? 1 : 0;
1347 for (j=0; j<R->cols; j++) {
1348 x = gretl_matrix_get(R, i, j);
1349 if (x != 0.0) {
1350 if (!started) {
1351 pputs(prn, " ");
1352 }
1353 if (x == 1.0) {
1354 if (started) {
1355 pputs(prn, " + ");
1356 }
1357 } else if (x == -1.0) {
1358 if (started) {
1359 pputs(prn, " - ");
1360 } else {
1361 pputc(prn, '-');
1362 }
1363 } else if (x > 0.0) {
1364 if (started) {
1365 pprintf(prn, " + %g*", x);
1366 } else {
1367 pprintf(prn, "%g*", x);
1368 }
1369 } else if (x < 0.0) {
1370 if (started) {
1371 pprintf(prn, " - %g*", -x);
1372 } else {
1373 pprintf(prn, "%g*", x);
1374 }
1375 }
1376 if (eqn > 0) {
1377 pprintf(prn, "%c[%d,%d]", letter, eqn, coeff);
1378 } else {
1379 pprintf(prn, "%c%d", letter, coeff);
1380 }
1381 started = 1;
1382 }
1383 if ((j + 1) % npar == 0) {
1384 eqn++;
1385 coeff = 1;
1386 } else {
1387 coeff++;
1388 }
1389 }
1390 pprintf(prn, " = %g\n", (q == NULL)? 0.0 : q->val[i]);
1391 }
1392 }
1393
1394 static int
add_term_to_restriction(rrow * r,double mult,int eq,int bnum,char letter,int i)1395 add_term_to_restriction (rrow *r, double mult, int eq, int bnum,
1396 char letter, int i)
1397 {
1398 int j;
1399
1400 for (j=0; j<i; j++) {
1401 if (r->letter != NULL && r->letter[j] != letter) {
1402 continue;
1403 }
1404 if (bnum == r->bnum[j] && (r->eq == NULL || eq == r->eq[j])) {
1405 /* additional reference to a previously referenced coeff */
1406 r->mult[j] += mult;
1407 r->nterms -= 1;
1408 return 0;
1409 }
1410 }
1411
1412 r->mult[i] = mult;
1413 r->bnum[i] = bnum;
1414
1415 if (r->eq != NULL) {
1416 r->eq[i] = eq;
1417 }
1418
1419 if (r->letter != NULL) {
1420 r->letter[i] = letter;
1421 }
1422
1423 return 0;
1424 }
1425
restriction_set_new(void * ptr,GretlObjType type,gretlopt opt)1426 static gretl_restriction *restriction_set_new (void *ptr,
1427 GretlObjType type,
1428 gretlopt opt)
1429 {
1430 gretl_restriction *rset;
1431
1432 rset = malloc(sizeof *rset);
1433 if (rset == NULL) return NULL;
1434
1435 rset->obj = ptr;
1436 rset->otype = type;
1437 rset->opt = opt;
1438
1439 rset->test = NADBL;
1440 rset->pval = NADBL;
1441 rset->lnl = NADBL;
1442 rset->bsum = NADBL;
1443 rset->bse = NADBL;
1444
1445 rset->g = rset->gmax = 0;
1446 rset->gb = rset->ga = 0;
1447 rset->R = NULL;
1448 rset->q = NULL;
1449 rset->Ra = NULL;
1450 rset->qa = NULL;
1451 rset->mask = NULL;
1452 rset->rows = NULL;
1453
1454 rset->rfunc = NULL;
1455
1456 rset->bmulti = 0;
1457 rset->amulti = 0;
1458 rset->bcols = 0;
1459 rset->acols = 0;
1460 rset->vecm = 0;
1461 rset->code = GRETL_STAT_NONE;
1462
1463 if (rset->otype == GRETL_OBJ_EQN) {
1464 MODEL *pmod = ptr;
1465
1466 rset->gmax = pmod->ncoeff;
1467 } else if (rset->otype == GRETL_OBJ_SYS) {
1468 rset->gmax = system_n_indep_vars(ptr);
1469 rset->bmulti = 1;
1470 } else if (rset->otype == GRETL_OBJ_VAR) {
1471 GRETL_VAR *var = ptr;
1472 int r = gretl_VECM_rank(var);
1473
1474 if (r == 0) {
1475 /* a plain VAR */
1476 rset->gmax = var->neqns * var->ncoeff;
1477 rset->bmulti = 1;
1478 } else {
1479 if (r > 1) {
1480 rset->bmulti = 1;
1481 rset->amulti = 1;
1482 }
1483 rset->vecm = VECM_B;
1484 rset->gmax = gretl_VECM_n_beta(var) * r;
1485 rset->gmax += gretl_VECM_n_alpha(var) * r;
1486 }
1487 }
1488
1489 #if RDEBUG
1490 fprintf(stderr, "restriction_set_new: rset->gmax = %d\n", rset->gmax);
1491 #endif
1492
1493 return rset;
1494 }
1495
rset_from_VECM(GRETL_VAR * var,int * err)1496 gretl_restriction *rset_from_VECM (GRETL_VAR *var, int *err)
1497 {
1498 const gretl_matrix *R = gretl_VECM_R_matrix(var);
1499 const gretl_matrix *q = gretl_VECM_q_matrix(var);
1500 const gretl_matrix *Ra = gretl_VECM_Ra_matrix(var);
1501 const gretl_matrix *qa = gretl_VECM_qa_matrix(var);
1502 gretl_restriction *rset;
1503
1504 if (R == NULL && q == NULL && Ra == NULL && qa == NULL) {
1505 /* no restriction in place */
1506 return NULL;
1507 }
1508
1509 rset = calloc(1, sizeof *rset);
1510 if (rset == NULL) {
1511 *err = E_ALLOC;
1512 return NULL;
1513 }
1514
1515 rset->obj = var;
1516 rset->otype = GRETL_OBJ_VAR;
1517 rset->opt = OPT_B;
1518
1519 rset->R = (gretl_matrix *) R;
1520 rset->q = (gretl_matrix *) q;
1521 rset->Ra = (gretl_matrix *) Ra;
1522 rset->qa = (gretl_matrix *) qa;
1523
1524 return rset;
1525 }
1526
list_get_bmax(const int * list)1527 static int list_get_bmax (const int *list)
1528 {
1529 int i, bmax = 0;
1530
1531 for (i=2; i<=list[0]; i++) {
1532 if (list[i] == LISTSEP) {
1533 break;
1534 }
1535 bmax++;
1536 }
1537
1538 return bmax;
1539 }
1540
1541 /* check that the coefficients referenced in a restriction are
1542 within bounds, relative to the equation or system that is
1543 to be restricted */
1544
bnum_out_of_bounds(const gretl_restriction * rset,int eq,int bnum,char letter)1545 static int bnum_out_of_bounds (const gretl_restriction *rset,
1546 int eq, int bnum, char letter)
1547 {
1548 int ret = 1;
1549
1550 if (rset->otype == GRETL_OBJ_VAR) {
1551 GRETL_VAR *var = rset->obj;
1552 int r = gretl_VECM_rank(var);
1553
1554 if (r == 0) {
1555 if (eq >= var->neqns) {
1556 gretl_errmsg_sprintf(_("Equation number (%d) is out of range"),
1557 eq + 1);
1558 } else if (bnum >= var->ncoeff) {
1559 gretl_errmsg_sprintf(_("Coefficient number (%d) is out of range"),
1560 bnum + 1);
1561 } else {
1562 ret = 0;
1563 }
1564 } else {
1565 /* note: @eq and @bnum are transposed here */
1566 if (eq >= gretl_VECM_rank(var)) {
1567 gretl_errmsg_sprintf(_("Coefficient number (%d) is out of range"),
1568 eq + 1);
1569 } else if ((letter == 'b' && bnum >= gretl_VECM_n_beta(var)) ||
1570 (letter == 'a' && bnum >= gretl_VECM_n_alpha(var))) {
1571 gretl_errmsg_sprintf(_("Equation number (%d) is out of range"),
1572 bnum + 1);
1573 } else {
1574 ret = 0;
1575 }
1576 }
1577 } else if (rset->otype == GRETL_OBJ_SYS) {
1578 equation_system *sys = rset->obj;
1579 const int *list = system_get_list(sys, eq);
1580
1581 if (list == NULL) {
1582 gretl_errmsg_sprintf(_("Equation number (%d) is out of range"),
1583 eq + 1);
1584 } else {
1585 int bmax = list_get_bmax(list);
1586
1587 if (bnum >= bmax) {
1588 gretl_errmsg_sprintf(_("Coefficient number (%d) out of range "
1589 "for equation %d"), bnum + 1, eq + 1);
1590 } else {
1591 ret = 0;
1592 }
1593 }
1594 } else {
1595 MODEL *pmod = rset->obj;
1596
1597 if (eq > 0) {
1598 gretl_errmsg_sprintf(_("Equation number (%d) is out of range"),
1599 eq + 1);
1600 } else if (bnum >= pmod->ncoeff || bnum < 0) {
1601 if (!gretl_errmsg_is_set()) {
1602 gretl_errmsg_sprintf(_("Coefficient number (%d) is out of range"),
1603 bnum + 1);
1604 }
1605 } else {
1606 ret = 0;
1607 }
1608 }
1609
1610 return ret;
1611 }
1612
1613 /* Parse a matrix specification for a restriction set. This should
1614 take the form "R = <matrix>" or "q = <matrix>", where <matrix> is
1615 the name of a pre-defined matrix. On this approach we want exactly
1616 one R definition and exactly one q definition, and no other sort
1617 of specification can be given.
1618 */
1619
read_matrix_line(const char * s,gretl_restriction * rset)1620 static int read_matrix_line (const char *s, gretl_restriction *rset)
1621 {
1622 const gretl_matrix *m;
1623 char mname[VNAMELEN];
1624 char job = *s;
1625 int err = 0;
1626
1627 if (rset->rows != NULL) {
1628 /* we already have "row-wise" restrictions */
1629 return E_PARSE;
1630 } else if (job == 'R' && rset->R != NULL) {
1631 /* duplicate entry for R matrix */
1632 return E_PARSE;
1633 } else if (job == 'q' && rset->q != NULL) {
1634 /* duplicate entry for q matrix */
1635 return E_PARSE;
1636 }
1637
1638 s++;
1639 while (isspace((unsigned char) *s)) s++;
1640 if (*s != '=') {
1641 return E_PARSE;
1642 }
1643
1644 s++;
1645 while (isspace((unsigned char) *s)) s++;
1646 if (gretl_scan_varname(s, mname) != 1) {
1647 return E_PARSE;
1648 }
1649
1650 m = get_matrix_by_name(mname);
1651 if (m == NULL) {
1652 return E_UNKVAR;
1653 }
1654
1655 if (job == 'R') {
1656 rset->R = gretl_matrix_copy(m);
1657 if (rset->R == NULL) {
1658 err = E_ALLOC;
1659 }
1660 } else if (job == 'q') {
1661 rset->q = gretl_matrix_copy(m);
1662 if (rset->q == NULL) {
1663 err = E_ALLOC;
1664 }
1665 }
1666
1667 return err;
1668 }
1669
1670 /* Try to read the name of a function to be called for testing
1671 a nonlinear restriction, as in "rfunc = <fname>". This is
1672 admissable only if it's not mixed with any other sort
1673 of restriction.
1674 */
1675
read_fncall_line(const char * s,gretl_restriction * rset)1676 static int read_fncall_line (const char *s, gretl_restriction *rset)
1677 {
1678 char fname[FN_NAMELEN];
1679
1680 if (rset->rows != NULL || rset->R != NULL || rset->q != NULL) {
1681 /* other stuff is in the way! */
1682 return E_PARSE;
1683 }
1684
1685 s += strspn(s, " ");
1686
1687 if (*s != '=') {
1688 return E_PARSE;
1689 }
1690
1691 s++;
1692 s += strspn(s, " ");
1693
1694 if (sscanf(s, "%31s", fname) != 1) {
1695 return E_PARSE;
1696 }
1697
1698 s += strlen(fname);
1699
1700 if (!string_is_blank(s)) {
1701 /* there should be no trailing junk */
1702 return E_PARSE;
1703 }
1704
1705 rset->rfunc = gretl_strdup(fname);
1706
1707 return 0;
1708 }
1709
1710 /* We're going to count operators to figure out how many
1711 terms we have in a restriction, so we'd better check
1712 first that we don't have any spurious pile-ups of
1713 operators (that is, '-', '+' and '*'), none of which
1714 can be repeated consecutively, or separated only by
1715 spaces.
1716 */
1717
1718 #define IS_ROP(c) (c=='-' || c == '+' || c == '*')
1719
restriction_op_check(const char * s)1720 static int restriction_op_check (const char *s)
1721 {
1722 int err = 0;
1723
1724 while (*s && !err) {
1725 if (IS_ROP(*s)) {
1726 s++;
1727 while (*s && !err) {
1728 if (*s == ' ') {
1729 ; /* OK so far */
1730 } else if (!(IS_ROP(*s))) {
1731 break;
1732 } else {
1733 err = E_PARSE;
1734 }
1735 s++;
1736 }
1737 }
1738 s++;
1739 }
1740
1741 return err;
1742 }
1743
parse_restriction_row(gretl_restriction * rset,const char * s,const DATASET * dset)1744 static int parse_restriction_row (gretl_restriction *rset,
1745 const char *s,
1746 const DATASET *dset)
1747 {
1748 rrow *row;
1749 int nt, sgn = 1;
1750 int i, j;
1751 int err = 0;
1752
1753 if (rset->R != NULL || rset->q != NULL) {
1754 /* can't mix row-wise spec with matrix spec */
1755 return E_PARSE;
1756 }
1757
1758 err = restriction_op_check(s);
1759 if (err) {
1760 return err;
1761 }
1762
1763 if (*s == '+' || *s == '-') {
1764 sgn = (*s == '+')? 1 : -1;
1765 s++;
1766 }
1767
1768 nt = 1 + count_ops(s);
1769
1770 #if RDEBUG
1771 fprintf(stderr, "restriction line: assuming %d terms\n", nt);
1772 #endif
1773
1774 row = restriction_set_add_row(rset, nt);
1775 if (row == NULL) {
1776 return E_ALLOC;
1777 }
1778
1779 j = 0;
1780
1781 for (i=0; i<nt; i++) {
1782 /* work on the left-hand side */
1783 char chunk[36];
1784 int len, bnum = 1, eq = 1;
1785 int numeric = 0;
1786 char letter = 'b';
1787 double mult;
1788
1789 len = strcspn(s, "+-=");
1790 if (len >= sizeof chunk) {
1791 err = 1;
1792 break;
1793 }
1794
1795 *chunk = '\0';
1796 strncat(chunk, s, len);
1797 s += len;
1798
1799 #if RDEBUG
1800 fprintf(stderr, " working on chunk %d, '%s'\n", i, chunk);
1801 #endif
1802
1803 err = parse_coeff_chunk(rset, chunk, &mult, &eq, &bnum,
1804 &letter, dset);
1805 if (err) {
1806 break;
1807 } else if (letter == '\0') {
1808 numeric = 1;
1809 } else if (bnum_out_of_bounds(rset, eq, bnum, letter)) {
1810 err = E_DATA;
1811 break;
1812 }
1813
1814 mult *= sgn;
1815
1816 if (numeric) {
1817 /* got a numeric constant */
1818 row->rhs -= mult;
1819 row->nterms -= 1;
1820 } else {
1821 add_term_to_restriction(row, mult, eq, bnum, letter, j++);
1822 }
1823
1824 if (*s == '+') {
1825 sgn = 1.0;
1826 s++;
1827 } else if (*s == '-') {
1828 sgn = -1.0;
1829 s++;
1830 }
1831 }
1832
1833 if (!err) {
1834 /* work on the right-hand side */
1835 double rhs = 0.0;
1836
1837 s += strspn(s, " ");
1838
1839 if (*s != '=') {
1840 err = E_PARSE;
1841 } else {
1842 /* we should find here a scalar constant or a
1843 variable of scalar type */
1844 rhs = generate_scalar(s+1, NULL, &err);
1845 if (!err) {
1846 if (na(rhs)) {
1847 err = E_DATA;
1848 } else {
1849 row->rhs += rhs;
1850 }
1851 }
1852 }
1853 }
1854
1855 #if RDEBUG
1856 fprintf(stderr, "restriction row: nterms = %d, rhs = %g\n",
1857 row->nterms, row->rhs);
1858 #endif
1859
1860 if (!err && row->nterms == 0) {
1861 gretl_errmsg_set(_("Invalid restriction"));
1862 err = E_PARSE;
1863 }
1864
1865 return err;
1866 }
1867
parse_strings_line(const char * line,gretl_restriction * rset,const DATASET * dset)1868 static int parse_strings_line (const char *line,
1869 gretl_restriction *rset,
1870 const DATASET *dset)
1871 {
1872 char aname[VNAMELEN];
1873 int err = 0;
1874
1875 if (sscanf(line, "%31s", aname) != 1) {
1876 err = E_PARSE;
1877 } else {
1878 gretl_array *a = get_strings_array_by_name(aname);
1879
1880 if (a == NULL) {
1881 err = E_INVARG;
1882 } else {
1883 char **S;
1884 int i, n;
1885
1886 S = gretl_array_get_strings(a, &n);
1887 for (i=0; i<n && !err; i++) {
1888 err = real_restriction_set_parse_line(rset, S[i], dset, 0);
1889 }
1890 }
1891 }
1892
1893 return err;
1894 }
1895
1896 static int
real_restriction_set_parse_line(gretl_restriction * rset,const char * line,const DATASET * dset,int first)1897 real_restriction_set_parse_line (gretl_restriction *rset,
1898 const char *line,
1899 const DATASET *dset,
1900 int first)
1901 {
1902 const char *s = line;
1903 int err = 0;
1904
1905 gretl_error_clear();
1906
1907 #if RDEBUG
1908 fprintf(stderr, "parse restriction line: got '%s'\n", line);
1909 #endif
1910
1911 /* if we have a function name specified already, nothing else
1912 should be supplied */
1913 if (rset->rfunc != NULL) {
1914 destroy_restriction_set(rset);
1915 return E_PARSE;
1916 }
1917
1918 if (!strncmp(s, "restrict", 8)) {
1919 if (strlen(line) == 8) {
1920 if (first) {
1921 return 0;
1922 } else {
1923 return E_PARSE;
1924 }
1925 }
1926 s += 8;
1927 while (isspace((unsigned char) *s)) s++;
1928 }
1929
1930 if (!strncmp(s, "rfunc", 5)) {
1931 /* nonlinear function given */
1932 err = read_fncall_line(s + 5, rset);
1933 } else if (*s == 'R' || *s == 'q') {
1934 /* restrictions given in matrix form */
1935 err = read_matrix_line(s, rset);
1936 } else if (!strncmp(s, "inject", 6)) {
1937 /* restrictions given as an array of strings */
1938 err = parse_strings_line(s + 6, rset, dset);
1939 } else {
1940 /* a regular linear restriction */
1941 err = parse_restriction_row(rset, s, dset);
1942 }
1943
1944 if (err) {
1945 destroy_restriction_set(rset);
1946 }
1947
1948 return err;
1949 }
1950
1951 int
restriction_set_parse_line(gretl_restriction * rset,const char * line,const DATASET * dset)1952 restriction_set_parse_line (gretl_restriction *rset, const char *line,
1953 const DATASET *dset)
1954 {
1955 if (rset->g > rset->gmax) {
1956 gretl_errmsg_sprintf(_("Too many restrictions (maximum is %d)"),
1957 rset->gmax);
1958 destroy_restriction_set(rset);
1959 return E_DATA;
1960 }
1961
1962 return real_restriction_set_parse_line(rset, line, dset, 0);
1963 }
1964
1965 /* set-up for a set of restrictions for a VAR (vecm, actually) */
1966
1967 gretl_restriction *
var_restriction_set_start(const char * line,GRETL_VAR * var)1968 var_restriction_set_start (const char *line, GRETL_VAR *var)
1969 {
1970 gretl_restriction *rset;
1971
1972 rset = restriction_set_new(var, GRETL_OBJ_VAR, OPT_NONE);
1973 if (rset == NULL) {
1974 gretl_errmsg_set(_("Out of memory!"));
1975 return NULL;
1976 }
1977
1978 gretl_error_clear();
1979
1980 if (real_restriction_set_parse_line(rset, line, NULL, 1)) {
1981 if (!gretl_errmsg_is_set()) {
1982 gretl_errmsg_sprintf(_("parse error in '%s'\n"), line);
1983 }
1984 return NULL;
1985 }
1986
1987 return rset;
1988 }
1989
1990 /* set-up for a set of (possibly) cross-equation restrictions, for a
1991 system of simultaneous equations */
1992
1993 gretl_restriction *
cross_restriction_set_start(const char * line,equation_system * sys)1994 cross_restriction_set_start (const char *line, equation_system *sys)
1995 {
1996 gretl_restriction *rset;
1997
1998 rset = restriction_set_new(sys, GRETL_OBJ_SYS, OPT_NONE);
1999 if (rset == NULL) {
2000 gretl_errmsg_set(_("Out of memory!"));
2001 return NULL;
2002 }
2003
2004 if (real_restriction_set_parse_line(rset, line, NULL, 1)) {
2005 gretl_errmsg_sprintf(_("parse error in '%s'\n"), line);
2006 return NULL;
2007 }
2008
2009 return rset;
2010 }
2011
2012 /* set-up for a set of restrictions on a single equation */
2013
2014 gretl_restriction *
eqn_restriction_set_start(const char * line,MODEL * pmod,const DATASET * dset,gretlopt opt)2015 eqn_restriction_set_start (const char *line, MODEL *pmod,
2016 const DATASET *dset,
2017 gretlopt opt)
2018 {
2019 gretl_restriction *rset;
2020
2021 rset = restriction_set_new(pmod, GRETL_OBJ_EQN, opt);
2022 if (rset == NULL) {
2023 gretl_errmsg_set(_("Out of memory!"));
2024 return NULL;
2025 }
2026
2027 if (real_restriction_set_parse_line(rset, line, dset, 1)) {
2028 gretl_errmsg_sprintf(_("parse error in '%s'\n"), line);
2029 return NULL;
2030 }
2031
2032 return rset;
2033 }
2034
2035 gretl_restriction *
restriction_set_start(const char * target,gretlopt opt,int * err)2036 restriction_set_start (const char *target, gretlopt opt, int *err)
2037 {
2038 gretl_restriction *rset = NULL;
2039 GretlObjType type = 0;
2040 void *ptr = NULL;
2041
2042 #if RDEBUG
2043 fprintf(stderr, "restriction_set_start: target='%s'\n", target);
2044 #endif
2045
2046 if (target != NULL) {
2047 /* get pointer to named object */
2048 if (!strcmp(target, "$system")) {
2049 ptr = get_anonymous_equation_system();
2050 if (ptr != NULL) {
2051 type = GRETL_OBJ_SYS;
2052 }
2053 } else {
2054 ptr = get_model_object_and_type(target, &type);
2055 }
2056 if (ptr == NULL) {
2057 gretl_errmsg_sprintf("'%s': unrecognized target", target);
2058 *err = E_DATA;
2059 }
2060 } else {
2061 /* get pointer to last-estimated model */
2062 ptr = get_last_model(&type);
2063 if (ptr == NULL) {
2064 *err = E_DATA;
2065 }
2066 }
2067
2068 if (!*err) {
2069 #if RDEBUG
2070 fprintf(stderr, " restriction: ptr = %p, type = %d\n", ptr, type);
2071 #endif
2072 if (type != GRETL_OBJ_EQN && type != GRETL_OBJ_SYS &&
2073 type != GRETL_OBJ_VAR) {
2074 *err = E_DATA;
2075 } else {
2076 rset = restriction_set_new(ptr, type, opt);
2077 if (rset == NULL) {
2078 *err = E_ALLOC;
2079 }
2080 }
2081 }
2082
2083 return rset;
2084 }
2085
print_restricted_estimates(MODEL * pmod,const DATASET * dset,gretl_matrix * vb,gretl_matrix * S,double s2,int nc,int rk,PRN * prn)2086 static int print_restricted_estimates (MODEL *pmod,
2087 const DATASET *dset,
2088 gretl_matrix *vb,
2089 gretl_matrix *S,
2090 double s2,
2091 int nc, int rk,
2092 PRN *prn)
2093 {
2094 const double *b = vb->val;
2095 char **names;
2096 double v, *se;
2097 int i;
2098
2099 names = strings_array_new_with_length(nc, VNAMELEN);
2100 if (names == NULL) {
2101 return E_ALLOC;
2102 }
2103
2104 se = malloc(nc * sizeof *se);
2105 if (se == NULL) {
2106 strings_array_free(names, nc);
2107 return E_ALLOC;
2108 }
2109
2110 pprintf(prn, "%s:\n\n", _("Restricted estimates"));
2111
2112 for (i=0; i<nc; i++) {
2113 v = gretl_matrix_get(S, i, i);
2114 if (fabs(v) < 1.0e-16) {
2115 se[i] = 0.0;
2116 } else {
2117 se[i] = sqrt(v);
2118 }
2119 gretl_model_get_param_name(pmod, dset, i, names[i]);
2120 }
2121
2122 print_coeffs(b, se, (const char **) names, nc, pmod->dfd + rk,
2123 pmod->ci, prn);
2124
2125 pputc(prn, '\n');
2126 pprintf(prn, " %s = %.*g\n", _("Standard error of the regression"),
2127 get_gretl_digits(), sqrt(s2));
2128
2129 strings_array_free(names, nc);
2130 free(se);
2131
2132 return 0;
2133 }
2134
2135 /* Used when generating restricted estimates or bootstrapping, which
2136 we do for single-equation OLS/WLS only. */
2137
rset_expand_R(gretl_restriction * rset,int k)2138 static int rset_expand_R (gretl_restriction *rset, int k)
2139 {
2140 gretl_matrix *R = NULL;
2141 double rij;
2142 int i, j, jj;
2143
2144 R = gretl_zero_matrix_new(rset->g, k);
2145 if (R == NULL) {
2146 return E_ALLOC;
2147 }
2148
2149 for (i=0; i<rset->g; i++) {
2150 jj = 0;
2151 for (j=0; j<k; j++) {
2152 if (rset->mask[j]) {
2153 rij = gretl_matrix_get(rset->R, i, jj);
2154 gretl_matrix_set(R, i, j, rij);
2155 jj++;
2156 }
2157 }
2158 }
2159
2160 gretl_matrix_replace(&rset->R, R);
2161
2162 return 0;
2163 }
2164
rmod_check_ifc(MODEL * rmod,const gretl_matrix * y)2165 static void rmod_check_ifc (MODEL *rmod, const gretl_matrix *y)
2166 {
2167 double x1 = 0.0, x2 = 0.0;
2168 double reldiff;
2169 int t, s = 0;
2170
2171 for (t=rmod->t1; t<=rmod->t2; t++) {
2172 if (!na(rmod->yhat[t])) {
2173 x1 += rmod->yhat[t];
2174 x2 += y->val[s++];
2175 }
2176 }
2177
2178 reldiff = fabs((x1 - x2) / x2);
2179 rmod->ifc = reldiff < 9.0e-16;
2180 }
2181
2182 /* respond to the --full option in case of a restriction on
2183 a model estimated via OLS: arrange to replace the "last model"
2184 with the restricted OLS estimates (and to print the model
2185 if wanted)
2186 */
2187
save_restricted_model(ExecState * state,MODEL * pmod,const DATASET * dset,const gretl_matrix * y,const gretl_matrix * b,const gretl_matrix * S,const gretl_matrix * u,double s2,gretl_restriction * rset,PRN * prn)2188 static int save_restricted_model (ExecState *state,
2189 MODEL *pmod,
2190 const DATASET *dset,
2191 const gretl_matrix *y,
2192 const gretl_matrix *b,
2193 const gretl_matrix *S,
2194 const gretl_matrix *u,
2195 double s2,
2196 gretl_restriction *rset,
2197 PRN *prn)
2198 {
2199 MODEL *rmod;
2200 int i, err = 0;
2201
2202 if (state == NULL) {
2203 return E_DATA;
2204 }
2205
2206 rmod = gretl_model_new();
2207 if (rmod == NULL) {
2208 return E_ALLOC;
2209 }
2210
2211 rmod->ci = OLS;
2212 rmod->ncoeff = gretl_vector_get_length(b);
2213 rmod->nobs = gretl_vector_get_length(y);
2214 rmod->full_n = dset->n;
2215 rmod->t1 = pmod->t1;
2216 rmod->t2 = pmod->t2;
2217
2218 rmod->dfn = pmod->dfn - rset->g;
2219 rmod->dfd = pmod->dfd + rset->g;
2220 rmod->ybar = pmod->ybar;
2221 rmod->sdy = pmod->sdy;
2222 rmod->tss = pmod->tss;
2223
2224 gretl_model_smpl_init(rmod, dset);
2225
2226 err = gretl_model_allocate_storage(rmod);
2227
2228 if (!err) {
2229 for (i=0; i<rmod->ncoeff; i++) {
2230 rmod->coeff[i] = b->val[i];
2231 }
2232 gretl_model_set_int(rmod, "restricted", 1);
2233 err = gretl_model_write_vcv(rmod, S);
2234 }
2235
2236 if (!err) {
2237 rmod->list = gretl_list_copy(pmod->list);
2238 if (rmod->list == NULL) {
2239 err = E_ALLOC;
2240 }
2241 }
2242
2243 if (!err) {
2244 double yt, ut;
2245 int t, s = 0;
2246
2247 rmod->sigma = sqrt(s2);
2248 rmod->ess = 0.0;
2249
2250 for (t=0; t<dset->n; t++) {
2251 if (!na(pmod->uhat[t])) {
2252 yt = gretl_vector_get(y, s);
2253 ut = gretl_vector_get(u, s);
2254 rmod->uhat[t] = ut;
2255 rmod->yhat[t] = yt - ut;
2256 rmod->ess += ut * ut;
2257 s++;
2258 }
2259 }
2260
2261 rmod_check_ifc(rmod, y); /* FIXME? */
2262
2263 if (rmod->ifc) {
2264 double tmp;
2265
2266 rmod->rsq = 1.0 - (rmod->ess / rmod->tss);
2267 if (rmod->rsq < 0.0) {
2268 rmod->rsq = 0.0;
2269 }
2270 if (rmod->dfn == 0) {
2271 rmod->adjrsq = 0.0;
2272 rmod->fstt = NADBL;
2273 } else {
2274 tmp = rmod->tss * rmod->dfd;
2275 rmod->adjrsq = 1 - (rmod->ess * (rmod->nobs - 1) / tmp);
2276 tmp = rmod->dfd / (double) rmod->dfn;
2277 rmod->fstt = tmp * (rmod->tss - rmod->ess) / rmod->ess;
2278 if (rmod->fstt < 0) {
2279 rmod->fstt = NADBL;
2280 }
2281 }
2282 }
2283
2284 rmod->ncoeff -= rset->g;
2285 ls_criteria(rmod);
2286 rmod->ncoeff += rset->g;
2287 }
2288
2289 if (err) {
2290 gretl_model_free(rmod);
2291 } else {
2292 set_model_id(rmod, OPT_NONE);
2293 gretl_exec_state_set_model(state, rmod);
2294 }
2295
2296 return err;
2297 }
2298
2299 /* generate full restricted estimates: this function is used
2300 only for single-equation models, estimated via OLS */
2301
do_restricted_estimates(ExecState * state,gretl_restriction * rset,const DATASET * dset,PRN * prn)2302 static int do_restricted_estimates (ExecState *state,
2303 gretl_restriction *rset,
2304 const DATASET *dset,
2305 PRN *prn)
2306 {
2307 MODEL *pmod = rset->obj;
2308 gretl_matrix_block *B;
2309 gretl_matrix *X, *y, *b, *S;
2310 gretl_matrix *u = NULL;
2311 int *xlist = NULL;
2312 double s2 = 0.0;
2313 int T = pmod->nobs;
2314 int k = pmod->ncoeff;
2315 int i, s, t;
2316 int yno, err = 0;
2317
2318 B = gretl_matrix_block_new(&X, T, k,
2319 &y, T, 1,
2320 &b, k, 1,
2321 &S, k, k,
2322 NULL);
2323
2324 if (B == NULL) {
2325 return E_ALLOC;
2326 }
2327
2328 if (gretl_matrix_cols(rset->R) != k) {
2329 err = rset_expand_R(rset, k);
2330 if (err) {
2331 goto bailout;
2332 }
2333 }
2334
2335 yno = gretl_model_get_depvar(pmod);
2336 xlist = gretl_model_get_x_list(pmod);
2337 if (xlist == NULL) {
2338 err = E_ALLOC;
2339 goto bailout;
2340 }
2341
2342 if (rset->opt & OPT_F) {
2343 /* full save wanted */
2344 u = gretl_matrix_alloc(T, 1);
2345 if (u == NULL) {
2346 err = E_ALLOC;
2347 goto bailout;
2348 }
2349 }
2350
2351 s = 0;
2352 for (t=pmod->t1; t<=pmod->t2; t++) {
2353 if (na(pmod->uhat[t])) {
2354 continue;
2355 }
2356 gretl_vector_set(y, s, dset->Z[yno][t]);
2357 for (i=0; i<k; i++) {
2358 gretl_matrix_set(X, s, i, dset->Z[xlist[i+1]][t]);
2359 }
2360 s++;
2361 }
2362
2363 #if RDEBUG
2364 gretl_matrix_print(rset->R, "R (for restricted estimates)");
2365 gretl_matrix_print(rset->q, "q (for restricted estimates)");
2366 #endif
2367
2368 err = gretl_matrix_restricted_ols(y, X, rset->R, rset->q,
2369 b, S, u, &s2);
2370
2371 if (!err) {
2372 if (rset->opt & OPT_F) {
2373 err = save_restricted_model(state, pmod, dset, y, b, S,
2374 u, s2, rset, prn);
2375 } else {
2376 print_restricted_estimates(pmod, dset, b, S, s2,
2377 k, rset->g, prn);
2378 }
2379 }
2380
2381 bailout:
2382
2383 gretl_matrix_block_destroy(B);
2384 gretl_matrix_free(u);
2385 free(xlist);
2386
2387 return err;
2388 }
2389
2390 /* print or save restricted estimates, single equation */
2391
print_or_save_result(ExecState * state,gretl_restriction * rset,const DATASET * dset,PRN * prn)2392 static void print_or_save_result (ExecState *state,
2393 gretl_restriction *rset,
2394 const DATASET *dset,
2395 PRN *prn)
2396 {
2397 MODEL *pmod = rset->obj;
2398 int robust, asym = 0;
2399 int quiet = (rset->opt & OPT_Q);
2400 int silent = (rset->opt & OPT_S);
2401
2402 robust = (pmod->opt & OPT_R);
2403
2404 if (ASYMPTOTIC_MODEL(pmod->ci) || rset->rfunc != NULL) {
2405 asym = 1;
2406 } else if (pmod->ci == PANEL && (pmod->opt & OPT_U)) {
2407 /* random effects */
2408 asym = 1;
2409 }
2410
2411 if (asym) {
2412 rset->code = GRETL_STAT_WALD_CHISQ;
2413 rset->pval = chisq_cdf_comp(rset->g, rset->test);
2414 if (!silent) {
2415 pprintf(prn, "\n%s: %s(%d) = %g, ", _("Test statistic"),
2416 (robust)? _("Robust chi^2"): "chi^2",
2417 rset->g, rset->test);
2418 }
2419 } else {
2420 rset->code = GRETL_STAT_F;
2421 rset->test /= rset->g;
2422 rset->pval = snedecor_cdf_comp(rset->g, pmod->dfd, rset->test);
2423 if (!silent) {
2424 pprintf(prn, "\n%s: %s(%d, %d) = %g, ", _("Test statistic"),
2425 (robust)? _("Robust F"): "F",
2426 rset->g, pmod->dfd, rset->test);
2427 }
2428 }
2429
2430 if (!na(rset->pval) && !silent) {
2431 pprintf(prn, _("with p-value = %g\n"), rset->pval);
2432 }
2433
2434 if (!silent) {
2435 pputc(prn, '\n');
2436 }
2437
2438 if (!(rset->opt & OPT_C)) {
2439 record_test_result(rset->test, rset->pval);
2440 }
2441
2442 if (pmod != NULL && pmod->ci == OLS &&
2443 dset != NULL && dset->Z != NULL) {
2444 if ((rset->opt & OPT_F) || !(quiet || silent)) {
2445 /* print and/or save restricted estimates */
2446 do_restricted_estimates(state, rset, dset, prn);
2447 }
2448 }
2449 }
2450
2451 /* execute the test, for a single equation */
2452
test_restriction_set(gretl_restriction * rset,PRN * prn)2453 static int test_restriction_set (gretl_restriction *rset, PRN *prn)
2454 {
2455 MODEL *pmod = rset->obj;
2456 gretl_matrix *vcv = NULL;
2457 gretl_vector *b = NULL;
2458 gretl_vector *br = NULL;
2459 gretl_matrix *RvR = NULL;
2460 int err, freeRvR = 1;
2461
2462 gretl_error_clear();
2463
2464 err = restriction_set_form_matrices(rset);
2465 if (err) {
2466 return err;
2467 }
2468
2469 #if RDEBUG
2470 gretl_matrix_print(rset->R, "R matrix");
2471 gretl_matrix_print(rset->q, "q vector");
2472 #endif
2473
2474 err = check_R_matrix(rset->R);
2475
2476 if (!err) {
2477 b = gretl_coeff_vector_from_model(pmod, rset->mask, &err);
2478 }
2479
2480 if (!err) {
2481 vcv = gretl_vcv_matrix_from_model(pmod, rset->mask, &err);
2482 }
2483
2484 if (err) {
2485 goto bailout;
2486 }
2487
2488 if (rset->g == 0) {
2489 rset->g = gretl_matrix_rows(rset->q);
2490 }
2491
2492 br = gretl_column_vector_alloc(rset->g);
2493 if (br == NULL) {
2494 err = E_ALLOC;
2495 goto bailout;
2496 }
2497
2498 #if RDEBUG
2499 gretl_matrix_print(vcv, "VCV matrix");
2500 gretl_matrix_print(b, "coeff vector");
2501 #endif
2502
2503 err = gretl_matrix_multiply(rset->R, b, br);
2504 if (err) {
2505 fprintf(stderr, "Failed: gretl_matrix_multiply(R, b, br)\n");
2506 goto bailout;
2507 }
2508
2509 #if RDEBUG
2510 gretl_matrix_print(br, "br");
2511 #endif
2512
2513 if (rset->opt & OPT_C) {
2514 rset->bsum = br->val[0];
2515 }
2516
2517 if (!gretl_is_zero_matrix(rset->q)) {
2518 err = gretl_matrix_subtract_from(br, rset->q);
2519 if (err) {
2520 fprintf(stderr, "Failed: gretl_matrix_subtract_from(br, q)\n");
2521 goto bailout;
2522 }
2523 }
2524
2525 if (gretl_is_identity_matrix(rset->R)) {
2526 #if RDEBUG
2527 fprintf(stderr, "R is identity matrix: taking shortcut\n");
2528 #endif
2529 RvR = vcv;
2530 freeRvR = 0;
2531 } else {
2532 RvR = gretl_matrix_alloc(rset->R->rows, rset->R->rows);
2533 if (RvR == NULL) {
2534 err = E_ALLOC;
2535 goto bailout;
2536 }
2537 gretl_matrix_qform(rset->R, GRETL_MOD_NONE, vcv,
2538 RvR, GRETL_MOD_NONE);
2539 #if RDEBUG
2540 gretl_matrix_print(RvR, "RvR");
2541 #endif
2542 if (rset->opt & OPT_C) {
2543 rset->bse = sqrt(RvR->val[0]);
2544 }
2545 }
2546
2547 err = gretl_invert_symmetric_matrix(RvR);
2548 if (err) {
2549 pputs(prn, _("Matrix inversion failed:\n"
2550 " restrictions may be inconsistent or redundant\n"));
2551 goto bailout;
2552 }
2553
2554 rset->test = gretl_scalar_qform(br, RvR, &err);
2555 if (err) {
2556 pputs(prn, _("Failed to compute test statistic\n"));
2557 goto bailout;
2558 }
2559
2560 bailout:
2561
2562 gretl_matrix_free(vcv);
2563 gretl_vector_free(b);
2564 gretl_vector_free(br);
2565
2566 if (freeRvR) {
2567 gretl_matrix_free(RvR);
2568 }
2569
2570 return err;
2571 }
2572
bootstrap_zero_restriction(gretl_restriction * rset,MODEL * pmod,const DATASET * dset,int method,PRN * prn)2573 static int bootstrap_zero_restriction (gretl_restriction *rset,
2574 MODEL *pmod,
2575 const DATASET *dset,
2576 int method,
2577 PRN *prn)
2578 {
2579 rrow *r = rset->rows[0];
2580 gretlopt bopt = OPT_P | OPT_R;
2581 int B = 0;
2582 int err = 0;
2583
2584 if (rset->opt & OPT_S) {
2585 /* silent */
2586 bopt |= OPT_S;
2587 }
2588
2589 if (method == BOOT_METHOD_PAIRS) {
2590 bopt |= OPT_X;
2591 } else if (method == BOOT_METHOD_WILD) {
2592 bopt |= OPT_W;
2593 } else if (method == BOOT_METHOD_PARAMETRIC) {
2594 bopt |= OPT_N;
2595 }
2596
2597 gretl_restriction_get_boot_params(&B, &bopt);
2598 err = bootstrap_analysis(pmod, r->bnum[0], B, 0, dset,
2599 bopt, prn);
2600 return err;
2601 }
2602
2603 /* Are we able to do a bootstrap test of the restrictions on
2604 the given model? Return non-zero if not. */
2605
check_bootstrap(MODEL * pmod,PRN * prn)2606 static int check_bootstrap (MODEL *pmod, PRN *prn)
2607 {
2608 int ok = bootstrap_ok(pmod->ci);
2609
2610 if (!ok) {
2611 pputs(prn, "Sorry, the bootstrap option is not supported for this test");
2612 pputc(prn, '\n');
2613 }
2614
2615 return (ok)? 0 : E_NOTIMP;
2616 }
2617
2618 /* Do we have a simple b[i] = 0? */
2619
is_simple_zero_restriction(gretl_restriction * rset)2620 static int is_simple_zero_restriction (gretl_restriction *rset)
2621 {
2622 if (rset->rows != NULL && rset->g == 1) {
2623 return (rset->rows[0]->nterms == 1 && rset->rows[0]->rhs == 0);
2624 } else {
2625 return 0;
2626 }
2627 }
2628
get_restriction_boot_method(int * method)2629 static int get_restriction_boot_method (int *method)
2630 {
2631 const char *s = get_optval_string(RESTRICT, OPT_B);
2632 int err = 0;
2633
2634 if (s != NULL) {
2635 if (!strcmp(s, "pairs")) {
2636 *method = BOOT_METHOD_PAIRS;
2637 } else if (!strcmp(s, "wild")) {
2638 *method = BOOT_METHOD_WILD;
2639 } else if (!strcmp(s, "parametric")) {
2640 *method = BOOT_METHOD_PARAMETRIC;
2641 } else if (strcmp(s, "residuals")) {
2642 gretl_errmsg_sprintf(_("field '%s' in command is invalid"), s);
2643 err = E_DATA;
2644 }
2645 }
2646
2647 return err;
2648 }
2649
do_single_equation_test(ExecState * state,gretl_restriction * rset,const DATASET * dset,PRN * prn)2650 static int do_single_equation_test (ExecState *state,
2651 gretl_restriction *rset,
2652 const DATASET *dset,
2653 PRN *prn)
2654 {
2655 int err = 0;
2656
2657 if (rset->rows != NULL) {
2658 err = restriction_set_make_mask(rset);
2659 if (err) {
2660 return err;
2661 }
2662 }
2663
2664 if (rset->opt & OPT_B) {
2665 /* bootstrapping, if possible */
2666 MODEL *pmod = rset->obj;
2667 int method = 1;
2668
2669 err = check_bootstrap(pmod, prn);
2670
2671 if (!err) {
2672 err = get_restriction_boot_method(&method);
2673 }
2674
2675 if (err) {
2676 return err;
2677 }
2678
2679 if (is_simple_zero_restriction(rset)) {
2680 err = bootstrap_zero_restriction(rset, pmod, dset,
2681 method, prn);
2682 } else {
2683 err = test_restriction_set(rset, prn);
2684 if (!err && rset->R->cols != pmod->ncoeff) {
2685 err = rset_expand_R(rset, pmod->ncoeff);
2686 }
2687 if (!err) {
2688 rset->test /= rset->g;
2689 err = bootstrap_test_restriction(pmod, rset->R, rset->q,
2690 rset->test, rset->g,
2691 dset, rset->opt,
2692 method, prn);
2693 }
2694 }
2695 } else {
2696 err = test_restriction_set(rset, prn);
2697 if (!err) {
2698 print_or_save_result(state, rset, dset, prn);
2699 }
2700 }
2701
2702 return err;
2703 }
2704
gretl_restricted_vecm(gretl_restriction * rset,const DATASET * dset,gretlopt opt,PRN * prn,int * err)2705 GRETL_VAR *gretl_restricted_vecm (gretl_restriction *rset,
2706 const DATASET *dset,
2707 gretlopt opt,
2708 PRN *prn,
2709 int *err)
2710 {
2711 GRETL_VAR *jvar = NULL;
2712
2713 #if RDEBUG
2714 fprintf(stderr, "gretl_restricted_vecm: starting\n");
2715 #endif
2716
2717 if (rset == NULL || rset->otype != GRETL_OBJ_VAR) {
2718 *err = E_DATA;
2719 return NULL;
2720 }
2721
2722 rset->opt |= opt;
2723
2724 *err = restriction_set_form_matrices(rset);
2725
2726 if (rset->rows != NULL && !(opt & OPT_S)) {
2727 print_restriction_set(rset, dset, opt, prn);
2728 }
2729
2730 if (!*err) {
2731 jvar = real_gretl_restricted_vecm(rset->obj, rset, dset,
2732 prn, err);
2733 }
2734
2735 destroy_restriction_set(rset);
2736
2737 if (jvar != NULL && *err != 0) {
2738 gretl_VAR_free(jvar);
2739 jvar = NULL;
2740 }
2741
2742 return jvar;
2743 }
2744
2745 #define RCOEFFNAME "restr__b"
2746
nonlinear_wald_test(gretl_restriction * rset,gretlopt opt,PRN * prn)2747 static int nonlinear_wald_test (gretl_restriction *rset, gretlopt opt,
2748 PRN *prn)
2749 {
2750 gretl_matrix *coeff = NULL;
2751 gretl_matrix *vcv = NULL;
2752 gretl_matrix *J = NULL;
2753 gretl_matrix *bread = NULL;
2754 gretl_matrix *ham = NULL;
2755 char fncall[64];
2756 int published = 0;
2757 int err = 0;
2758
2759 if (get_user_function_by_name(rset->rfunc) == NULL) {
2760 gretl_errmsg_sprintf(_("The symbol '%s' is undefined\n"), rset->rfunc);
2761 return E_UNKVAR;
2762 }
2763
2764 if (rset->otype == GRETL_OBJ_EQN) {
2765 MODEL *pmod = rset->obj;
2766
2767 coeff = gretl_coeff_vector_from_model(pmod, NULL, &err);
2768 if (!err) {
2769 vcv = gretl_vcv_matrix_from_model(pmod, NULL, &err);
2770 }
2771 } else if (rset->otype == GRETL_OBJ_SYS) {
2772 equation_system *sys = rset->obj;
2773
2774 coeff = equation_system_get_matrix(sys, M_COEFF, &err);
2775 if (!err) {
2776 vcv = equation_system_get_matrix(sys, M_VCV, &err);
2777 }
2778 } else {
2779 return E_NOTIMP; /* relax this later? */
2780 }
2781
2782 if (!err) {
2783 /* "publish" the coeff matrix temporarily */
2784 err = private_matrix_add(coeff, RCOEFFNAME);
2785 if (!err) {
2786 published = 1;
2787 }
2788 }
2789
2790 if (!err) {
2791 /* formulate function call string and make the call:
2792 should get a vector */
2793 sprintf(fncall, "%s(%s)", rset->rfunc, RCOEFFNAME);
2794 bread = generate_matrix(fncall, NULL, &err);
2795 }
2796
2797 #if RDEBUG
2798 fprintf(stderr, "nonlinear_wald_test: fncall = '%s'\n", fncall);
2799 gretl_matrix_print(bread, "'bread'");
2800 #endif
2801
2802 if (!err) {
2803 rset->g = gretl_vector_get_length(bread);
2804 if (rset->g == 0) {
2805 err = E_NONCONF;
2806 }
2807 }
2808
2809 if (!err) {
2810 J = user_fdjac(coeff, fncall, 0.0, NULL, &err);
2811 }
2812
2813 #if RDEBUG
2814 gretl_matrix_print(J, "J");
2815 #endif
2816
2817 if (!err) {
2818 ham = gretl_matrix_alloc(J->rows, J->rows);
2819 if (ham == NULL) {
2820 err = E_ALLOC;
2821 } else {
2822 err = gretl_matrix_qform(J, GRETL_MOD_NONE, vcv,
2823 ham, GRETL_MOD_NONE);
2824 }
2825 }
2826
2827 if (!err) {
2828 err = gretl_invpd(ham);
2829 }
2830
2831 if (!err) {
2832 rset->test = gretl_scalar_qform(bread, ham, &err);
2833 }
2834
2835 if (!err) {
2836 print_or_save_result(NULL, rset, NULL, prn);
2837 }
2838
2839 gretl_matrix_free(vcv);
2840 gretl_matrix_free(J);
2841 gretl_matrix_free(bread);
2842 gretl_matrix_free(ham);
2843
2844 if (published) {
2845 destroy_private_matrices();
2846 } else {
2847 gretl_matrix_free(coeff);
2848 }
2849
2850 return err;
2851 }
2852
2853 /* Respond to "end restrict": in the case of a single equation, go
2854 ahead and do the test; for a VECM, hand off to the driver in
2855 var.c; for non-VAR systems, form the restriction matrices R
2856 and q and attach these to the equation system.
2857 */
2858
2859 int
gretl_restriction_finalize_full(ExecState * state,gretl_restriction * rset,const DATASET * dset,gretlopt opt,PRN * prn)2860 gretl_restriction_finalize_full (ExecState *state,
2861 gretl_restriction *rset,
2862 const DATASET *dset,
2863 gretlopt opt,
2864 PRN *prn)
2865 {
2866 int silent, err = 0;
2867
2868 #if RDEBUG
2869 fprintf(stderr, "gretl_restriction_finalize_full: starting\n");
2870 #endif
2871
2872 if (rset == NULL) {
2873 return 1;
2874 }
2875
2876 rset->opt |= opt;
2877 silent = (rset->opt & OPT_S);
2878
2879 if (rset->otype == GRETL_OBJ_EQN && (opt & OPT_F)) {
2880 /* --full option is for OLS only */
2881 MODEL *pmod = rset->obj;
2882
2883 if (pmod->ci != OLS) {
2884 return E_BADOPT;
2885 }
2886 }
2887
2888 if (rset->rfunc != NULL) {
2889 err = nonlinear_wald_test(rset, opt, prn);
2890 destroy_restriction_set(rset);
2891 return err;
2892 }
2893
2894 if (rset->otype != GRETL_OBJ_EQN) {
2895 err = restriction_set_form_matrices(rset);
2896 if (err) {
2897 destroy_restriction_set(rset);
2898 return err;
2899 }
2900 }
2901
2902 if (!silent && rset->rows != NULL) {
2903 print_restriction_set(rset, dset, opt, prn);
2904 }
2905
2906 if (rset->otype == GRETL_OBJ_VAR) {
2907 if (gretl_VECM_rank(rset->obj) > 0) {
2908 err = gretl_VECM_test(rset->obj, rset, dset, rset->opt, prn);
2909 } else {
2910 err = gretl_VAR_test(rset->obj, rset, dset, rset->opt, prn);
2911 }
2912 } else if (rset->otype == GRETL_OBJ_SYS) {
2913 equation_system *sys = rset->obj;
2914
2915 if (sys->flags & SYSTEM_ROBUST) {
2916 /* for now, we'll only do a Wald test */
2917 err = system_wald_test(sys, rset->R, rset->q, opt, prn);
2918 } else if (rset->opt & OPT_W) {
2919 err = system_wald_test(sys, rset->R, rset->q, opt, prn);
2920 } else {
2921 system_set_restriction_matrices(sys, rset->R, rset->q);
2922 rset->R = NULL;
2923 rset->q = NULL;
2924 }
2925 } else {
2926 /* single-equation model */
2927 err = do_single_equation_test(state, rset, dset, prn);
2928 }
2929
2930 if (!(rset->opt & OPT_C)) {
2931 destroy_restriction_set(rset);
2932 }
2933
2934 return err;
2935 }
2936
gretl_restriction_finalize(gretl_restriction * rset,const DATASET * dset,gretlopt opt,PRN * prn)2937 int gretl_restriction_finalize (gretl_restriction *rset,
2938 const DATASET *dset,
2939 gretlopt opt,
2940 PRN *prn)
2941 {
2942 return gretl_restriction_finalize_full(NULL, rset, dset,
2943 opt, prn);
2944 }
2945
2946 /**
2947 * gretl_sum_test:
2948 * @list: list of variables to use.
2949 * @pmod: pointer to model.
2950 * @dset: information on the data set.
2951 * @opt: option flags (OPT_Q gives quiet operation).
2952 * @prn: gretl printing struct.
2953 *
2954 * Calculates the sum of the coefficients, relative to the given model,
2955 * for the variables given in @list. Prints this estimate along
2956 * with its standard error, unless OPT_Q is given.
2957 *
2958 * Returns: 0 on successful completion, error code on error.
2959 */
2960
2961 int
gretl_sum_test(const int * list,MODEL * pmod,DATASET * dset,gretlopt opt,PRN * prn)2962 gretl_sum_test (const int *list, MODEL *pmod, DATASET *dset,
2963 gretlopt opt, PRN *prn)
2964 {
2965 gretl_restriction *r;
2966 char line[MAXLEN];
2967 char bstr[36];
2968 int quiet = (opt & OPT_Q);
2969 int i, len, err = 0;
2970
2971 if (list[0] < 2) {
2972 pprintf(prn, _("Invalid input\n"));
2973 return E_DATA;
2974 }
2975
2976 if (!command_ok_for_model(COEFFSUM, 0, pmod)) {
2977 return E_NOTIMP;
2978 }
2979
2980 if (exact_fit_check(pmod, prn)) {
2981 return 0;
2982 }
2983
2984 r = restriction_set_new(pmod, GRETL_OBJ_EQN, OPT_Q | OPT_C);
2985 if (r == NULL) {
2986 return 1;
2987 }
2988
2989 *line = '\0';
2990 len = 0;
2991
2992 for (i=1; i<=list[0]; i++) {
2993 sprintf(bstr, "b[%s]", dset->varname[list[i]]);
2994 len += strlen(bstr) + 4;
2995 if (len >= MAXLEN - 1) {
2996 err = E_PARSE;
2997 break;
2998 }
2999 strcat(line, bstr);
3000 if (i < list[0]) {
3001 strcat(line, " + ");
3002 } else {
3003 strcat(line, " = 0");
3004 }
3005 }
3006
3007 if (!err) {
3008 err = real_restriction_set_parse_line(r, line, dset, 1);
3009 }
3010
3011 if (!err) {
3012 err = gretl_restriction_finalize(r, dset, OPT_NONE, NULL);
3013 }
3014
3015 if (!err) {
3016 double test;
3017
3018 if (!quiet) {
3019 pprintf(prn, "\n%s: ", _("Variables"));
3020 for (i=1; i<=list[0]; i++) {
3021 pprintf(prn, "%s ", dset->varname[list[i]]);
3022 }
3023 pprintf(prn, "\n %s = %g\n", _("Sum of coefficients"), r->bsum);
3024 }
3025
3026 test = sqrt(r->test);
3027 if (r->bsum < 0) {
3028 test = -test;
3029 }
3030
3031 if (r->code == GRETL_STAT_F) {
3032 if (!quiet) {
3033 pprintf(prn, " %s = %g\n", _("Standard error"), r->bse);
3034 pprintf(prn, " t(%d) = %g ", pmod->dfd, test);
3035 pprintf(prn, _("with p-value = %g\n"), r->pval);
3036 }
3037 record_test_result(test, r->pval);
3038 } else if (r->code == GRETL_STAT_WALD_CHISQ) {
3039 if (!quiet) {
3040 pprintf(prn, " %s = %g\n", _("Standard error"), r->bse);
3041 r->pval = normal_pvalue_2(test);
3042 pprintf(prn, " z = %g ", test);
3043 pprintf(prn, _("with p-value = %g\n"), r->pval);
3044 }
3045 record_test_result(test, r->pval);
3046 }
3047
3048 destroy_restriction_set(r);
3049 }
3050
3051 return err;
3052 }
3053
3054 static int restrict_B;
3055 static gretlopt rboot_opt;
3056
gretl_restriction_set_boot_params(int B,gretlopt opt)3057 int gretl_restriction_set_boot_params (int B, gretlopt opt)
3058 {
3059 int err = 0;
3060
3061 rboot_opt = opt;
3062
3063 if (B > 0) {
3064 restrict_B = B;
3065 } else {
3066 err = E_DATA;
3067 }
3068
3069 return err;
3070 }
3071
gretl_restriction_get_boot_params(int * pB,gretlopt * popt)3072 void gretl_restriction_get_boot_params (int *pB, gretlopt *popt)
3073 {
3074 *pB = restrict_B;
3075 *popt |= rboot_opt;
3076
3077 /* these are ad hoc values */
3078 restrict_B = 0;
3079 rboot_opt = OPT_NONE;
3080 }
3081
gretl_restriction_get_options(const gretl_restriction * rset)3082 gretlopt gretl_restriction_get_options (const gretl_restriction *rset)
3083 {
3084 return (rset != NULL)? rset->opt : OPT_NONE;
3085 }
3086
gretl_restriction_get_type(const gretl_restriction * rset)3087 GretlObjType gretl_restriction_get_type (const gretl_restriction *rset)
3088 {
3089 return (rset != NULL)? rset->otype : GRETL_OBJ_NULL;
3090 }
3091
rset_get_R_matrix(const gretl_restriction * rset)3092 const gretl_matrix *rset_get_R_matrix (const gretl_restriction *rset)
3093 {
3094 return (rset != NULL)? rset->R : NULL;
3095 }
3096
rset_get_q_matrix(const gretl_restriction * rset)3097 const gretl_matrix *rset_get_q_matrix (const gretl_restriction *rset)
3098 {
3099 return (rset != NULL)? rset->q : NULL;
3100 }
3101
rset_get_Ra_matrix(const gretl_restriction * rset)3102 const gretl_matrix *rset_get_Ra_matrix (const gretl_restriction *rset)
3103 {
3104 return (rset != NULL)? rset->Ra : NULL;
3105 }
3106
rset_get_qa_matrix(const gretl_restriction * rset)3107 const gretl_matrix *rset_get_qa_matrix (const gretl_restriction *rset)
3108 {
3109 return (rset != NULL)? rset->qa : NULL;
3110 }
3111
rset_VECM_bcols(const gretl_restriction * rset)3112 int rset_VECM_bcols (const gretl_restriction *rset)
3113 {
3114 return (rset == NULL || rset->R == NULL)? 0 : rset->R->cols;
3115 }
3116
rset_VECM_acols(const gretl_restriction * rset)3117 int rset_VECM_acols (const gretl_restriction *rset)
3118 {
3119 return (rset == NULL || rset->Ra == NULL)? 0 : rset->Ra->cols;
3120 }
3121
rset_add_results(gretl_restriction * rset,double test,double pval,double lnl)3122 void rset_add_results (gretl_restriction *rset,
3123 double test, double pval,
3124 double lnl)
3125 {
3126 rset->test = test;
3127 rset->pval = pval;
3128 rset->lnl = lnl;
3129 }
3130
rset_record_LR_result(gretl_restriction * rset)3131 void rset_record_LR_result (gretl_restriction *rset)
3132 {
3133 record_LR_test_result(rset->test, rset->pval, rset->lnl);
3134 }
3135