1 /*
2 * gretl -- Gnu Regression, Econometrics and Time-series Library
3 * Copyright (C) 2001 Allin Cottrell and Riccardo "Jack" Lucchetti
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, either version 3 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program. If not, see <http://www.gnu.org/licenses/>.
17 *
18 */
19
20 #define FULL_XML_HEADERS
21
22 #include "libgretl.h"
23 #include "gretl_xml.h"
24 #include "matrix_extra.h"
25 #include "libset.h"
26 #include "gretl_func.h"
27 #include "gretl_cmatrix.h"
28 #include "gretl_panel.h"
29 #include "gretl_array.h"
30 #include "qr_estimate.h"
31
32 /**
33 * SECTION:gretl_model
34 * @short_description: handling of the MODEL struct
35 * @title: Model structure
36 * @include: libgretl.h
37 *
38 * Provides underlying functionality for gretl's MODEL datatype.
39 */
40
41 #define MDEBUG 0
42
43 struct model_data_item_ {
44 char *key;
45 void *ptr;
46 int type;
47 size_t size;
48 void (*destructor) (void *);
49 };
50
51 struct ModelTest_ {
52 int type;
53 int order;
54 char *param;
55 unsigned char teststat;
56 int dfn;
57 double dfd;
58 double value;
59 double pvalue;
60 double crit;
61 double alpha;
62 gretlopt opt;
63 };
64
65 typedef struct CoeffSep_ CoeffSep;
66
67 #define CSLEN 64
68
69 struct CoeffSep_ {
70 char str[CSLEN];
71 int pos;
72 };
73
74 #define PNAMELEN 16 /* for parameter names */
75
76 static gretl_bundle *bundlize_test (const ModelTest *src);
77
gretl_test_init(ModelTest * test,ModelTestType ttype)78 static void gretl_test_init (ModelTest *test, ModelTestType ttype)
79 {
80 test->type = ttype;
81 test->order = 0;
82 test->param = NULL;
83 test->teststat = GRETL_STAT_NONE;
84 test->dfn = 0;
85 test->dfd = 0;
86 test->value = test->pvalue = NADBL;
87 test->crit = test->alpha = NADBL;
88 test->opt = OPT_NONE;
89 }
90
test_type_key(ModelTestType t)91 static const char *test_type_key (ModelTestType t)
92 {
93 if (t == GRETL_TEST_ADD) {
94 return "add_test";
95 } else if (t == GRETL_TEST_ARCH) {
96 return "arch_test";
97 } else if (t == GRETL_TEST_AUTOCORR) {
98 return "autocorr_test";
99 } else if (t == GRETL_TEST_CHOW ||
100 t == GRETL_TEST_CHOWDUM) {
101 return "chow_test";
102 } else if (t == GRETL_TEST_CUSUM) {
103 return "cusum_test";
104 } else if (t == GRETL_TEST_QLR) {
105 return "qlr_test";
106 } else if (t == GRETL_TEST_GROUPWISE) {
107 return "grpwise_test";
108 } else if (t == GRETL_TEST_LOGS) {
109 return "logs_test";
110 } else if (t == GRETL_TEST_NORMAL) {
111 return "normality_test";
112 } else if (t == GRETL_TEST_OMIT) {
113 return "omit_test";
114 } else if (t == GRETL_TEST_RESET) {
115 return "reset_test";
116 } else if (t == GRETL_TEST_SQUARES) {
117 return "squares_test";
118 } else if (t == GRETL_TEST_WHITES) {
119 return "whites_test";
120 } else if (t == GRETL_TEST_SARGAN) {
121 return "sargan_test";
122 } else if (t == GRETL_TEST_IV_HAUSMAN ||
123 t == GRETL_TEST_PANEL_HAUSMAN) {
124 return "hausman_test";
125 } else if (t == GRETL_TEST_PANEL_F ||
126 t == GRETL_TEST_PANEL_WELCH) {
127 return "fixed_effects_F";
128 } else if (t == GRETL_TEST_PANEL_BP ||
129 t == GRETL_TEST_BP) {
130 return "bp_test";
131 } else if (t == GRETL_TEST_PANEL_TIMEDUM) {
132 return "timedum_test";
133 } else if (t == GRETL_TEST_HET_1) {
134 return "het1_test";
135 } else if (t == GRETL_TEST_COMFAC) {
136 return "comfac_test";
137 } else if (t == GRETL_TEST_INDEP) {
138 return "independence_test";
139 } else if (t == GRETL_TEST_RE) {
140 return "rho_test";
141 } else if (t == GRETL_TEST_WITHIN_F) {
142 return "within_F";
143 } else if (t == GRETL_TEST_RE_WALD) {
144 return "re_wald_test";
145 } else if (t == GRETL_TEST_XDEPEND) {
146 return "cross_sectional_dependence_test";
147 } else if (t == GRETL_TEST_PANEL_AR) {
148 return "panel_ar_test";
149 } else {
150 fprintf(stderr, "test_type_key(): type %d has no key!\n", t);
151 return NULL;
152 }
153 }
154
free_model_data_item(model_data_item * item)155 static void free_model_data_item (model_data_item *item)
156 {
157 if (item->destructor != NULL) {
158 (*item->destructor)(item->ptr);
159 } else {
160 free(item->ptr);
161 }
162 free(item->key);
163 free(item);
164 }
165
create_data_item(const char * key,void * ptr,GretlType type,size_t size,void (* destructor)(void *))166 static model_data_item *create_data_item (const char *key, void *ptr,
167 GretlType type, size_t size,
168 void (*destructor) (void *))
169 {
170 model_data_item *item = malloc(sizeof *item);
171
172 if (item != NULL) {
173 item->key = gretl_strdup(key);
174 if (item->key == NULL) {
175 free(item);
176 item = NULL;
177 } else {
178 item->ptr = ptr;
179 item->type = type;
180 item->size = size;
181 item->destructor = destructor;
182 }
183 }
184
185 return item;
186 }
187
188 static model_data_item *
replicate_data_item(const model_data_item * orig)189 replicate_data_item (const model_data_item *orig)
190 {
191 model_data_item *item = malloc(sizeof *item);
192 int err = 0;
193
194 if (item != NULL) {
195 item->key = gretl_strdup(orig->key);
196 if (item->key == NULL) {
197 free(item);
198 item = NULL;
199 }
200 }
201
202 if (item != NULL) {
203 if (orig->type == GRETL_TYPE_MATRIX) {
204 item->ptr = gretl_matrix_copy(orig->ptr);
205 } else if (orig->type == GRETL_TYPE_ARRAY) {
206 item->ptr = gretl_array_copy(orig->ptr, &err);
207 } else {
208 item->ptr = malloc(orig->size);
209 }
210 if (item->ptr == NULL) {
211 free(item->key);
212 free(item);
213 item = NULL;
214 }
215 }
216
217 if (item != NULL) {
218 if (orig->type != GRETL_TYPE_MATRIX &&
219 orig->type != GRETL_TYPE_ARRAY) {
220 memcpy(item->ptr, orig->ptr, orig->size);
221 }
222 item->type = orig->type;
223 item->size = orig->size;
224 item->destructor = orig->destructor;
225 }
226
227 return item;
228 }
229
230 /**
231 * gretl_model_set_data_with_destructor:
232 * @pmod: pointer to #MODEL.
233 * @key: key string for data, used in retrieval.
234 * @ptr: data-pointer to be attached to model.
235 * @type: type of data to set.
236 * @size: size of data in bytes.
237 * @destructor: pointer to function that should be used to free
238 * the data-pointer in question.
239 *
240 * Attaches data to @pmod: the data can be retrieved later using
241 * gretl_model_get_data(). Note that the data are not "physically"
242 * copied to the model; simply, @ptr is recorded on the model.
243 * This means that the data referenced by the pointer now in
244 * effect belong to @pmod. When @pmod is cleared with clear_model(),
245 * @destructor will be invoked with @ptr as its single argument.
246 * If a simple "free" is OK for freeing the data, you can use
247 * gretl_model_set_data() instead.
248 *
249 * The @size is needed in case the model is copied with
250 * copy_model(), in which case the target of the copying
251 * operation receives a newly allocated copy of the data in
252 * question.
253 *
254 * Returns: 0 on success, 1 on failure.
255 */
256
gretl_model_set_data_with_destructor(MODEL * pmod,const char * key,void * ptr,GretlType type,size_t size,void (* destructor)(void *))257 int gretl_model_set_data_with_destructor (MODEL *pmod, const char *key, void *ptr,
258 GretlType type, size_t size,
259 void (*destructor) (void *))
260 {
261 model_data_item **items;
262 model_data_item *item;
263 int i, n;
264
265 for (i=0; i<pmod->n_data_items; i++) {
266 item = pmod->data_items[i];
267 if (!strcmp(key, item->key)) {
268 /* there's a pre-existing item of this name */
269 if (item->destructor != NULL) {
270 (*item->destructor)(item->ptr);
271 } else {
272 free(item->ptr);
273 }
274 item->type = type;
275 item->ptr = ptr;
276 item->size = size;
277 item->destructor = destructor;
278 /* handled */
279 return 0;
280 }
281 }
282
283 n = pmod->n_data_items + 1;
284
285 items = realloc(pmod->data_items, n * sizeof *items);
286 if (items == NULL) {
287 return 1;
288 }
289
290 pmod->data_items = items;
291
292 item = create_data_item(key, ptr, type, size, destructor);
293 if (item == NULL) {
294 return 1;
295 }
296
297 pmod->data_items[n - 1] = item;
298 pmod->n_data_items += 1;
299
300 return 0;
301 }
302
303 /**
304 * gretl_model_set_data:
305 * @pmod: pointer to #MODEL.
306 * @key: key string for data, used in retrieval.
307 * @ptr: data-pointer to be attached to model.
308 * @type: type of the data to set.
309 * @size: size of data in bytes.
310 *
311 * Attaches data to @pmod: the data can be retrieved later using
312 * gretl_model_get_data(). Note that the data are not "physically"
313 * copied to the model; simply, @ptr is recorded on the model.
314 * This means that the data referenced by the pointer now in
315 * effect belong to @pmod. The data pointer will be freed when
316 * @pmod is cleared with clear_model(). If the data has deep
317 * structure that requires special treatment on freeing, use
318 * gretl_model_set_data_with_destructor() instead.
319 *
320 * The @size is needed in case the model is copied with
321 * copy_model(), in which case the target of the copying
322 * operation receives a newly allocated copy of the data in
323 * question.
324 *
325 * Returns: 0 on success, 1 on failure.
326 */
327
gretl_model_set_data(MODEL * pmod,const char * key,void * ptr,GretlType type,size_t size)328 int gretl_model_set_data (MODEL *pmod, const char *key, void *ptr,
329 GretlType type, size_t size)
330 {
331 return gretl_model_set_data_with_destructor(pmod, key, ptr, type,
332 size, NULL);
333 }
334
matrix_free_callback(void * p)335 static void matrix_free_callback (void *p)
336 {
337 gretl_matrix_free((gretl_matrix *) p);
338 }
339
340 /**
341 * gretl_model_set_matrix_as_data:
342 * @pmod: pointer to #MODEL.
343 * @key: key string, used in retrieval.
344 * @m: matrix to attach.
345 *
346 * Attaches @m to @pmod as data, recoverable via the key @key
347 * using gretl_model_get_data().
348 *
349 * Returns: 0 on success, 1 on failure.
350 */
351
gretl_model_set_matrix_as_data(MODEL * pmod,const char * key,gretl_matrix * m)352 int gretl_model_set_matrix_as_data (MODEL *pmod, const char *key,
353 gretl_matrix *m)
354 {
355 return gretl_model_set_data_with_destructor(pmod, key, (void *) m,
356 GRETL_TYPE_MATRIX, 0,
357 matrix_free_callback);
358 }
359
array_free_callback(void * p)360 static void array_free_callback (void *p)
361 {
362 gretl_array_destroy((gretl_array *) p);
363 }
364
365 /**
366 * gretl_model_set_array_as_data:
367 * @pmod: pointer to #MODEL.
368 * @key: key string, used in retrieval.
369 * @A: array to attach.
370 *
371 * Attaches @A to @pmod as data, recoverable via the key @key
372 * using gretl_model_get_data().
373 *
374 * Returns: 0 on success, 1 on failure.
375 */
376
gretl_model_set_array_as_data(MODEL * pmod,const char * key,gretl_array * A)377 int gretl_model_set_array_as_data (MODEL *pmod, const char *key,
378 gretl_array *A)
379 {
380 return gretl_model_set_data_with_destructor(pmod, key, (void *) A,
381 GRETL_TYPE_ARRAY, 0,
382 array_free_callback);
383 }
384
385 /**
386 * gretl_model_set_list_as_data:
387 * @pmod: pointer to #MODEL.
388 * @key: key string, used in retrieval.
389 * @list: list to attach.
390 *
391 * Attaches @list to @pmod as data, recoverable via the key @key
392 * using gretl_model_get_list(). Note that the model takes
393 * ownership of the supplied list.
394 *
395 * Returns: 0 on success, 1 on failure.
396 */
397
gretl_model_set_list_as_data(MODEL * pmod,const char * key,int * list)398 int gretl_model_set_list_as_data (MODEL *pmod, const char *key, int *list)
399 {
400 size_t size = (list[0] + 1) * sizeof *list;
401
402 return gretl_model_set_data_with_destructor(pmod, key, (void *) list,
403 GRETL_TYPE_LIST, size,
404 NULL);
405 }
406
407 /**
408 * gretl_model_set_string_as_data:
409 * @pmod: pointer to #MODEL.
410 * @key: key string, used in retrieval.
411 * @str: string to attach.
412 *
413 * Attaches @str to @pmod as data, recoverable via the key @key
414 * using gretl_model_get_data().
415 *
416 * Returns: 0 on success, 1 on failure.
417 */
418
gretl_model_set_string_as_data(MODEL * pmod,const char * key,char * str)419 int gretl_model_set_string_as_data (MODEL *pmod, const char *key, char *str)
420 {
421 size_t size = strlen(str) + 1;
422
423 return gretl_model_set_data_with_destructor(pmod, key, (void *) str,
424 GRETL_TYPE_STRING, size,
425 NULL);
426 }
427
428 /**
429 * gretl_model_set_int:
430 * @pmod: pointer to #MODEL.
431 * @key: key string, used in retrieval.
432 * @val: integer value to set.
433 *
434 * Records an integer value on a model: the value can be retrieved
435 * later using gretl_model_get_int(), using the appropriate @key.
436 *
437 * Returns: 0 on success, 1 on failure.
438 */
439
gretl_model_set_int(MODEL * pmod,const char * key,int val)440 int gretl_model_set_int (MODEL *pmod, const char *key, int val)
441 {
442 int *valp;
443 int err;
444
445 /* if value is already set, reset it */
446 valp = gretl_model_get_data(pmod, key);
447 if (valp != NULL) {
448 *valp = val;
449 return 0;
450 }
451
452 valp = malloc(sizeof *valp);
453 if (valp == NULL) return 1;
454
455 *valp = val;
456
457 err = gretl_model_set_data(pmod, key, valp, GRETL_TYPE_INT,
458 sizeof(int));
459 if (err) {
460 free(valp);
461 }
462
463 return err;
464 }
465
466 /**
467 * gretl_model_set_double:
468 * @pmod: pointer to model.
469 * @key: key string, used in retrieval.
470 * @val: double-precision value to set.
471 *
472 * Records a floating-point value on @pmod: the value can be
473 * retrieved later using gretl_model_get_double() with the
474 * appropriate @key.
475 *
476 * Returns: 0 on success, 1 on failure.
477 */
478
gretl_model_set_double(MODEL * pmod,const char * key,double val)479 int gretl_model_set_double (MODEL *pmod, const char *key, double val)
480 {
481 double *valp;
482 int err;
483
484 /* if value is already set, reset it */
485 valp = gretl_model_get_data(pmod, key);
486 if (valp != NULL) {
487 *valp = val;
488 return 0;
489 }
490
491 valp = malloc(sizeof *valp);
492 if (valp == NULL) return 1;
493
494 *valp = val;
495
496 err = gretl_model_set_data(pmod, key, valp, GRETL_TYPE_DOUBLE,
497 sizeof(double));
498 if (err) free(valp);
499
500 return err;
501 }
502
vcv_info_new(void)503 static VCVInfo *vcv_info_new (void)
504 {
505 VCVInfo *vi;
506
507 vi = malloc(sizeof *vi);
508
509 if (vi != NULL) {
510 vi->vmaj = vi->vmin = 0;
511 vi->order = vi->flags = 0;
512 vi->bw = NADBL;
513 }
514
515 return vi;
516 }
517
518 /**
519 * gretl_model_set_full_vcv_info:
520 *
521 * Returns: 0 on success, 1 on failure.
522 */
523
gretl_model_set_full_vcv_info(MODEL * pmod,int vmaj,int vmin,int order,int flags,double bw)524 int gretl_model_set_full_vcv_info (MODEL *pmod, int vmaj, int vmin,
525 int order, int flags, double bw)
526 {
527 VCVInfo *vi;
528 int prev = 0;
529 int err = 0;
530
531 vi = gretl_model_get_data(pmod, "vcv_info");
532
533 if (vi == NULL) {
534 vi = vcv_info_new();
535 if (vi == NULL) {
536 return E_ALLOC;
537 }
538 } else {
539 prev = 1;
540 }
541
542 vi->vmaj = vmaj;
543 vi->vmin = vmin;
544 vi->order = order;
545 vi->flags = flags;
546 vi->bw = bw;
547
548 if (!prev) {
549 err = gretl_model_set_data(pmod, "vcv_info", vi,
550 GRETL_TYPE_STRUCT,
551 sizeof *vi);
552 }
553
554 return err;
555 }
556
557 /**
558 * gretl_model_set_vcv_info:
559 * @pmod: pointer to model.
560 *
561 * Returns: 0 on success, 1 on failure.
562 */
563
gretl_model_set_vcv_info(MODEL * pmod,int vmaj,int vmin)564 int gretl_model_set_vcv_info (MODEL *pmod, int vmaj, int vmin)
565 {
566 return gretl_model_set_full_vcv_info(pmod, vmaj, vmin,
567 0, 0, 0);
568 }
569
570 /**
571 * gretl_model_get_vcv_type:
572 * @pmod: pointer to model.
573 *
574 * Returns: index of variance-covariance matrix type.
575 */
576
gretl_model_get_vcv_type(const MODEL * pmod)577 int gretl_model_get_vcv_type (const MODEL *pmod)
578 {
579 VCVInfo *vi = gretl_model_get_data(pmod, "vcv_info");
580
581 if (vi != NULL) {
582 return vi->vmaj;
583 } else {
584 return 0;
585 }
586 }
587
588 /**
589 * gretl_model_get_hc_version:
590 * @pmod: pointer to model.
591 *
592 * Returns: the "HC" (HCCME) variant employed in @pmod,
593 * or -1 if this does not apply.
594 */
595
gretl_model_get_hc_version(const MODEL * pmod)596 int gretl_model_get_hc_version (const MODEL *pmod)
597 {
598 VCVInfo *vi = gretl_model_get_data(pmod, "vcv_info");
599
600 if (vi != NULL && vi->vmaj == VCV_HC) {
601 return vi->vmin;
602 } else {
603 return -1;
604 }
605 }
606
607 /**
608 * gretl_model_get_cluster_var:
609 * @pmod: pointer to model.
610 *
611 * Returns: the dataset index of the clustering variable used for
612 * the variance-covariance matrix in @pmod, or 0 if there
613 * is no such variable.
614 */
615
gretl_model_get_cluster_var(const MODEL * pmod)616 int gretl_model_get_cluster_var (const MODEL *pmod)
617 {
618 VCVInfo *vi = gretl_model_get_data(pmod, "vcv_info");
619
620 if (vi != NULL && vi->vmaj == VCV_CLUSTER) {
621 return vi->vmin;
622 } else {
623 return 0;
624 }
625 }
626
627 /**
628 * gretl_model_get_data_full:
629 * @pmod: pointer to model.
630 * @key: key string.
631 * @copied: location to receive flag indicating whether the
632 * return value is an allocated copy of the original data.
633 * @type: location to receive data type.
634 * @sz: location to receive the size of the data.
635 *
636 * Returns: the data pointer identified by @key, or %NULL on failure.
637 * If a non-zero value is written to @copied this indicates that the
638 * return value is a copy of the original (and therefore it is the
639 * caller's responsibility to free the data when it is no longer
640 * required).
641 */
642
gretl_model_get_data_full(const MODEL * pmod,const char * key,GretlType * type,int * copied,size_t * sz)643 void *gretl_model_get_data_full (const MODEL *pmod, const char *key,
644 GretlType *type, int *copied,
645 size_t *sz)
646 {
647 void *ret = NULL;
648 GretlType itype = 0;
649 size_t isize = 0;
650 int alloced = 0;
651 int i, found = 0;
652
653 if (pmod == NULL) {
654 return NULL;
655 }
656
657 for (i=0; i<pmod->n_data_items; i++) {
658 if (!strcmp(key, pmod->data_items[i]->key)) {
659 ret = pmod->data_items[i]->ptr;
660 itype = pmod->data_items[i]->type;
661 isize = pmod->data_items[i]->size;
662 found = 1;
663 break;
664 }
665 }
666
667 if (!found && pmod->tests != NULL) {
668 const ModelTest *test;
669 const char *tkey;
670 gretl_bundle *b;
671
672 for (i=0; i<pmod->ntests; i++) {
673 test = &pmod->tests[i];
674 tkey = test_type_key(test->type);
675 if (tkey != NULL && !strcmp(key, tkey)) {
676 b = bundlize_test(test);
677 if (b != NULL) {
678 ret = b;
679 itype = GRETL_TYPE_BUNDLE;
680 alloced = 1;
681 }
682 break;
683 }
684 }
685 }
686
687 if (ret != NULL) {
688 if (type != NULL) {
689 *type = itype;
690 }
691 if (sz != NULL) {
692 *sz = isize;
693 }
694 if (copied != NULL) {
695 *copied = alloced;
696 }
697 }
698
699 return ret;
700 }
701
702 /**
703 * gretl_model_get_data:
704 * @pmod: pointer to model.
705 * @key: key string.
706 *
707 * Returns: the data pointer identified by @key, or %NULL on failure.
708 */
709
gretl_model_get_data(const MODEL * pmod,const char * key)710 void *gretl_model_get_data (const MODEL *pmod, const char *key)
711 {
712 return gretl_model_get_data_full(pmod, key, NULL, NULL, NULL);
713 }
714
715 /**
716 * gretl_model_get_int:
717 * @pmod: pointer to model.
718 * @key: key string.
719 *
720 * Returns: the integer value identified by @key, or 0 on failure.
721 */
722
gretl_model_get_int(const MODEL * pmod,const char * key)723 int gretl_model_get_int (const MODEL *pmod, const char *key)
724 {
725 int *valp = NULL;
726 int i;
727
728 if (pmod == NULL) {
729 return 0;
730 }
731
732 for (i=0; i<pmod->n_data_items; i++) {
733 if (pmod->data_items[i]->type != GRETL_TYPE_INT) {
734 continue;
735 }
736 if (!strcmp(key, pmod->data_items[i]->key)) {
737 valp = (int *) pmod->data_items[i]->ptr;
738 return *valp;
739 }
740 }
741
742 return 0;
743 }
744
745 /**
746 * gretl_model_get_double:
747 * @pmod: pointer to model.
748 * @key: key string.
749 *
750 * Returns: the double-precision value identified by @key, or
751 * #NADBL on failure.
752 */
753
gretl_model_get_double(const MODEL * pmod,const char * key)754 double gretl_model_get_double (const MODEL *pmod, const char *key)
755 {
756 double *valp = NULL;
757 int i;
758
759 if (pmod == NULL) {
760 return NADBL;
761 }
762
763 for (i=0; i<pmod->n_data_items; i++) {
764 if (pmod->data_items[i]->type != GRETL_TYPE_DOUBLE) {
765 continue;
766 }
767 if (!strcmp(key, pmod->data_items[i]->key)) {
768 valp = (double *) pmod->data_items[i]->ptr;
769 return *valp;
770 }
771 }
772
773 return NADBL;
774 }
775
776 /**
777 * gretl_model_get_double_default:
778 * @pmod: pointer to model.
779 * @key: key string.
780 * @deflt: default value
781 *
782 * Returns: the double-precision value identified by @key, or
783 * @deflt if there is no such value.
784 */
785
gretl_model_get_double_default(const MODEL * pmod,const char * key,double deflt)786 double gretl_model_get_double_default (const MODEL *pmod,
787 const char *key,
788 double deflt)
789 {
790 double *valp = NULL;
791 double ret = deflt;
792 int i;
793
794 if (pmod == NULL) {
795 return NADBL;
796 }
797
798 for (i=0; i<pmod->n_data_items; i++) {
799 if (pmod->data_items[i]->type != GRETL_TYPE_DOUBLE) {
800 continue;
801 }
802 if (!strcmp(key, pmod->data_items[i]->key)) {
803 valp = (double *) pmod->data_items[i]->ptr;
804 ret = *valp;
805 break;
806 }
807 }
808
809 return ret;
810 }
811
812 /**
813 * gretl_model_get_list:
814 * @pmod: pointer to model.
815 * @key: key string.
816 *
817 * Returns: the list of integers identified by @key, or
818 * %NULL on failure.
819 */
820
gretl_model_get_list(const MODEL * pmod,const char * key)821 int *gretl_model_get_list (const MODEL *pmod, const char *key)
822 {
823 int *list = NULL;
824 int i;
825
826 if (pmod == NULL) {
827 return NULL;
828 }
829
830 for (i=0; i<pmod->n_data_items; i++) {
831 if (pmod->data_items[i]->type != GRETL_TYPE_LIST) {
832 continue;
833 }
834 if (!strcmp(key, pmod->data_items[i]->key)) {
835 list = (int *) pmod->data_items[i]->ptr;
836 break;
837 }
838 }
839
840 return list;
841 }
842
843 /**
844 * gretl_model_set_coeff_separator:
845 * @pmod: pointer to model.
846 * @s: informative string (or %NULL).
847 * @pos: position in the array of coefficients.
848 *
849 * Arranges for the insertion of the given string (or a blank line if
850 * @s is %NULL) at the given position in the array of coefficients,
851 * when the model is printed. The extra line is printed before
852 * coefficient @pos, where @pos is interpreted as a zero-based
853 * index.
854 *
855 * Returns: 0 on success, %E_ALLOC on failure.
856 */
857
gretl_model_set_coeff_separator(MODEL * pmod,const char * s,int pos)858 int gretl_model_set_coeff_separator (MODEL *pmod, const char *s, int pos)
859 {
860 CoeffSep *cs = malloc(sizeof *cs);
861 int err = 0;
862
863 if (cs == NULL) {
864 return E_ALLOC;
865 }
866
867 cs->str[0] = '\0';
868 if (s != NULL) {
869 strncat(cs->str, s, CSLEN - 1);
870 }
871 cs->pos = pos;
872
873 err = gretl_model_set_data(pmod, "coeffsep", cs, GRETL_TYPE_STRUCT,
874 sizeof *cs);
875 if (err) {
876 free(cs);
877 }
878
879 return err;
880 }
881
882 /**
883 * gretl_model_get_coeff_separator:
884 * @pmod: pointer to model.
885 * @ps: location to receive string, if any, or NULL.
886 * @ppos: location to receive position, if any, or NULL.
887 *
888 * Retrieves information that has been set on @pmod regarding the
889 * insertion of an extra line when printing the coefficients, if any.
890 *
891 * Returns: 1 if such information is present, 0 otherwise.
892 */
893
gretl_model_get_coeff_separator(const MODEL * pmod,char ** ps,int * ppos)894 int gretl_model_get_coeff_separator (const MODEL *pmod, char **ps, int *ppos)
895 {
896 CoeffSep *cs = gretl_model_get_data(pmod, "coeffsep");
897
898 if (cs == NULL) {
899 return 0;
900 }
901
902 if (ps != NULL) {
903 *ps = gretl_strdup(cs->str);
904 }
905 if (ppos != NULL) {
906 *ppos = cs->pos;
907 }
908
909 return 1;
910 }
911
make_cname(const char * orig,char * cname)912 static void make_cname (const char *orig, char *cname)
913 {
914 char *p;
915
916 if (orig == NULL || *orig == 0) {
917 return;
918 }
919
920 p = strrchr(orig, '_');
921
922 if (p == NULL) {
923 strcpy(cname, orig);
924 } else {
925 unsigned char c = (unsigned char) *(p + 1);
926
927 if (isdigit(c)) {
928 int lag = atoi(++p);
929
930 sprintf(cname, "ut^2(-%d)", lag);
931 }
932 }
933 }
934
plain_ar_coeff_name(char * targ,const MODEL * pmod,int i)935 static void plain_ar_coeff_name (char *targ, const MODEL *pmod, int i)
936 {
937 int k = i - pmod->ncoeff;
938
939 if (pmod->arinfo != NULL && pmod->arinfo->arlist != NULL &&
940 k >= 0 && k < pmod->arinfo->arlist[0]) {
941 sprintf(targ, "u_%d", pmod->arinfo->arlist[k+1]);
942 } else {
943 strcpy(targ, "unknown");
944 }
945 }
946
947 /**
948 * gretl_model_get_param_name:
949 * @pmod: pointer to model.
950 * @dset: dataset information.
951 * @i: index number for parameter, zero-based, corresponding
952 * to position in the %coeff array in @pmod.
953 * @targ: string into which to write param name.
954 *
955 * Writes the appropriate parameter name into @targ, which
956 * should be at least #VNAMELEN bytes long. Usually this is
957 * the name of a variable in the dataset, but sometimes it is
958 * a special string (e.g. for nonlinear models).
959 *
960 * Returns: @targ.
961 */
962
gretl_model_get_param_name(const MODEL * pmod,const DATASET * dset,int i,char * targ)963 char *gretl_model_get_param_name (const MODEL *pmod,
964 const DATASET *dset,
965 int i, char *targ)
966 {
967 *targ = '\0';
968
969 if (pmod != NULL) {
970 int j = i + 2;
971 int k = -1;
972
973 #if 0
974 /* should we do this? */
975 if (pmod->dataset != NULL) {
976 dset = pmod->dataset;
977 }
978 #endif
979
980 if (pmod->aux == AUX_ARCH) {
981 make_cname(dset->varname[pmod->list[j]], targ);
982 } else if (pmod->ci == PANEL && (pmod->opt & OPT_H)) {
983 /* --unit-weights */
984 strcpy(targ, dset->varname[pmod->list[j]]);
985 } else if (NONLIST_MODEL(pmod->ci) ||
986 pmod->ci == ARMA || pmod->ci == PANEL ||
987 pmod->ci == DPANEL || pmod->ci == GARCH ||
988 pmod->ci == BIPROBIT) {
989 k = i;
990 } else if (pmod->ci == MPOLS && pmod->params != NULL) {
991 k = i;
992 } else if ((pmod->ci == PROBIT || pmod->ci == LOGIT ||
993 pmod->ci == HECKIT) && pmod->params != NULL) {
994 k = i;
995 } else if (pmod->ci == AR && i >= pmod->ncoeff) {
996 plain_ar_coeff_name(targ, pmod, i);
997 } else if (pmod->ci == ARCH && i >= pmod->ncoeff) {
998 sprintf(targ, "alpha(%d)", i - pmod->ncoeff);
999 } else if (pmod->list == NULL || j > pmod->list[0] ||
1000 pmod->list[j] >= dset->v) {
1001 k = i;
1002 } else {
1003 strcpy(targ, dset->varname[pmod->list[j]]);
1004 }
1005
1006 if (k >= 0) {
1007 if (pmod->params != NULL) {
1008 strcpy(targ, pmod->params[k]);
1009 } else {
1010 strcpy(targ, "unknown");
1011 }
1012 }
1013 }
1014
1015 return targ;
1016 }
1017
gretl_model_get_param_names(const MODEL * pmod,const DATASET * dset,int * err)1018 gretl_array *gretl_model_get_param_names (const MODEL *pmod,
1019 const DATASET *dset,
1020 int *err)
1021 {
1022 gretl_array *names = NULL;
1023
1024 if (pmod == NULL) {
1025 return NULL;
1026 }
1027
1028 names = gretl_array_new(GRETL_TYPE_STRINGS,
1029 pmod->ncoeff, err);
1030
1031 if (names != NULL) {
1032 char targ[VNAMELEN];
1033 int i;
1034
1035 for (i=0; i<pmod->ncoeff; i++) {
1036 gretl_model_get_param_name(pmod, dset, i, targ);
1037 gretl_array_set_string(names, i, targ, 1);
1038 }
1039 }
1040
1041 return names;
1042 }
1043
1044 /**
1045 * gretl_model_get_param_number:
1046 * @pmod: pointer to model.
1047 * @dset: dataset information.
1048 * @s: name of model parameter.
1049 *
1050 * Returns: the zero-based index of the coefficient in @pmod
1051 * corresponding to @s, or -1 if @s is not the name
1052 * of a parameter.
1053 */
1054
1055 int
gretl_model_get_param_number(const MODEL * pmod,const DATASET * dset,const char * s)1056 gretl_model_get_param_number (const MODEL *pmod, const DATASET *dset,
1057 const char *s)
1058 {
1059 char pname[VNAMELEN], tmp[VNAMELEN];
1060 int i, idx = -1;
1061
1062 if (pmod == NULL || s == NULL) {
1063 return -1;
1064 }
1065
1066 if (!strcmp(s, "0")) {
1067 strcpy(pname, "const");
1068 } else {
1069 strcpy(pname, s);
1070 }
1071
1072 for (i=0; i<pmod->ncoeff; i++) {
1073 gretl_model_get_param_name(pmod, dset, i, tmp);
1074 if (!strcmp(pname, tmp)) {
1075 idx = i;
1076 break;
1077 }
1078 }
1079
1080 return idx;
1081 }
1082
free_coeff_intervals(CoeffIntervals * cf)1083 void free_coeff_intervals (CoeffIntervals *cf)
1084 {
1085 int i;
1086
1087 free(cf->coeff);
1088 free(cf->maxerr);
1089
1090 if (cf->names != NULL) {
1091 for (i=0; i<cf->ncoeff; i++) {
1092 free(cf->names[i]);
1093 }
1094 free(cf->names);
1095 }
1096
1097 free(cf);
1098 }
1099
1100 /**
1101 * reset_coeff_intervals:
1102 * @cf: pointer to confidence intervals struct.
1103 * @alpha: nominal non-coverage, as decimal.
1104 *
1105 * Recomputes the intervals in @cf using the given value
1106 * of @alpha.
1107 *
1108 * Returns: 0 on success, non-zero if @alpha is out of
1109 * bounds.
1110 */
1111
reset_coeff_intervals(CoeffIntervals * cf,double alpha)1112 int reset_coeff_intervals (CoeffIntervals *cf, double alpha)
1113 {
1114 double se, tbak = cf->t;
1115 int i;
1116
1117 if (alpha <= 0 || alpha >= 1) {
1118 return E_DATA;
1119 }
1120
1121 if (cf->asy) {
1122 cf->t = normal_cdf_inverse(1 - alpha / 2);
1123 } else {
1124 cf->t = student_cdf_inverse(cf->df, 1 - alpha / 2);
1125 }
1126
1127 for (i=0; i<cf->ncoeff; i++) {
1128 if (cf->maxerr[i] > 0) {
1129 se = cf->maxerr[i] / tbak;
1130 cf->maxerr[i] = se * cf->t;
1131 }
1132 }
1133
1134 cf->alpha = alpha;
1135
1136 return 0;
1137 }
1138
1139 /**
1140 * gretl_model_get_coeff_intervals:
1141 * @pmod: pointer to gretl model.
1142 * @dset: dataset information.
1143 *
1144 * Save the 95 percent confidence intervals for the parameter
1145 * estimates in @pmod.
1146 *
1147 * Returns: pointer to #CONFINT struct containing the results.
1148 */
1149
1150 CoeffIntervals *
gretl_model_get_coeff_intervals(const MODEL * pmod,const DATASET * dset)1151 gretl_model_get_coeff_intervals (const MODEL *pmod,
1152 const DATASET *dset)
1153 {
1154 CoeffIntervals *cf;
1155 char pname[24];
1156 int i, err = 0;
1157
1158 cf = malloc(sizeof *cf);
1159 if (cf == NULL) {
1160 return NULL;
1161 }
1162
1163 cf->ncoeff = pmod->ncoeff;
1164 cf->df = pmod->dfd;
1165 cf->ifc = pmod->ifc;
1166
1167 cf->coeff = NULL;
1168 cf->maxerr = NULL;
1169 cf->names = NULL;
1170
1171 cf->coeff = malloc(cf->ncoeff * sizeof *cf->coeff);
1172 if (cf->coeff == NULL) {
1173 err = 1;
1174 goto bailout;
1175 }
1176
1177 cf->maxerr = malloc(cf->ncoeff * sizeof *cf->maxerr);
1178 if (cf->maxerr == NULL) {
1179 err = 1;
1180 goto bailout;
1181 }
1182
1183 cf->names = malloc(cf->ncoeff * sizeof *cf->names);
1184 if (cf->names == NULL) {
1185 err = 1;
1186 goto bailout;
1187 }
1188
1189 cf->alpha = .05;
1190
1191 if (ASYMPTOTIC_MODEL(pmod->ci)) {
1192 cf->asy = 1;
1193 cf->t = normal_cdf_inverse(0.975);
1194 } else {
1195 cf->asy = 0;
1196 cf->t = tcrit95(pmod->dfd);
1197 }
1198
1199 for (i=0; i<cf->ncoeff && !err; i++) {
1200 cf->coeff[i] = pmod->coeff[i];
1201 cf->maxerr[i] = (pmod->sderr[i] > 0)? cf->t * pmod->sderr[i] : 0.0;
1202 gretl_model_get_param_name(pmod, dset, i, pname);
1203 cf->names[i] = gretl_strdup(pname);
1204 if (cf->names[i] == NULL) {
1205 int j;
1206
1207 for (j=0; j<i; j++) {
1208 free(cf->names[i]);
1209 }
1210 free(cf->names);
1211 cf->names = NULL;
1212 err = 1;
1213 }
1214 }
1215
1216 bailout:
1217
1218 if (err) {
1219 free_coeff_intervals(cf);
1220 cf = NULL;
1221 }
1222
1223 return cf;
1224 }
1225
gretl_is_simple_OLS(const MODEL * pmod)1226 int gretl_is_simple_OLS (const MODEL *pmod)
1227 {
1228 if (pmod->ci == OLS &&
1229 pmod->list != NULL &&
1230 pmod->list[0] == 3 &&
1231 pmod->list[2] == 0) {
1232 return 1;
1233 } else {
1234 return 0;
1235 }
1236 }
1237
gretl_is_arima_model(const MODEL * pmod)1238 int gretl_is_arima_model (const MODEL *pmod)
1239 {
1240 int d = gretl_model_get_int(pmod, "arima_d");
1241 int D = gretl_model_get_int(pmod, "arima_D");
1242
1243 return (d > 0 || D > 0);
1244 }
1245
gretl_is_between_model(const MODEL * pmod)1246 int gretl_is_between_model (const MODEL *pmod)
1247 {
1248 if (pmod->ci == PANEL && (pmod->opt & OPT_B)) {
1249 return 1;
1250 } else {
1251 return 0;
1252 }
1253 }
1254
1255 /* "regular" here means pooled OLS, fixed effects or
1256 random effects
1257 */
1258
gretl_is_regular_panel_model(const MODEL * pmod)1259 int gretl_is_regular_panel_model (const MODEL *pmod)
1260 {
1261 if (pmod->ci == OLS && gretl_model_get_int(pmod, "pooled")) {
1262 return 1;
1263 } else if (pmod->ci != PANEL) {
1264 return 0;
1265 } else {
1266 gretlopt opt = pmod->opt & (OPT_P | OPT_F | OPT_U);
1267
1268 return opt != 0;
1269 }
1270 }
1271
1272 #define arma_included(m,i) (m == NULL || m[i] == '1')
1273
1274 /**
1275 * arma_model_nonseasonal_AR_order:
1276 * @pmod: pointer to gretl model.
1277 *
1278 * Returns: the non-seasonal autoregressive order of @pmod, or 0 if
1279 * @pmod is not an ARMA model.
1280 */
1281
arma_model_nonseasonal_AR_order(const MODEL * pmod)1282 int arma_model_nonseasonal_AR_order (const MODEL *pmod)
1283 {
1284 int p = 0;
1285
1286 if (pmod->ci == ARMA) {
1287 p = pmod->list[1];
1288 }
1289
1290 return p;
1291 }
1292
1293 /**
1294 * arma_model_nonseasonal_MA_order:
1295 * @pmod: pointer to gretl model.
1296 *
1297 * Returns: the non-seasonal moving-average order of @pmod, or 0 if
1298 * @pmod is not an ARMA model.
1299 */
1300
arma_model_nonseasonal_MA_order(const MODEL * pmod)1301 int arma_model_nonseasonal_MA_order (const MODEL *pmod)
1302 {
1303 int q = 0;
1304
1305 if (pmod->ci == ARMA) {
1306 if (gretl_is_arima_model(pmod)) {
1307 q = pmod->list[3];
1308 } else {
1309 q = pmod->list[2];
1310 }
1311 }
1312
1313 return q;
1314 }
1315
1316 /**
1317 * arma_model_max_AR_lag:
1318 * @pmod: pointer to gretl model.
1319 *
1320 * Returns: the maximum autoregressive lag in @pmod, or 0 if
1321 * @pmod is not an ARMA model. The maximum AR lag takes into
1322 * account any differencing (seasonal and/or non-seasonal) in
1323 * an ARIMA model.
1324 */
1325
arma_model_max_AR_lag(const MODEL * pmod)1326 int arma_model_max_AR_lag (const MODEL *pmod)
1327 {
1328 int pmax = 0;
1329
1330 if (pmod->ci == ARMA) {
1331 int p = arma_model_nonseasonal_AR_order(pmod);
1332 int P = gretl_model_get_int(pmod, "arma_P");
1333 int s = gretl_model_get_int(pmod, "arma_pd");
1334 int d = gretl_model_get_int(pmod, "arima_d");
1335 int D = gretl_model_get_int(pmod, "arima_D");
1336
1337 pmax = p + s * P;
1338 pmax += d + s * D;
1339 }
1340
1341 return pmax;
1342 }
1343
1344 /**
1345 * arma_model_max_MA_lag:
1346 * @pmod: pointer to gretl model.
1347 *
1348 * Returns: the maximum moving-average lag in @pmod, or 0 if
1349 * @pmod is not an ARMA model.
1350 */
1351
arma_model_max_MA_lag(const MODEL * pmod)1352 int arma_model_max_MA_lag (const MODEL *pmod)
1353 {
1354 int qmax = 0;
1355
1356 if (pmod->ci == ARMA) {
1357 int q, Q, pd;
1358
1359 q = arma_model_nonseasonal_MA_order(pmod);
1360 Q = gretl_model_get_int(pmod, "arma_Q");
1361
1362 if (Q == 0) {
1363 qmax = q;
1364 } else {
1365 pd = gretl_model_get_int(pmod, "arma_pd");
1366 qmax = Q * pd + q;
1367 }
1368 }
1369
1370 return qmax;
1371 }
1372
1373 /* From Box and Jenkins, 1976, pp 506-7, "Program 4": a clever
1374 algorithm for "unscrambling the coefficients", or in effect
1375 producing reduced-form AR coefficients that take into account any
1376 differencing, seasonal and/or non-seasonal.
1377 */
1378
ar_coeff_integrate(double * c0,int d,int D,int s,int pmax)1379 static int ar_coeff_integrate (double *c0, int d, int D, int s, int pmax)
1380 {
1381 int pstar = pmax + d + s * D;
1382 double *c1;
1383 int i, j, pp;
1384
1385 c1 = malloc((pstar + 1) * sizeof *c1);
1386 if (c1 == NULL) {
1387 return E_ALLOC;
1388 }
1389
1390 for (j=0; j<=pstar; j++) {
1391 c1[j] = 0.0;
1392 }
1393
1394 pp = pmax;
1395
1396 for (i=0; D>0 && i<D; i++) {
1397 for (j=0; j<=pstar; j++) {
1398 if (j < s) {
1399 c1[j] = c0[j];
1400 } else if (j <= pp) {
1401 c1[j] = c0[j] - c0[j-s];
1402 } else if (j <= pp + s) {
1403 c1[j] = -c0[j-s];
1404 }
1405 }
1406 pp += s;
1407 for (j=0; j<=pstar; j++) {
1408 c0[j] = c1[j];
1409 }
1410 }
1411
1412 for (i=0; d>0 && i<d; i++) {
1413 for (j=0; j<=pstar; j++) {
1414 if (j < 1) {
1415 c1[j] = c0[j];
1416 } else if (j <= pp) {
1417 c1[j] = c0[j] - c0[j-1];
1418 } else if (j <= pp + 1) {
1419 c1[j] = -c0[j-1];
1420 }
1421 }
1422 pp += 1;
1423 for (j=0; j<=pstar; j++) {
1424 c0[j] = c1[j];
1425 }
1426 }
1427
1428 free(c1);
1429
1430 return 0;
1431 }
1432
arma_included_lags(int k,const char * mask)1433 static int arma_included_lags (int k, const char *mask)
1434 {
1435 int i, nk = k;
1436
1437 if (mask != NULL) {
1438 nk = 0;
1439 for (i=0; i<k; i++) {
1440 if (mask[i] == '1') nk++;
1441 }
1442 }
1443
1444 return nk;
1445 }
1446
arma_AR_lags(const MODEL * pmod)1447 static int arma_AR_lags (const MODEL *pmod)
1448 {
1449 const char *pmask = gretl_model_get_data(pmod, "pmask");
1450 int p = arma_model_nonseasonal_AR_order(pmod);
1451
1452 return arma_included_lags(p, pmask);
1453 }
1454
arma_MA_lags(const MODEL * pmod)1455 static int arma_MA_lags (const MODEL *pmod)
1456 {
1457 const char *qmask = gretl_model_get_data(pmod, "qmask");
1458 int q = arma_model_nonseasonal_MA_order(pmod);
1459
1460 return arma_included_lags(q, qmask);
1461 }
1462
1463 /**
1464 * arma_model_AR_MA_coeffs:
1465 * @pmod: pointer to gretl model.
1466 * @phi_star: pointer to receive AR coeff vector.
1467 * @theta_star: pointer to receive MA coeff vector.
1468 * @opt: use OPT_I to get the "integrated" variant, which
1469 * handles differencing.
1470 *
1471 * Creates consolidated versions of the AR and MA coefficient vectors
1472 * from @pmod. If @pmod includes seasonal ARMA terms, the vectors are
1473 * suitably expanded to include the interactions between seasonal and
1474 * non-seasonal terms. If OPT_I is given and the dependent variable
1475 * has been differenced, the AR coefficients are integrated to account
1476 * for the differencing. These are the \Phi^* and \Theta^* as used by
1477 * Box and Jenkins for forecasting.
1478 *
1479 * Returns: 0 on success, non-zero on error.
1480 */
1481
arma_model_AR_MA_coeffs(const MODEL * pmod,gretl_vector ** phi_star,gretl_vector ** theta_star,gretlopt opt)1482 int arma_model_AR_MA_coeffs (const MODEL *pmod,
1483 gretl_vector **phi_star,
1484 gretl_vector **theta_star,
1485 gretlopt opt)
1486 {
1487 gretl_vector *ac = NULL;
1488 gretl_vector *mc = NULL;
1489 int *ainfo;
1490 int err = 0;
1491
1492 if (pmod->ci != ARMA) {
1493 err = E_DATA;
1494 } else if ((ainfo = gretl_model_get_data(pmod, "ainfo")) == NULL) {
1495 fprintf(stderr, "AR_MA_coeffs: no 'ainfo' available!\n");
1496 err = E_DATA;
1497 } else {
1498 const double *phi = NULL, *Phi = NULL;
1499 const double *theta = NULL, *Theta = NULL;
1500 const char *pmask = gretl_model_get_data(pmod, "pmask");
1501 const char *qmask = gretl_model_get_data(pmod, "qmask");
1502 int p = ainfo[1];
1503 int q = ainfo[2];
1504 int P = ainfo[3];
1505 int Q = ainfo[4];
1506 int np = ainfo[5];
1507 int nq = ainfo[6];
1508 int d = ainfo[7];
1509 int D = ainfo[8];
1510 int s = ainfo[9];
1511 int pmax = p + s * P;
1512 int pstar = pmax + d + s * D;
1513 int qmax = q + s * Q;
1514 double x, y;
1515 int i, j, k, ii;
1516
1517 if (pstar > 0) {
1518 /* we have some AR terms */
1519 ac = gretl_zero_matrix_new(pstar + 1, 1);
1520 if (ac == NULL) {
1521 err = E_ALLOC;
1522 }
1523 }
1524
1525 if (!err && qmax > 0) {
1526 /* we have some MA terms */
1527 mc = gretl_zero_matrix_new(qmax + 1, 1);
1528 if (mc == NULL) {
1529 gretl_vector_free(ac);
1530 ac = NULL;
1531 err = E_ALLOC;
1532 }
1533 }
1534
1535 if (!err) {
1536 phi = pmod->coeff + pmod->ifc; /* non-seasonal AR coeffs */
1537 Phi = phi + np; /* seasonal AR coeffs */
1538 theta = Phi + P; /* non-seasonal MA coeffs */
1539 Theta = theta + nq; /* seasonal MA coeffs */
1540 }
1541
1542 if (ac != NULL) {
1543 for (j=0; j<=P; j++) {
1544 x = (j == 0)? -1 : Phi[j-1];
1545 k = 0;
1546 for (i=0; i<=p; i++) {
1547 y = (i == 0)? -1 :
1548 (arma_included(pmask, i-1))? phi[k++] : 0;
1549 ii = i + s * j;
1550 ac->val[ii] -= x * y;
1551 }
1552 }
1553 if ((opt & OPT_I) && (D > 0 || d > 0)) {
1554 ar_coeff_integrate(ac->val, d, D, s, pmax);
1555 }
1556 }
1557
1558 if (mc != NULL) {
1559 for (j=0; j<=Q; j++) {
1560 x = (j == 0)? 1 : Theta[j-1];
1561 k = 0;
1562 for (i=0; i<=q; i++) {
1563 y = (i == 0)? 1 :
1564 (arma_included(qmask, i-1))? theta[k++] : 0;
1565 ii = i + s * j;
1566 mc->val[ii] += x * y;
1567 }
1568 }
1569 }
1570 }
1571
1572 if (!err) {
1573 *phi_star = ac;
1574 *theta_star = mc;
1575 }
1576
1577 return err;
1578 }
1579
1580 /*
1581 arma_model_spectrum computes
1582
1583 s2 C(exp(i*omega)) * C(exp(-i*omega))
1584 ---- * ----------------------------------
1585 2*pi A(exp(i*omega)) * A(exp(-i*omega))
1586
1587 for omega that goes from 0 to pi in T steps. The
1588 @phi matrix should contain the AR parameters and
1589 the @theta matrix the MA parameters.
1590 */
1591
arma_model_spectrum(gretl_matrix * phi,gretl_matrix * theta,double s2,int T,int * err)1592 static gretl_matrix *arma_model_spectrum (gretl_matrix *phi,
1593 gretl_matrix *theta,
1594 double s2, int T,
1595 int *err)
1596 {
1597 gretl_vector *ret = NULL;
1598 double num, den, nre_t, nim_t, dre_t, dim_t, xt;
1599 double c, s;
1600 double ar, ma;
1601 double scale;
1602 int nar, nma;
1603 int i, n, t;
1604
1605 ret = gretl_matrix_alloc(T, 2);
1606
1607 if (ret == NULL) {
1608 *err = E_ALLOC;
1609 return ret;
1610 }
1611
1612 nar = phi != NULL ? phi->rows : 0;
1613 nma = theta != NULL ? theta->rows : 0;
1614
1615 n = MAX(nar, nma);
1616 scale = s2/M_2PI;
1617
1618 for (t=0; t<T; t++) {
1619 xt = t * M_PI / (T-1);
1620 ar = nar > 0 ? gretl_vector_get(phi, 0) : 1.0;
1621 ma = nma > 0 ? gretl_vector_get(theta, 0) : 1.0;
1622 nre_t = ma * ma;
1623 dre_t = ar * ar;
1624 nim_t = dim_t = 0;
1625
1626 for (i=1; i<n; i++) {
1627 ar = i < nar ? gretl_vector_get(phi, i) : 0.0;
1628 ma = i < nma ? gretl_vector_get(theta, i) : 0.0;
1629 c = cos(i * xt);
1630 s = sin(i * xt);
1631 nre_t += c * ma;
1632 nim_t += s * ma;
1633 dre_t += c * ar;
1634 dim_t += s * ar;
1635 }
1636
1637 num = nre_t * nre_t + nim_t * nim_t;
1638 den = dre_t * dre_t + dim_t * dim_t;
1639 gretl_matrix_set(ret, t, 0, xt);
1640 gretl_matrix_set(ret, t, 1, scale * num/den);
1641 }
1642
1643 return ret;
1644 }
1645
arima_diff(const double * x,int t,int * delta,int k,int * err)1646 static double arima_diff (const double *x, int t,
1647 int *delta, int k, int *err)
1648 {
1649 double dxt = x[t];
1650 int i, p;
1651
1652 if (na(dxt)) {
1653 *err = E_MISSDATA;
1654 }
1655
1656 for (i=0; i<k && !*err; i++) {
1657 if (delta[i] != 0) {
1658 p = t - i - 1;
1659 if (p < 0 || na(x[p])) {
1660 dxt = NADBL;
1661 *err = E_MISSDATA;
1662 } else {
1663 dxt -= delta[i] * x[p];
1664 }
1665 }
1666 }
1667
1668 return dxt;
1669 }
1670
1671 /* get_arma_yvec: here we reconstruct the dependent variable from an
1672 ARMA model, differencing as we go in the ARIMA case. If the model
1673 includes any regressors besides the constant we have to net out
1674 their effect in order to produce a periodogram that is comparable
1675 with the ARMA spectrum. A complication is that in the ARIMAX case
1676 the regressors may also have to be differenced: this is indicated
1677 by the presence of the "xdiff" flag on the model.
1678 */
1679
get_arma_yvec(const MODEL * pmod,const DATASET * dset,int * err)1680 static gretl_vector *get_arma_yvec (const MODEL *pmod,
1681 const DATASET *dset,
1682 int *err)
1683 {
1684 gretl_vector *y = NULL;
1685 const double *beta = NULL;
1686 double yt, xti;
1687 int *xlist = NULL;
1688 int *delta = NULL;
1689 int yno = gretl_model_get_depvar(pmod);
1690 int T = pmod->nobs;
1691 int i, j, t, t1 = pmod->t1;
1692 int s, d, D, nx = 0;
1693 int k = 0, xdiff = 0;
1694
1695 if (yno < 1 || yno >= dset->v) {
1696 *err = E_DATA;
1697 return NULL;
1698 }
1699
1700 xlist = gretl_model_get_x_list(pmod);
1701
1702 if (xlist != NULL) {
1703 for (i=1; i<=xlist[0]; i++) {
1704 if (xlist[i] >= dset->v) {
1705 *err = E_DATA;
1706 break;
1707 } else if (xlist[i] != 0) {
1708 nx++;
1709 }
1710 }
1711 if (!*err && nx > 0) {
1712 beta = arma_model_get_x_coeffs(pmod);
1713 }
1714 }
1715
1716 if (*err) {
1717 free(xlist);
1718 return NULL;
1719 }
1720
1721 d = gretl_model_get_int(pmod, "arima_d");
1722 D = gretl_model_get_int(pmod, "arima_D");
1723
1724 if (d > 0 || D > 0) {
1725 delta = arima_delta_coeffs(d, D, dset->pd);
1726 if (delta == NULL) {
1727 *err = E_ALLOC;
1728 } else {
1729 k = d + dset->pd * D;
1730 xdiff = gretl_model_get_int(pmod, "xdiff");
1731 }
1732 }
1733
1734 if (!*err) {
1735 y = gretl_matrix_alloc(T, 1);
1736 if (y == NULL) {
1737 *err = E_ALLOC;
1738 }
1739 }
1740
1741 if (!*err) {
1742 for (t=0; t<T && !*err; t++) {
1743 s = t + t1;
1744 if (delta != NULL) {
1745 yt = arima_diff(dset->Z[yno], s, delta, k, err);
1746 } else {
1747 yt = dset->Z[yno][s];
1748 }
1749 if (beta != NULL) {
1750 /* ARMAX: net out X_t\beta */
1751 for (i=1, j=0; i<=xlist[0]; i++) {
1752 if (xlist[i] != 0) {
1753 if (xdiff) {
1754 xti = arima_diff(dset->Z[xlist[i]], s,
1755 delta, k, err);
1756 } else {
1757 xti = dset->Z[xlist[i]][s];
1758 }
1759 yt -= beta[j++] * xti;
1760 }
1761 }
1762 }
1763 gretl_vector_set(y, t, yt);
1764 }
1765 }
1766
1767 free(xlist);
1768 free(delta);
1769
1770 if (*err) {
1771 gretl_matrix_free(y);
1772 y = NULL;
1773 }
1774
1775 return y;
1776 }
1777
arma_spectrum_plot_data(const MODEL * pmod,const DATASET * dset,int * err)1778 gretl_matrix *arma_spectrum_plot_data (const MODEL *pmod,
1779 const DATASET *dset,
1780 int *err)
1781 {
1782 gretl_matrix *pdata = NULL;
1783 gretl_matrix *phi = NULL;
1784 gretl_matrix *theta = NULL;
1785 gretl_matrix *spec = NULL;
1786 gretl_matrix *pergm = NULL;
1787 gretl_matrix *y = NULL;
1788 int i, T = pmod->nobs;
1789 int grid = (T-1)/2;
1790 double s2 = pmod->sigma * pmod->sigma;
1791
1792 *err = arma_model_AR_MA_coeffs(pmod, &phi, &theta, OPT_NONE);
1793 if (*err) {
1794 return NULL;
1795 }
1796
1797 /* different sign convention from forecasting */
1798 if (phi != NULL) {
1799 gretl_matrix_multiply_by_scalar(phi, -1.0);
1800 }
1801
1802 spec = arma_model_spectrum(phi, theta, s2, grid, err);
1803 if (*err) {
1804 goto bailout;
1805 }
1806
1807 y = get_arma_yvec(pmod, dset, err);
1808
1809 if (!*err) {
1810 pergm = gretl_matrix_fft(y, 0, err);
1811 }
1812
1813 if (!*err) {
1814 pdata = gretl_matrix_alloc(grid, 4);
1815 if (pdata == NULL) {
1816 *err = E_ALLOC;
1817 } else {
1818 for (i=0; i<grid; i++) {
1819 gretl_matrix_set(pdata, i, 0, gretl_matrix_get(spec, i, 0));
1820 gretl_matrix_set(pdata, i, 1, gretl_matrix_get(spec, i, 1));
1821 gretl_matrix_set(pdata, i, 2, gretl_matrix_get(pergm, i+1, 0));
1822 gretl_matrix_set(pdata, i, 3, gretl_matrix_get(pergm, i+1, 1));
1823 }
1824 }
1825 }
1826
1827 bailout:
1828
1829 gretl_matrix_free(phi);
1830 gretl_matrix_free(theta);
1831 gretl_matrix_free(spec);
1832 gretl_matrix_free(pergm);
1833 gretl_matrix_free(y);
1834
1835 return pdata;
1836 }
1837
1838 /**
1839 * regarma_model_AR_coeffs:
1840 * @pmod: pointer to gretl model.
1841 * @phi0: pointer to receive AR coeff vector.
1842 * @pp: pointer to receive length of @phi0.
1843 *
1844 * Creates a consolidated version of the AR coefficients from @pmod.
1845 * If @pmod includes seasonal AR terms the vector is suitably expanded
1846 * to include the interactions between seasonal and non-seasonal
1847 * terms, but it is not integrated with respect to any differencing of
1848 * the dependent variable.
1849 *
1850 * Returns: 0 on success, non-zero on error.
1851 */
1852
regarma_model_AR_coeffs(const MODEL * pmod,double ** phi0,int * pp)1853 int regarma_model_AR_coeffs (const MODEL *pmod,
1854 double **phi0,
1855 int *pp)
1856 {
1857 const char *pmask = gretl_model_get_data(pmod, "pmask");
1858 int p = arma_model_nonseasonal_AR_order(pmod);
1859 int np = arma_included_lags(p, pmask);
1860 int P = gretl_model_get_int(pmod, "arma_P");
1861 int s = gretl_model_get_int(pmod, "arma_pd");
1862 const double *phi = NULL, *Phi = NULL;
1863 double *ac = NULL;
1864 double x, y;
1865 int i, j, k, ii, pmax;
1866 int err = 0;
1867
1868 pmax = p + s * P;
1869
1870 #if 0
1871 fprintf(stderr, "regarma_model_AR_coeffs: pmax = %d\n", pmax);
1872 #endif
1873
1874 if (pmax == 0) {
1875 *pp = 0;
1876 return 0;
1877 }
1878
1879 ac = malloc((pmax + 1) * sizeof *ac);
1880 if (ac == NULL) {
1881 return E_ALLOC;
1882 }
1883
1884 phi = pmod->coeff + pmod->ifc;
1885 Phi = phi + np;
1886
1887 for (i=0; i<=pmax; i++) {
1888 ac[i] = 0.0;
1889 }
1890
1891 for (j=0; j<=P; j++) {
1892 x = (j == 0)? -1 : Phi[j-1];
1893 k = 0;
1894 for (i=0; i<=p; i++) {
1895 y = (i == 0)? -1 :
1896 (arma_included(pmask, i-1))? phi[k++] : 0;
1897 ii = i + s * j;
1898 ac[ii] -= x * y;
1899 }
1900 }
1901
1902 *phi0 = ac;
1903 *pp = pmax;
1904
1905 return err;
1906 }
1907
1908 /**
1909 * arima_delta_coeffs:
1910 * @d: order of non-seasonal differencing (<= 2)
1911 * @D: order of seasonal differencing (<= 2)
1912 * @s: seasonal periodicity
1913 *
1914 * Returns: array of d + s * D coefficients of the lag operator in
1915 * the expansion of (1-L)^d * (1-L^s)^D. These are given in the
1916 * negative; for example, if d = 1 then c[0] = 1.
1917 */
1918
arima_delta_coeffs(int d,int D,int s)1919 int *arima_delta_coeffs (int d, int D, int s)
1920 {
1921 int i, k = d + s * D;
1922 int *c = malloc(k * sizeof *c);
1923
1924 if (c == NULL) {
1925 return NULL;
1926 }
1927
1928 for (i=0; i<k; i++) {
1929 c[i] = 0;
1930 }
1931
1932 if (d == 1) {
1933 c[0] = 1;
1934 } else if (d == 2) {
1935 c[0] = 2;
1936 c[1] = -1;
1937 }
1938
1939 if (D > 0) {
1940 c[s-1] += 1;
1941 if (d > 0) {
1942 c[s] -= 1;
1943 }
1944 if (d == 2) {
1945 c[s] -= 1;
1946 c[s+1] += 1;
1947 }
1948 if (D == 2) {
1949 c[s-1] += 1;
1950 c[2*s-1] -= 1;
1951 if (d > 0) {
1952 c[s] -= 1;
1953 c[2*s] += 1;
1954 }
1955 if (d == 2) {
1956 c[s] -= 1;
1957 c[2*s] += 1;
1958 c[s+1] += 1;
1959 c[2*s+1] -= 1;
1960 }
1961 }
1962 }
1963
1964 return c;
1965 }
1966
1967 /**
1968 * arma_model_get_x_coeffs:
1969 * @pmod: pointer to gretl model.
1970 *
1971 * Returns: pointer to the array of coefficients on the exogenous
1972 * regressors in @pmod, or %NULL if the model is not ARMA or if there
1973 * are no such regressors.
1974 */
1975
arma_model_get_x_coeffs(const MODEL * pmod)1976 const double *arma_model_get_x_coeffs (const MODEL *pmod)
1977 {
1978 const double *xc = NULL;
1979
1980 if (pmod->ci == ARMA && gretl_model_get_int(pmod, "armax")) {
1981 xc = pmod->coeff;
1982 xc += pmod->ifc;
1983 xc += arma_AR_lags(pmod);
1984 xc += arma_MA_lags(pmod);
1985 xc += gretl_model_get_int(pmod, "arma_P");
1986 xc += gretl_model_get_int(pmod, "arma_Q");
1987 }
1988
1989 return xc;
1990 }
1991
1992 /**
1993 * arma_model_get_n_arma_coeffs:
1994 * @pmod: pointer to gretl model.
1995 *
1996 * Returns: the sum of the numbers of AR and MA coefficients in
1997 * @pmod.
1998 */
1999
arma_model_get_n_arma_coeffs(const MODEL * pmod)2000 int arma_model_get_n_arma_coeffs (const MODEL *pmod)
2001 {
2002 int npq = 0;
2003
2004 if (pmod->ci == ARMA) {
2005 npq += arma_AR_lags(pmod);
2006 npq += arma_MA_lags(pmod);
2007 npq += gretl_model_get_int(pmod, "arma_P");
2008 npq += gretl_model_get_int(pmod, "arma_Q");
2009 }
2010
2011 return npq;
2012 }
2013
2014 /**
2015 * gretl_model_ahat_vec:
2016 * @pmod: pointer to gretl model.
2017 * @err: location to receive error code.
2018 *
2019 * If @pmod was estimated on panel data via fixed or random
2020 * effects, provides a column vector holding the "ahat" values
2021 * or estimated individual effects for the individuals
2022 * included in the active dataset when the model was estimated.
2023 * Some of these values may be NA if some units were not
2024 * actually included in estimation.
2025 *
2026 * Returns: allocated column vector on success or NULL on failure.
2027 */
2028
gretl_model_ahat_vec(const MODEL * pmod,int * err)2029 gretl_matrix *gretl_model_ahat_vec (const MODEL *pmod,
2030 int *err)
2031 {
2032 gretl_vector *a_vec = NULL;
2033 double *ahat;
2034 int i, s, t, N, T;
2035
2036 ahat = gretl_model_get_data(pmod, "ahat");
2037 T = gretl_model_get_int(pmod, "panel_T");
2038
2039 if (ahat == NULL || T == 0) {
2040 *err = E_BADSTAT;
2041 return NULL;
2042 }
2043
2044 /* number of units included in model dataset */
2045 N = pmod->full_n / T;
2046
2047 a_vec = gretl_column_vector_alloc(N);
2048 if (a_vec == NULL) {
2049 *err = E_ALLOC;
2050 return NULL;
2051 }
2052
2053 for (i=0; i<N; i++) {
2054 a_vec->val[i] = NADBL;
2055 for (t=0; t<T; t++) {
2056 s = i * T + t;
2057 if (!na(ahat[s])) {
2058 /* pick up first non-missing value */
2059 a_vec->val[i] = ahat[s];
2060 break;
2061 }
2062 }
2063 }
2064
2065 return a_vec;
2066 }
2067
arbond_get_depvar(const MODEL * pmod)2068 static int arbond_get_depvar (const MODEL *pmod)
2069 {
2070 int i;
2071
2072 for (i=1; i<=pmod->list[0]; i++) {
2073 if (pmod->list[i] == LISTSEP && i < pmod->list[0]) {
2074 return pmod->list[i+1];
2075 }
2076 }
2077
2078 return 0;
2079 }
2080
arma_depvar_pos(const MODEL * pmod)2081 static int arma_depvar_pos (const MODEL *pmod)
2082 {
2083 int sep1 = 0, sep2 = 0;
2084 int seasonal = 0;
2085 int arima = 0;
2086 int i, dvpos;
2087
2088 for (i=1; i<pmod->list[0]; i++) {
2089 if (pmod->list[i] == LISTSEP) {
2090 if (sep1 == 0) {
2091 sep1 = i;
2092 } else {
2093 sep2 = i;
2094 }
2095 }
2096 }
2097
2098 if (sep2) {
2099 dvpos = sep2 + 1;
2100 } else if (sep1) {
2101 dvpos = sep1 + 1;
2102 } else {
2103 if (gretl_model_get_int(pmod, "arma_P") ||
2104 gretl_model_get_int(pmod, "arima_D") ||
2105 gretl_model_get_int(pmod, "arma_Q")) {
2106 seasonal = 1;
2107 }
2108 if (gretl_model_get_int(pmod, "arima_d") ||
2109 gretl_model_get_int(pmod, "arima_D")) {
2110 arima = 1;
2111 }
2112 if (arima) {
2113 dvpos = (seasonal)? 9 : 5;
2114 } else {
2115 dvpos = (seasonal)? 7 : 4;
2116 }
2117 }
2118
2119 /* safety: should be impossible */
2120 if (dvpos == LISTSEP) {
2121 fprintf(stderr, "internal error in arma_depvar_pos\n");
2122 dvpos = 0;
2123 }
2124
2125 return dvpos;
2126 }
2127
2128 /**
2129 * gretl_model_get_depvar:
2130 * @pmod: pointer to gretl model.
2131 *
2132 * Returns: the ID number of the dependent variable in @pmod.
2133 */
2134
gretl_model_get_depvar(const MODEL * pmod)2135 int gretl_model_get_depvar (const MODEL *pmod)
2136 {
2137 int dv = gretl_model_get_int(pmod, "yno");
2138
2139 if (dv > 0) {
2140 return dv;
2141 }
2142
2143 if (pmod != NULL && pmod->list != NULL) {
2144 if (pmod->ci == GARCH) {
2145 dv = pmod->list[4];
2146 } else if (pmod->ci == ARMA) {
2147 dv = pmod->list[arma_depvar_pos(pmod)];
2148 } else if (pmod->ci == DPANEL) {
2149 dv = arbond_get_depvar(pmod);
2150 } else {
2151 dv = pmod->list[1];
2152 }
2153 }
2154
2155 return dv;
2156 }
2157
2158 /**
2159 * gretl_model_get_depvar_name:
2160 * @pmod: pointer to gretl model.
2161 * @dset: dataset information.
2162 *
2163 * Returns: the name of the dependent variable in @pmod.
2164 */
2165
gretl_model_get_depvar_name(const MODEL * pmod,const DATASET * dset)2166 const char *gretl_model_get_depvar_name (const MODEL *pmod,
2167 const DATASET *dset)
2168 {
2169 int dv;
2170
2171 if (pmod->depvar != NULL) {
2172 return pmod->depvar;
2173 }
2174
2175 dv = gretl_model_get_int(pmod, "yno");
2176
2177 if (dv == 0) {
2178 if (pmod != NULL && pmod->list != NULL) {
2179 if (pmod->ci == GARCH) {
2180 dv = pmod->list[4];
2181 } else if (pmod->ci == ARMA) {
2182 dv = pmod->list[arma_depvar_pos(pmod)];
2183 } else if (pmod->ci == DPANEL) {
2184 dv = arbond_get_depvar(pmod);
2185 } else {
2186 dv = pmod->list[1];
2187 }
2188 }
2189 }
2190
2191 return dset->varname[dv];
2192 }
2193
ordered_model(const MODEL * pmod)2194 static int ordered_model (const MODEL *pmod)
2195 {
2196 if (pmod->ci == PROBIT || pmod->ci == LOGIT) {
2197 if (gretl_model_get_int(pmod, "ordered")) {
2198 return 1;
2199 }
2200 }
2201
2202 return 0;
2203 }
2204
longlist_model(const MODEL * pmod)2205 static int longlist_model (const MODEL *pmod)
2206 {
2207 if (pmod->ci == LOGIT ||
2208 pmod->ci == NEGBIN ||
2209 pmod->ci == DURATION ||
2210 pmod->ci == PANEL) {
2211 return 1;
2212 } else if (pmod->ci == PROBIT && (pmod->opt & OPT_E)) {
2213 /* random-effects probit */
2214 return 1;
2215 } else {
2216 return 0;
2217 }
2218 }
2219
2220 /**
2221 * gretl_model_get_x_list:
2222 * @pmod: pointer to gretl model.
2223 *
2224 * Returns: an allocated copy of the list of independent
2225 * variables included in @pmod, or %NULL on failure.
2226 */
2227
gretl_model_get_x_list(const MODEL * pmod)2228 int *gretl_model_get_x_list (const MODEL *pmod)
2229 {
2230 int *list = NULL;
2231 int i, nx;
2232
2233 if (pmod == NULL || pmod->errcode || pmod->list == NULL) {
2234 return NULL;
2235 }
2236
2237 if (pmod->ci == ARMA) {
2238 int start = arma_depvar_pos(pmod);
2239
2240 nx = pmod->list[0] - start + pmod->ifc;
2241 if (nx > 0) {
2242 list = gretl_list_new(nx);
2243 if (list != NULL) {
2244 if (pmod->ifc) {
2245 list[1] = 0;
2246 for (i=2; i<=list[0]; i++) {
2247 list[i] = pmod->list[i + start - 1];
2248 }
2249 } else {
2250 for (i=1; i<=list[0]; i++) {
2251 list[i] = pmod->list[i + start];
2252 }
2253 }
2254 }
2255 }
2256 } else if (pmod->ci == GARCH) {
2257 nx = pmod->list[0] - 4;
2258 if (nx > 0) {
2259 list = gretl_list_new(nx);
2260 if (list != NULL) {
2261 for (i=1; i<=list[0]; i++) {
2262 list[i] = pmod->list[i + 4];
2263 }
2264 }
2265 }
2266 } else if (pmod->ci == DPANEL) {
2267 int ifc = gretl_model_get_int(pmod, "ifc");
2268 int vi, sep = 0;
2269
2270 nx = 0;
2271 for (i=2; i<=pmod->list[0]; i++) {
2272 if (pmod->list[i] == LISTSEP) {
2273 sep++;
2274 if (sep == 1) {
2275 i += 2;
2276 } else if (sep == 2) {
2277 break;
2278 }
2279 }
2280 if (sep == 1 && i <= pmod->list[0]) {
2281 vi = pmod->list[i];
2282 if (vi > 0 || ifc) {
2283 /* don't include const if it was dropped */
2284 list = gretl_list_append_term(&list, vi);
2285 if (list == NULL) {
2286 return NULL;
2287 }
2288 }
2289 }
2290 }
2291 } else if (pmod->ci == BIPROBIT) {
2292 /* not sure what we ought to do here, but for now we'll
2293 return the list of regressors for the first equation
2294 */
2295 nx = 0;
2296 for (i=3; i<=pmod->list[0]; i++) {
2297 if (pmod->list[i] == LISTSEP) {
2298 nx = i - 3;
2299 break;
2300 }
2301 }
2302 if (nx > 0) {
2303 list = gretl_list_new(nx);
2304 if (list != NULL) {
2305 for (i=1; i<=nx; i++) {
2306 list[i] = pmod->list[i+2];
2307 }
2308 }
2309 }
2310 } else if (pmod->ci == MIDASREG) {
2311 int *lfx = gretl_model_get_list(pmod, "lfxlist");
2312
2313 if (lfx != NULL) {
2314 list = gretl_list_copy(lfx);
2315 }
2316 } else if (ordered_model(pmod)) {
2317 nx = pmod->list[0] - 1;
2318 list = gretl_list_new(nx);
2319 if (list != NULL) {
2320 for (i=1; i<=list[0]; i++) {
2321 list[i] = pmod->list[i + 1];
2322 }
2323 }
2324 } else if (!NONLIST_MODEL(pmod->ci)) {
2325 if (pmod->ci == HECKIT) {
2326 nx = gretl_model_get_int(pmod, "base-coeffs");
2327 } else if (longlist_model(pmod)) {
2328 /* models in which the array of coefficients
2329 is (or may be) longer than the list of
2330 regressors
2331 */
2332 nx = pmod->list[0] - 1;
2333 } else {
2334 nx = pmod->ncoeff;
2335 }
2336 if (nx > 0) {
2337 if (nx > pmod->list[0] - 1) {
2338 /* don't read off the end of pmod->list! */
2339 nx = pmod->list[0] - 1;
2340 }
2341 list = gretl_list_new(nx);
2342 if (list != NULL) {
2343 for (i=1; i<=list[0]; i++) {
2344 list[i] = pmod->list[i + 1];
2345 }
2346 }
2347 }
2348 }
2349
2350 return list;
2351 }
2352
2353 /**
2354 * gretl_model_get_secondary_list:
2355 * @pmod: model to examine.
2356 *
2357 * Retrieve an allocated copy of the secondary list from
2358 * @pmod: e.g. the list of instruments in the case of %IVREG, or
2359 * the selection equation list for %HECKIT.
2360 *
2361 * Returns: allocated list or %NULL on error.
2362 */
2363
gretl_model_get_secondary_list(const MODEL * pmod)2364 int *gretl_model_get_secondary_list (const MODEL *pmod)
2365 {
2366 int *list = NULL;
2367 int pos;
2368
2369 pos = gretl_list_separator_position(pmod->list);
2370
2371 if (pos > 0) {
2372 list = gretl_list_copy_from_pos(pmod->list, pos + 1);
2373 }
2374
2375 return list;
2376 }
2377
2378 /**
2379 * gretl_model_get_y_list:
2380 * @pmod: model to examine.
2381 *
2382 * Retrieve an allocated copy of the list of dependent variables
2383 * for @pmod: in almost all cases this will have a single
2384 * element; an exception is biprobit.
2385 *
2386 * Returns: allocated list or %NULL on error.
2387 */
2388
gretl_model_get_y_list(const MODEL * pmod)2389 int *gretl_model_get_y_list (const MODEL *pmod)
2390 {
2391 int *list = NULL;
2392
2393 if (pmod->list == NULL) {
2394 return NULL;
2395 }
2396
2397 if (pmod->ci == BIPROBIT) {
2398 list = gretl_list_new(2);
2399 if (list != NULL) {
2400 list[1] = pmod->list[1];
2401 list[2] = pmod->list[2];
2402 }
2403 } else {
2404 list = gretl_list_new(1);
2405 if (list != NULL) {
2406 list[1] = gretl_model_get_depvar(pmod);
2407 }
2408 }
2409
2410 return list;
2411 }
2412
2413 /**
2414 * gretl_model_allocate_storage:
2415 * @pmod: pointer to model.
2416 *
2417 * Allocates space for coefficients and standard errors,
2418 * residuals and fitted values in @pmod. The sizes of
2419 * the arrays are based on the %ncoeff and %full_n
2420 * members of @pmod, which must be set first. The
2421 * residuals and fitted values are initialized to
2422 * gretl's missing value.
2423 *
2424 * Returns: 0 on success, %E_ALLOC on error.
2425 */
2426
gretl_model_allocate_storage(MODEL * pmod)2427 int gretl_model_allocate_storage (MODEL *pmod)
2428 {
2429 int k = pmod->ncoeff;
2430 int T = pmod->full_n;
2431
2432 if (k > 0) {
2433 int i;
2434
2435 pmod->coeff = malloc(k * sizeof *pmod->coeff);
2436 if (pmod->coeff == NULL) {
2437 return E_ALLOC;
2438 }
2439 pmod->sderr = malloc(k * sizeof *pmod->sderr);
2440 if (pmod->sderr == NULL) {
2441 return E_ALLOC;
2442 }
2443 for (i=0; i<k; i++) {
2444 pmod->coeff[i] = pmod->sderr[i] = NADBL;
2445 }
2446 }
2447
2448 if (T > 0) {
2449 int t;
2450
2451 pmod->uhat = malloc(T * sizeof *pmod->uhat);
2452 if (pmod->uhat == NULL) {
2453 return E_ALLOC;
2454 }
2455 pmod->yhat = malloc(T * sizeof *pmod->yhat);
2456 if (pmod->yhat == NULL) {
2457 return E_ALLOC;
2458 }
2459 for (t=0; t<T; t++) {
2460 pmod->uhat[t] = pmod->yhat[t] = NADBL;
2461 }
2462 }
2463
2464 return 0;
2465 }
2466
2467 /**
2468 * gretl_model_new_vcv:
2469 * @pmod: pointer to model.
2470 * @nelem: pointer to receive number of elements in
2471 * the packed array, or %NULL;
2472 *
2473 * Allocates space for a packed coefficient covariance matrix
2474 * in @pmod (if such space is not already allocated). Sets
2475 * all entries in the array to zero.
2476 *
2477 * Returns: 0 on success, %E_ALLOC on error.
2478 */
2479
gretl_model_new_vcv(MODEL * pmod,int * nelem)2480 int gretl_model_new_vcv (MODEL *pmod, int *nelem)
2481 {
2482 int nv = pmod->ncoeff;
2483 int nxpx = (nv * nv + nv) / 2;
2484 int i, err = 0;
2485
2486 if (pmod->vcv == NULL) {
2487 /* not already allocated */
2488 pmod->vcv = malloc(nxpx * sizeof *pmod->vcv);
2489 if (pmod->vcv == NULL) {
2490 err = E_ALLOC;
2491 }
2492 }
2493
2494 if (pmod->vcv != NULL) {
2495 for (i=0; i<nxpx; i++) {
2496 pmod->vcv[i] = 0.0;
2497 }
2498 if (nelem != NULL) {
2499 *nelem = nxpx;
2500 }
2501 }
2502
2503 return err;
2504 }
2505
vcv_get_se(double vii,int restricted)2506 static double vcv_get_se (double vii, int restricted)
2507 {
2508 if (restricted) {
2509 vii = fabs(vii) < 1.0e-17 ? 0.0 : vii;
2510 }
2511
2512 return (na(vii) || vii < 0)? NADBL : sqrt(vii);
2513 }
2514
2515 /**
2516 * gretl_model_write_vcv:
2517 * @pmod: pointer to model.
2518 * @V: covariance matrix.
2519 *
2520 * Write the covariance matrix @V into the model @pmod, using the
2521 * special packed format that is required by the MODEL struct,
2522 * and set the standard errors to the square root of the diagonal
2523 * elements of this matrix.
2524 *
2525 * Returns: 0 on success, non-zero code on error.
2526 */
2527
gretl_model_write_vcv(MODEL * pmod,const gretl_matrix * V)2528 int gretl_model_write_vcv (MODEL *pmod, const gretl_matrix *V)
2529 {
2530 int i, j, k, n;
2531 double x, *tmp;
2532 int err = 0;
2533
2534 if (gretl_is_null_matrix(V)) {
2535 return 0; /* no-op */
2536 }
2537
2538 if (V->cols != V->rows) {
2539 return E_NONCONF;
2540 }
2541
2542 k = V->rows;
2543 n = (k * k + k) / 2;
2544
2545 /* reallocate model vcv in case it's wrongly sized */
2546 tmp = realloc(pmod->vcv, n * sizeof *tmp);
2547 if (tmp == NULL) {
2548 err = E_ALLOC;
2549 } else {
2550 pmod->vcv = tmp;
2551 }
2552
2553 /* same for standard errors array */
2554 if (!err) {
2555 tmp = realloc(pmod->sderr, k * sizeof *tmp);
2556 if (tmp == NULL) {
2557 err = E_ALLOC;
2558 } else {
2559 pmod->sderr = tmp;
2560 }
2561 }
2562
2563 if (!err) {
2564 int restricted, idx = 0;
2565
2566 restricted = gretl_model_get_int(pmod, "restricted");
2567
2568 for (i=0; i<k; i++) {
2569 for (j=i; j<k; j++) {
2570 x = gretl_matrix_get(V, i, j);
2571 pmod->vcv[idx++] = x;
2572 if (i == j) {
2573 pmod->sderr[i] = vcv_get_se(x, restricted);
2574 }
2575 }
2576 }
2577 }
2578
2579 return err;
2580 }
2581
make_cluster_vals(MODEL * pmod,int cvar,const DATASET * dset,int * err)2582 static gretl_matrix *make_cluster_vals (MODEL *pmod, int cvar,
2583 const DATASET *dset,
2584 int *err)
2585 {
2586 gretl_matrix *cvals = NULL;
2587 double ct, *cdata;
2588 int s, t;
2589
2590 cdata = malloc(pmod->nobs * sizeof *cdata);
2591
2592 if (cdata == NULL) {
2593 *err = E_ALLOC;
2594 } else {
2595 s = 0;
2596 for (t=pmod->t1; t<=pmod->t2 && !*err; t++) {
2597 if (!model_missing(pmod, t)) {
2598 ct = dset->Z[cvar][t];
2599 if (na(ct)) {
2600 *err = E_MISSDATA;
2601 } else {
2602 cdata[s++] = ct;
2603 }
2604 }
2605 }
2606 }
2607
2608 if (!*err) {
2609 cvals = gretl_matrix_values(cdata, pmod->nobs, OPT_S, err);
2610 if (!*err && gretl_vector_get_length(cvals) < 2) {
2611 gretl_matrix_free(cvals);
2612 cvals = NULL;
2613 *err = E_DATA;
2614 }
2615 }
2616
2617 free(cdata);
2618
2619 return cvals;
2620 }
2621
model_make_clustered_GG(MODEL * pmod,int ci,const gretl_matrix * G,gretl_matrix * GG,const DATASET * dset,int * pcvar,int * pnc)2622 static int model_make_clustered_GG (MODEL *pmod, int ci,
2623 const gretl_matrix *G,
2624 gretl_matrix *GG,
2625 const DATASET *dset,
2626 int *pcvar,
2627 int *pnc)
2628 {
2629 gretl_matrix *cvals = NULL;
2630 gretl_matrix *Gi = NULL;
2631 const double *cZ = NULL;
2632 const char *cname;
2633 int cvar, n_c = 0;
2634 int k = G->cols;
2635 int i, j, t;
2636 int err = 0;
2637
2638 cname = get_optval_string(ci, OPT_C);
2639 if (cname == NULL) {
2640 return E_PARSE;
2641 }
2642
2643 cvar = current_series_index(dset, cname);
2644
2645 if (cvar < 1 || cvar >= dset->v) {
2646 err = E_UNKVAR;
2647 } else {
2648 cvals = make_cluster_vals(pmod, cvar, dset, &err);
2649 }
2650
2651 if (!err) {
2652 Gi = gretl_matrix_alloc(1, k);
2653 if (Gi == NULL) {
2654 gretl_matrix_free(cvals);
2655 err = E_ALLOC;
2656 }
2657 }
2658
2659 if (err) {
2660 return err;
2661 }
2662
2663 n_c = gretl_vector_get_length(cvals);
2664 cZ = dset->Z[cvar];
2665
2666 for (i=0; i<n_c; i++) {
2667 double cvi = cvals->val[i];
2668 int s = 0;
2669
2670 gretl_matrix_zero(Gi);
2671
2672 for (t=pmod->t1; t<=pmod->t2; t++) {
2673 if (!model_missing(pmod, t)) {
2674 if (cZ[t] == cvi) {
2675 for (j=0; j<k; j++) {
2676 Gi->val[j] += gretl_matrix_get(G, s, j);
2677 }
2678 }
2679 s++;
2680 }
2681 }
2682
2683 gretl_matrix_multiply_mod(Gi, GRETL_MOD_TRANSPOSE,
2684 Gi, GRETL_MOD_NONE,
2685 GG, GRETL_MOD_CUMULATE);
2686 }
2687
2688 if (!err) {
2689 *pcvar = cvar;
2690 *pnc = n_c;
2691 }
2692
2693 gretl_matrix_free(cvals);
2694 gretl_matrix_free(Gi);
2695
2696 return err;
2697 }
2698
do_biprobit_adjustment(MODEL * pmod,gretl_matrix * V)2699 static int do_biprobit_adjustment (MODEL *pmod,
2700 gretl_matrix *V)
2701 {
2702 void (*biprobit_adj) (MODEL *, gretl_matrix *);
2703
2704 biprobit_adj = get_plugin_function("biprobit_adjust_vcv");
2705
2706 if (biprobit_adj == NULL) {
2707 return E_FOPEN;
2708 } else {
2709 (*biprobit_adj) (pmod, V);
2710 return 0;
2711 }
2712 }
2713
2714 /**
2715 * gretl_model_add_QML_vcv:
2716 * @pmod: pointer to model.
2717 * @ci: command index for model.
2718 * @H: inverse of the (negative) Hessian, k x k.
2719 * @G: score matrix, T x k.
2720 * @dset: pointer to dataset (can be NULL if not doing
2721 * clustering).
2722 * @opt: may include OPT_C for cluster-robust variant.
2723 * @pV: location to receive the full covariance matrix, or NULL.
2724 *
2725 * Write a QML covariance matrix into the model @pmod, and set
2726 * the standard errors to the square root of the diagonal
2727 * elements of this matrix. The @ci argument, specifying the
2728 * estimator for @pmod, is required only if @opt includes
2729 * OPT_C; otherwise it is ignored.
2730 *
2731 * Returns: 0 on success, non-zero code on error.
2732 */
2733
gretl_model_add_QML_vcv(MODEL * pmod,int ci,const gretl_matrix * H,const gretl_matrix * G,const DATASET * dset,gretlopt opt,gretl_matrix ** pV)2734 int gretl_model_add_QML_vcv (MODEL *pmod, int ci,
2735 const gretl_matrix *H,
2736 const gretl_matrix *G,
2737 const DATASET *dset,
2738 gretlopt opt,
2739 gretl_matrix **pV)
2740 {
2741 gretl_matrix *GG = NULL;
2742 gretl_matrix *V = NULL;
2743 int cvar = 0, n_c = 0;
2744 int k = H->rows;
2745 int err = 0;
2746
2747 V = gretl_matrix_alloc(k, k);
2748
2749 if (V == NULL) {
2750 return E_ALLOC;
2751 }
2752
2753 if (!err) {
2754 if (opt & OPT_C) {
2755 /* clustered */
2756 GG = gretl_zero_matrix_new(k, k);
2757 if (GG == NULL) {
2758 err = E_ALLOC;
2759 } else {
2760 err = model_make_clustered_GG(pmod, ci, G, GG,
2761 dset, &cvar, &n_c);
2762 }
2763 } else if (opt & OPT_N) {
2764 /* Newey-West */
2765 GG = newey_west_OPG(G, &err);
2766 } else {
2767 /* regular QML using OPG */
2768 GG = gretl_matrix_XTX_new(G);
2769 if (GG == NULL) {
2770 err = E_ALLOC;
2771 }
2772 }
2773 }
2774
2775 if (!err) {
2776 err = gretl_matrix_qform(H, GRETL_MOD_NONE, GG,
2777 V, GRETL_MOD_NONE);
2778 if (!err && (opt & OPT_C)) {
2779 /* clustering: use stata-style df adjustment */
2780 double dfc = n_c / (n_c - 1.0);
2781
2782 gretl_matrix_multiply_by_scalar(V, dfc);
2783 }
2784 }
2785
2786 if (!err) {
2787 if (ci == BIPROBIT) {
2788 do_biprobit_adjustment(pmod, V);
2789 }
2790 err = gretl_model_write_vcv(pmod, V);
2791 }
2792
2793 if (!err) {
2794 if (opt & OPT_C) {
2795 gretl_model_set_int(pmod, "n_clusters", n_c);
2796 gretl_model_set_vcv_info(pmod, VCV_CLUSTER, cvar);
2797 pmod->opt |= OPT_C;
2798 } else {
2799 MLVCVType vt = (opt & OPT_N)? ML_HAC : ML_QML;
2800
2801 gretl_model_set_vcv_info(pmod, VCV_ML, vt);
2802 }
2803 pmod->opt |= OPT_R;
2804 }
2805
2806 gretl_matrix_free(GG);
2807
2808 if (pV != NULL) {
2809 *pV = V;
2810 } else {
2811 gretl_matrix_free(V);
2812 }
2813
2814 return err;
2815 }
2816
2817 /**
2818 * gretl_model_add_hessian_vcv:
2819 * @pmod: pointer to model.
2820 * @H: inverse of the (negative) Hessian.
2821 *
2822 * Write @H into the model @pmod as its covariance matrix, and
2823 * set the standard errors to the square roots of the diagonal
2824 * elements of this matrix.
2825 *
2826 * Returns: 0 on success, non-zero code on error.
2827 */
2828
gretl_model_add_hessian_vcv(MODEL * pmod,const gretl_matrix * H)2829 int gretl_model_add_hessian_vcv (MODEL *pmod,
2830 const gretl_matrix *H)
2831 {
2832 int err = gretl_model_write_vcv(pmod, H);
2833
2834 if (!err) {
2835 gretl_model_set_vcv_info(pmod, VCV_ML, ML_HESSIAN);
2836 }
2837
2838 return err;
2839 }
2840
2841 /**
2842 * gretl_model_add_OPG_vcv:
2843 * @pmod: pointer to model.
2844 * @G: T x k gradient matrix.
2845 * @pV: location to receive the full covariance matrix, or NULL.
2846 *
2847 * Compute (G'G)^{-1}, write this into @pmod as its covariance matrix,
2848 * and set the standard errors to the square roots of the diagonal
2849 * elements of this matrix.
2850 *
2851 * Returns: 0 on success, non-zero code on error.
2852 */
2853
gretl_model_add_OPG_vcv(MODEL * pmod,const gretl_matrix * G,gretl_matrix ** pV)2854 int gretl_model_add_OPG_vcv (MODEL *pmod,
2855 const gretl_matrix *G,
2856 gretl_matrix **pV)
2857 {
2858 gretl_matrix *GG = gretl_matrix_XTX_new(G);
2859 int err = 0;
2860
2861 if (GG == NULL) {
2862 return E_ALLOC;
2863 }
2864
2865 err = gretl_invert_symmetric_matrix(GG);
2866
2867 if (!err) {
2868 err = gretl_model_write_vcv(pmod, GG);
2869 if (!err) {
2870 gretl_model_set_vcv_info(pmod, VCV_ML, ML_OP);
2871 }
2872 } else {
2873 /* diagnostic: recreate and print G'G */
2874 gretl_matrix_multiply_mod(G, GRETL_MOD_TRANSPOSE,
2875 G, GRETL_MOD_NONE,
2876 GG, GRETL_MOD_NONE);
2877 gretl_matrix_print(GG, "non-p.d. G'G");
2878 }
2879
2880 if (pV != NULL) {
2881 *pV = GG;
2882 } else {
2883 gretl_matrix_free(GG);
2884 }
2885
2886 return err;
2887 }
2888
copy_vcv_subset(const MODEL * pmod)2889 static double *copy_vcv_subset (const MODEL *pmod)
2890 {
2891 double *V;
2892 int nc = pmod->ncoeff;
2893 int k = pmod->list[0] - 1;
2894 int n = k * (k + 1) / 2;
2895 int i, j;
2896
2897 V = malloc(n * sizeof *V);
2898 if (V == NULL) {
2899 return NULL;
2900 }
2901
2902 for (i=0; i<k; i++) {
2903 for (j=0; j<=i; j++) {
2904 V[ijton(i, j, k)] = pmod->vcv[ijton(i, j, nc)];
2905 }
2906 }
2907
2908 return V;
2909 }
2910
2911 /**
2912 * gretl_model_get_vcv:
2913 * @pmod: pointer to model.
2914 * @dset: dataset information.
2915 *
2916 * Supplies the caller with a copy of the variance-covariance
2917 * matrix for the parameter estimates in @pmod, in a format
2918 * suitable for printing. See also free_vcv(). To get the
2919 * covariance matrix as a gretl_matrix, see
2920 * gretl_vcv_matrix_from_model().
2921 *
2922 * Returns: #VMatrix struct or %NULL on error.
2923 */
2924
gretl_model_get_vcv(MODEL * pmod,const DATASET * dset)2925 VMatrix *gretl_model_get_vcv (MODEL *pmod, const DATASET *dset)
2926 {
2927 char varname[VNAMELEN];
2928 int i, k = pmod->ncoeff;
2929 int special = 0;
2930 VMatrix *vcv;
2931
2932 vcv = vmatrix_new();
2933
2934 if (vcv == NULL) {
2935 return NULL;
2936 }
2937
2938 /* special for fixed effects panel model: strip out the
2939 vcv elements for per-unit dummies, if present
2940 */
2941 if (pmod->ci == PANEL) {
2942 int k2 = pmod->list[0] - 1;
2943
2944 if (k > k2) {
2945 k = k2;
2946 special = 1;
2947 }
2948 }
2949
2950 vcv->names = strings_array_new(k);
2951 if (vcv->names == NULL) {
2952 free(vcv);
2953 return NULL;
2954 }
2955
2956 for (i=0; i<k; i++) {
2957 gretl_model_get_param_name(pmod, dset, i, varname);
2958 vcv->names[i] = gretl_strdup(varname);
2959 if (vcv->names[i] == NULL) {
2960 free_vmatrix(vcv);
2961 return NULL;
2962 }
2963 }
2964
2965 if (pmod->vcv == NULL && makevcv(pmod, pmod->sigma)) {
2966 free_vmatrix(vcv);
2967 return NULL;
2968 }
2969
2970 if (special) {
2971 /* copy subset of vcv */
2972 vcv->vec = copy_vcv_subset(pmod);
2973 } else {
2974 /* copy full vcv */
2975 vcv->vec = copyvec(pmod->vcv, k * (k + 1) / 2);
2976 }
2977
2978 if (vcv->vec == NULL) {
2979 free_vmatrix(vcv);
2980 return NULL;
2981 }
2982
2983 vcv->ci = pmod->ci;
2984 vcv->dim = k;
2985 vcv->t1 = pmod->t1;
2986 vcv->t2 = pmod->t2;
2987
2988 return vcv;
2989 }
2990
2991 /**
2992 * gretl_model_get_vcv_element:
2993 * @pmod: pointer to model.
2994 * @i: row (0-based).
2995 * @j: column(0-based).
2996 * @np: the number of parameters represented in @pmod's
2997 * covariance matrix.
2998 *
2999 * Returns: the (@i, @j) element of @pmod's covariance matrix,
3000 * or #NADBL on failure.
3001 */
3002
gretl_model_get_vcv_element(const MODEL * pmod,int i,int j,int np)3003 double gretl_model_get_vcv_element (const MODEL *pmod,
3004 int i, int j,
3005 int np)
3006 {
3007 if (pmod->vcv == NULL) {
3008 return NADBL;
3009 } else {
3010 int k = ijton(i, j, np);
3011
3012 return pmod->vcv[k];
3013 }
3014 }
3015
3016 /**
3017 * gretl_model_write_coeffs:
3018 * @pmod: pointer to model.
3019 * @b: array of coefficients.
3020 * @k: number of elements in @b.
3021 *
3022 * Write the coefficients @b into the model @pmod, whose
3023 * coefficient array is resized appropriately if need be.
3024 *
3025 * Returns: 0 on success, non-zero code on error.
3026 */
3027
gretl_model_write_coeffs(MODEL * pmod,double * b,int k)3028 int gretl_model_write_coeffs (MODEL *pmod, double *b, int k)
3029 {
3030 size_t sz = k * sizeof(double);
3031 int err = 0;
3032
3033 if (pmod->coeff == NULL || pmod->ncoeff != k) {
3034 double *tmp = realloc(pmod->coeff, sz);
3035
3036 if (tmp == NULL) {
3037 err = E_ALLOC;
3038 } else {
3039 pmod->coeff = tmp;
3040 }
3041 }
3042
3043 if (!err) {
3044 memcpy(pmod->coeff, b, sz);
3045 pmod->ncoeff = k;
3046 }
3047
3048 return err;
3049 }
3050
3051 /**
3052 * impose_model_smpl:
3053 * @pmod: pointer to model.
3054 * @dset: dataset information.
3055 *
3056 * Sets on @dset the sample range (starting and ending
3057 * observations) that was in effect when @pmod was estimated.
3058 * This is not always the same as the data range over which
3059 * @pmod was actually estimated (e.g. in case of
3060 * autoregressive models, where observations are dropped
3061 * to allow for lags).
3062 */
3063
impose_model_smpl(const MODEL * pmod,DATASET * dset)3064 void impose_model_smpl (const MODEL *pmod, DATASET *dset)
3065 {
3066 dset->t1 = pmod->smpl.t1;
3067 dset->t2 = pmod->smpl.t2;
3068 #if MDEBUG
3069 fprintf(stderr, "impose_model_smpl: set t1=%d, t2=%d\n",
3070 dset->t1, dset->t2);
3071 #endif
3072 }
3073
3074 #ifdef HAVE_X12A
3075
maybe_delete_x12_file(const MODEL * pmod)3076 static void maybe_delete_x12_file (const MODEL *pmod)
3077 {
3078 char *fname = NULL;
3079
3080 fname = gretl_model_get_data(pmod, "x12a_output");
3081 if (fname != NULL) {
3082 gretl_remove(fname);
3083 }
3084 }
3085
3086 #endif
3087
destroy_all_data_items(MODEL * pmod)3088 static void destroy_all_data_items (MODEL *pmod)
3089 {
3090 int i;
3091
3092 if (pmod->n_data_items == 0) {
3093 return;
3094 }
3095
3096 #ifdef HAVE_X12A
3097 maybe_delete_x12_file(pmod);
3098 #endif
3099
3100 for (i=0; i<pmod->n_data_items; i++) {
3101 free_model_data_item(pmod->data_items[i]);
3102 }
3103
3104 free(pmod->data_items);
3105 pmod->data_items = NULL;
3106 }
3107
discard_model_data_item(MODEL * pmod,const char * key,int free_data)3108 static int discard_model_data_item (MODEL *pmod, const char *key,
3109 int free_data)
3110 {
3111 model_data_item *junk = NULL;
3112 int i, targ = 0;
3113 int err = 0;
3114
3115 for (i=0; i<pmod->n_data_items; i++) {
3116 if (!strcmp(key, pmod->data_items[i]->key)) {
3117 junk = pmod->data_items[i];
3118 targ = i;
3119 break;
3120 }
3121 }
3122
3123 if (junk == NULL) {
3124 err = 1;
3125 } else {
3126 int n_items = pmod->n_data_items - 1;
3127
3128 if (n_items == 0) {
3129 free(pmod->data_items);
3130 pmod->data_items = NULL;
3131 } else {
3132 model_data_item **items;
3133
3134 for (i=targ; i<n_items; i++) {
3135 pmod->data_items[i] = pmod->data_items[i+1];
3136 }
3137 items = realloc(pmod->data_items, n_items * sizeof *items);
3138 if (items != NULL) {
3139 pmod->data_items = items;
3140 }
3141 }
3142
3143 pmod->n_data_items -= 1;
3144
3145 if (free_data) {
3146 /* deep free the data item */
3147 free_model_data_item(junk);
3148 } else {
3149 /* just free the item, not the actual data */
3150 free(junk->key);
3151 free(junk);
3152 }
3153 }
3154
3155 return err;
3156 }
3157
3158 /**
3159 * gretl_model_destroy_data_item:
3160 * @pmod: pointer to model.
3161 * @key: key string.
3162 *
3163 * Looks up the data pointer, attached to @pmod, that is
3164 * identified by @key, and if a pointer is found, frees
3165 * it (or applies the destructor function that was set for
3166 * the item, if any) and removes it from the model's list of
3167 * data items. If you want to remove the item from the
3168 * model's list without freeing the underlying data pointer,
3169 * use gretl_model_detach_data_item().
3170 *
3171 * Returns: 0 on success, 1 on failure (pointer not found).
3172 */
3173
gretl_model_destroy_data_item(MODEL * pmod,const char * key)3174 int gretl_model_destroy_data_item (MODEL *pmod, const char *key)
3175 {
3176 return discard_model_data_item(pmod, key, 1);
3177 }
3178
3179 /**
3180 * gretl_model_detach_data_item:
3181 * @pmod: pointer to model.
3182 * @key: key string.
3183 *
3184 * Looks up the data item, attached to @pmod, that is
3185 * identified by @key, and if an item is found, removes
3186 * it from the model's list of such items. The data
3187 * pointer associated with @key is not touched. If you
3188 * want the underlying resources associated with @key to be
3189 * freed, use gretl_model_destroy_data_item().
3190 *
3191 * Returns: 0 on success, 1 on failure (key not found).
3192 */
3193
gretl_model_detach_data_item(MODEL * pmod,const char * key)3194 int gretl_model_detach_data_item (MODEL *pmod, const char *key)
3195 {
3196 return discard_model_data_item(pmod, key, 0);
3197 }
3198
3199 /**
3200 * gretl_model_set_auxiliary:
3201 * @pmod: pointer to model.
3202 * @aux: code indicating a model's function in an
3203 * auxiliary role (typically, in relation to a hypothesis
3204 * test on another model).
3205 *
3206 * Sets an auxiliary code on @pmod, which may be relevant
3207 * for how the model is printed.
3208 */
3209
gretl_model_set_auxiliary(MODEL * pmod,ModelAuxCode aux)3210 void gretl_model_set_auxiliary (MODEL *pmod, ModelAuxCode aux)
3211 {
3212 pmod->aux = aux;
3213 }
3214
3215 /**
3216 * gretl_model_add_arinfo:
3217 * @pmod: pointer to model.
3218 * @nterms: number of autoregressive coefficients.
3219 *
3220 * Performs initial setup for structure to hold info
3221 * on autoregressive coefficients.
3222 *
3223 * Returns: 0 on success, 1 on error.
3224 */
3225
gretl_model_add_arinfo(MODEL * pmod,int nterms)3226 int gretl_model_add_arinfo (MODEL *pmod, int nterms)
3227 {
3228 int i;
3229
3230 pmod->arinfo = malloc(sizeof *pmod->arinfo);
3231 if (pmod->arinfo == NULL) {
3232 return 1;
3233 }
3234
3235 if (nterms == 0) {
3236 pmod->arinfo->arlist = NULL;
3237 pmod->arinfo->rho = NULL;
3238 pmod->arinfo->sderr = NULL;
3239 return 0;
3240 }
3241
3242 pmod->arinfo->arlist = gretl_list_new(nterms);
3243 if (pmod->arinfo->arlist == NULL) {
3244 free(pmod->arinfo);
3245 pmod->arinfo = NULL;
3246 return 1;
3247 }
3248
3249 pmod->arinfo->rho = malloc(nterms * sizeof *pmod->arinfo->rho);
3250 if (pmod->arinfo->rho == NULL) {
3251 free(pmod->arinfo->arlist);
3252 free(pmod->arinfo);
3253 pmod->arinfo = NULL;
3254 return 1;
3255 }
3256
3257 pmod->arinfo->sderr = malloc(nterms * sizeof *pmod->arinfo->sderr);
3258 if (pmod->arinfo->sderr == NULL) {
3259 free(pmod->arinfo->arlist);
3260 free(pmod->arinfo->rho);
3261 free(pmod->arinfo);
3262 pmod->arinfo = NULL;
3263 return 1;
3264 }
3265
3266 for (i=0; i<nterms; i++) {
3267 pmod->arinfo->sderr[i] = pmod->arinfo->rho[i] = NADBL;
3268 }
3269
3270 return 0;
3271 }
3272
model_stats_init(MODEL * pmod)3273 static void model_stats_init (MODEL *pmod)
3274 {
3275 int i;
3276
3277 pmod->ess = pmod->tss = NADBL;
3278 pmod->sigma = NADBL;
3279 pmod->rsq = pmod->adjrsq = NADBL;
3280 pmod->fstt = pmod->chisq = NADBL;
3281 pmod->lnL = NADBL;
3282 pmod->ybar = pmod->sdy = NADBL;
3283 pmod->dw = pmod->rho = NADBL;
3284
3285 for (i=0; i<C_MAX; i++) {
3286 pmod->criterion[i] = NADBL;
3287 }
3288 }
3289
gretl_model_init_pointers(MODEL * pmod)3290 static void gretl_model_init_pointers (MODEL *pmod)
3291 {
3292 pmod->list = NULL;
3293 pmod->submask = NULL;
3294 pmod->missmask = NULL;
3295 pmod->coeff = NULL;
3296 pmod->sderr = NULL;
3297 pmod->yhat = NULL;
3298 pmod->uhat = NULL;
3299 pmod->xpx = NULL;
3300 pmod->vcv = NULL;
3301 pmod->arinfo = NULL;
3302 pmod->name = NULL;
3303 pmod->depvar = NULL;
3304 pmod->params = NULL;
3305 pmod->tests = NULL;
3306 pmod->dataset = NULL;
3307 pmod->data_items = NULL;
3308 }
3309
3310 /**
3311 * gretl_model_init:
3312 * @pmod: pointer to model.
3313 * @dset: pointer to dataset.
3314 *
3315 * Initializes a gretl #MODEL, including setting its pointer members
3316 * to %NULL. This initialization should be done if the caller has
3317 * declared a #MODEL struct directly, rather than obtaining a pointer to
3318 * #MODEL using gretl_model_new() (in which case the initialization is
3319 * done automatically).
3320 */
3321
gretl_model_init(MODEL * pmod,const DATASET * dset)3322 void gretl_model_init (MODEL *pmod, const DATASET *dset)
3323 {
3324 if (pmod == NULL) return;
3325
3326 #if MDEBUG
3327 fprintf(stderr, "gretl_model_init: pmod at %p\n", (void *) pmod);
3328 #endif
3329
3330 pmod->ID = 0;
3331 pmod->refcount = 0;
3332 pmod->ci = 0;
3333 pmod->opt = OPT_NONE;
3334 pmod->full_n = 0;
3335 pmod->t1 = pmod->t2 = 0;
3336 pmod->nobs = 0;
3337
3338 if (dset != NULL) {
3339 pmod->smpl.t1 = dset->t1;
3340 pmod->smpl.t2 = dset->t2;
3341 pmod->smpl.rseed = dset->rseed;
3342 } else {
3343 pmod->smpl.t1 = 0;
3344 pmod->smpl.t2 = 0;
3345 pmod->smpl.rseed = 0;
3346 }
3347
3348 pmod->ncoeff = 0;
3349 pmod->ntests = 0;
3350 pmod->nparams = 0;
3351 pmod->errcode = 0;
3352 pmod->ifc = 0;
3353 pmod->nwt = 0;
3354 pmod->aux = AUX_NONE;
3355 pmod->esttime = 0;
3356
3357 model_stats_init(pmod);
3358
3359 gretl_model_init_pointers(pmod);
3360 pmod->n_data_items = 0;
3361 }
3362
3363 /**
3364 * gretl_model_smpl_init:
3365 * @pmod: pointer to model.
3366 * @dset: dataset information.
3367 *
3368 * Records the start and end of the current sample range in
3369 * the model @pmod, which may be necessary for future reference
3370 * if a hypothesis test is to be performed. Note that this
3371 * sample range may not be the same as the data range over
3372 * which the model is actually estimated (for example, in the
3373 * case of autoregressive models where observations have to
3374 * be dropped to allow for lags).
3375 */
3376
gretl_model_smpl_init(MODEL * pmod,const DATASET * dset)3377 void gretl_model_smpl_init (MODEL *pmod, const DATASET *dset)
3378 {
3379 pmod->smpl.t1 = dset->t1;
3380 pmod->smpl.t2 = dset->t2;
3381 pmod->smpl.rseed = dset->rseed;
3382 #if MDEBUG
3383 fprintf(stderr, "gretl_model_smpl_init: set t1=%d, t2=%d\n",
3384 dset->t1, dset->t2);
3385 #endif
3386 }
3387
3388 /**
3389 * gretl_model_new:
3390 *
3391 * Allocates memory for a gretl #MODEL struct and initializes the struct,
3392 * using gretl_model_init().
3393 *
3394 * Returns: pointer to model (or %NULL if allocation fails).
3395 */
3396
gretl_model_new(void)3397 MODEL *gretl_model_new (void)
3398 {
3399 MODEL *pmod = malloc(sizeof *pmod);
3400
3401 #if MDEBUG
3402 fprintf(stderr, "gretl_model_new: model at %p\n", (void *) pmod);
3403 if (pmod != NULL) fprintf(stderr, " calling gretl_model_init\n");
3404 #endif
3405
3406 if (pmod != NULL) {
3407 gretl_model_init(pmod, NULL);
3408 }
3409
3410 return pmod;
3411 }
3412
3413 /**
3414 * gretl_model_array_new:
3415 * @n: number of models in array.
3416 *
3417 * Allocates memory for an array of @n gretl #MODEL structs and
3418 * initializes each model using gretl_model_init().
3419 *
3420 * Returns: pointer to models array (or %NULL if allocation fails).
3421 */
3422
gretl_model_array_new(int n)3423 MODEL **gretl_model_array_new (int n)
3424 {
3425 MODEL **models;
3426 int i, j;
3427
3428 models = malloc(n * sizeof *models);
3429
3430 if (models != NULL) {
3431 for (i=0; i<n; i++) {
3432 models[i] = gretl_model_new();
3433 if (models[i] == NULL) {
3434 for (j=0; j<i; j++) {
3435 free(models[i]);
3436 }
3437 free(models);
3438 models = NULL;
3439 break;
3440 }
3441 }
3442 }
3443
3444 return models;
3445 }
3446
3447 /**
3448 * allocate_working_model:
3449 *
3450 * Allocates memory for gretl #MODEL struct and initializes it.
3451 * The model is "protected" against deletion.
3452 *
3453 * Returns: pointer to model (or %NULL if allocation fails).
3454 */
3455
allocate_working_model(void)3456 MODEL *allocate_working_model (void)
3457 {
3458 MODEL *model = gretl_model_new();
3459
3460 if (model != NULL) {
3461 int err = gretl_model_protect(model);
3462
3463 if (err) {
3464 gretl_model_free(model);
3465 model = NULL;
3466 }
3467 }
3468
3469 return model;
3470 }
3471
3472 /**
3473 * gretl_model_array_destroy:
3474 * @models: array of gretl models.
3475 * @n: number of models in array.
3476 *
3477 * Frees all resources associated with an array of models, which
3478 * should have been obtained via gretl_model_array_new().
3479 */
3480
gretl_model_array_destroy(MODEL ** models,int n)3481 void gretl_model_array_destroy (MODEL **models, int n)
3482 {
3483 int i;
3484
3485 if (models != NULL) {
3486 for (i=0; i<n; i++) {
3487 clear_model(models[i]);
3488 free(models[i]);
3489 }
3490 free(models);
3491 }
3492 }
3493
destroy_working_model(MODEL * model)3494 void destroy_working_model (MODEL *model)
3495 {
3496 if (model != NULL) {
3497 gretl_model_free_on_exit(model);
3498 }
3499 }
3500
clear_ar_info(MODEL * pmod)3501 static void clear_ar_info (MODEL *pmod)
3502 {
3503 if (pmod->arinfo->arlist) {
3504 free(pmod->arinfo->arlist);
3505 }
3506 if (pmod->arinfo->rho) {
3507 free(pmod->arinfo->rho);
3508 }
3509 if (pmod->arinfo->sderr) {
3510 free(pmod->arinfo->sderr);
3511 }
3512
3513 free(pmod->arinfo);
3514 pmod->arinfo = NULL;
3515 }
3516
3517 #if MDEBUG > 1
3518
3519 static void
debug_print_model_info(const MODEL * pmod,const char * msg)3520 debug_print_model_info (const MODEL *pmod, const char *msg)
3521 {
3522 fprintf(stderr, "%s:\n"
3523 " pmod = %p\n"
3524 " pmod->list = %p\n"
3525 " pmod->submask = %p\n"
3526 " pmod->missmask = %p\n"
3527 " pmod->coeff = %p\n"
3528 " pmod->sderr = %p\n"
3529 " pmod->yhat = %p\n"
3530 " pmod->uhat = %p\n"
3531 " pmod->xpx = %p\n"
3532 " pmod->vcv = %p\n"
3533 " pmod->name = %p\n"
3534 " pmod->depvar = %p\n"
3535 " pmod->params = %p\n"
3536 " pmod->arinfo = %p\n"
3537 " pmod->tests = %p\n"
3538 " pmod->data_items = %p\n", msg,
3539 (void *) pmod, (void *) pmod->list,
3540 (void *) pmod->submask, (void *) pmod->missmask,
3541 (void *) pmod->coeff, (void *) pmod->sderr,
3542 (void *) pmod->yhat, (void *) pmod->uhat,
3543 (void *) pmod->xpx,
3544 (void *) pmod->vcv, (void *) pmod->name,
3545 (void *) pmod->depvar, (void *) pmod->params,
3546 (void *) pmod->arinfo, (void *) pmod->tests,
3547 (void *) pmod->data_items);
3548 }
3549
3550 #endif /* MDEBUG */
3551
3552 /**
3553 * gretl_model_destroy_tests:
3554 * @pmod: pointer to model.
3555 *
3556 * Clears any hypothesis test structs that have been attached
3557 * to @pmod.
3558 */
3559
gretl_model_destroy_tests(MODEL * pmod)3560 void gretl_model_destroy_tests (MODEL *pmod)
3561 {
3562 if (pmod != NULL && pmod->ntests > 0) {
3563 int i;
3564
3565 for (i=0; i<pmod->ntests; i++) {
3566 if (pmod->tests[i].param != NULL) {
3567 free(pmod->tests[i].param);
3568 }
3569 }
3570 free(pmod->tests);
3571 pmod->tests = NULL;
3572 pmod->ntests = 0;
3573 }
3574 }
3575
3576 /**
3577 * clear_model:
3578 * @pmod: pointer to model.
3579 *
3580 * Clears a gretl #MODEL, freeing all allocated storage and setting
3581 * pointer members to %NULL. Also frees any data pointers attached
3582 * via gretl_model_set_data(). The model pointer itself is not
3583 * freed, so this function may be called on a #MODEL which has been
3584 * declared directly by the caller; in that case the caller should
3585 * pass the address of the #MODEL in question.
3586 */
3587
clear_model(MODEL * pmod)3588 void clear_model (MODEL *pmod)
3589 {
3590 if (pmod != NULL) {
3591 #if MDEBUG
3592 fprintf(stderr, "clear model: model at %p\n", (void *) pmod);
3593 #endif
3594
3595 #if MDEBUG > 1
3596 debug_print_model_info(pmod, "Doing clear_model");
3597 #endif
3598 if (pmod->list != NULL) free(pmod->list);
3599 if (pmod->missmask != NULL) free(pmod->missmask);
3600 if (pmod->coeff != NULL) free(pmod->coeff);
3601 if (pmod->sderr != NULL) free(pmod->sderr);
3602 if (pmod->yhat != NULL) free(pmod->yhat);
3603 if (pmod->uhat != NULL) free(pmod->uhat);
3604 if (pmod->xpx != NULL) free(pmod->xpx);
3605 if (pmod->vcv != NULL) free(pmod->vcv);
3606 if (pmod->name != NULL) free(pmod->name);
3607 if (pmod->depvar != NULL) free(pmod->depvar);
3608
3609 if (pmod->submask != NULL) {
3610 free_subsample_mask(pmod->submask);
3611 }
3612
3613 if (pmod->arinfo != NULL) {
3614 clear_ar_info(pmod);
3615 }
3616 if (pmod->params != NULL) {
3617 strings_array_free(pmod->params, pmod->nparams);
3618 }
3619
3620 destroy_dataset(pmod->dataset);
3621 gretl_model_destroy_tests(pmod);
3622 destroy_all_data_items(pmod);
3623 }
3624
3625 /* this may be redundant */
3626 #if MDEBUG
3627 fprintf(stderr, "clear_model, calling gretl_model_init\n");
3628 #endif
3629 #if 1
3630 gretl_model_init(pmod, NULL);
3631 #endif
3632 }
3633
clear_model_xpx(MODEL * pmod)3634 void clear_model_xpx (MODEL *pmod)
3635 {
3636 if (pmod->xpx != NULL) {
3637 free(pmod->xpx);
3638 pmod->xpx = NULL;
3639 }
3640 }
3641
3642 /**
3643 * gretl_model_free:
3644 * @pmod: pointer to #MODEL.
3645 *
3646 * Free allocated content of @pmod then the pointer itself.
3647 */
3648
gretl_model_free(MODEL * pmod)3649 void gretl_model_free (MODEL *pmod)
3650 {
3651 if (pmod != NULL) {
3652 clear_model(pmod);
3653 #if MDEBUG
3654 fprintf(stderr, "gretl_model_free: freeing at %p\n", (void *) pmod);
3655 #endif
3656 free(pmod);
3657 }
3658 }
3659
3660 /**
3661 * gretl_model_free_on_exit:
3662 * @pmod: pointer to #MODEL.
3663 *
3664 * Free allocated content of @pmod then the pointer itself,
3665 * without regard to the model's reference count.
3666 */
3667
gretl_model_free_on_exit(MODEL * pmod)3668 void gretl_model_free_on_exit (MODEL *pmod)
3669 {
3670 #if MDEBUG
3671 fprintf(stderr, "gretl_model_free_on_exit: pmod at %p\n", (void *) pmod);
3672 #endif
3673 if (pmod != NULL) {
3674 remove_model_from_stack_on_exit(pmod);
3675 clear_model(pmod);
3676 free(pmod);
3677 }
3678 }
3679
copy_test(ModelTest * targ,const ModelTest * src)3680 static void copy_test (ModelTest *targ, const ModelTest *src)
3681 {
3682 targ->type = src->type;
3683
3684 if (src->param != NULL && *src->param != '\0') {
3685 targ->param = gretl_strdup(src->param);
3686 } else {
3687 targ->param = NULL;
3688 }
3689
3690 targ->teststat = src->teststat;
3691 targ->dfn = src->dfn;
3692 targ->dfd = src->dfd;
3693 targ->order = src->order;
3694 targ->value = src->value;
3695 targ->pvalue = src->pvalue;
3696 targ->crit = src->crit;
3697 targ->alpha = src->alpha;
3698 targ->opt = src->opt;
3699 }
3700
copy_model_tests(MODEL * targ,const MODEL * src)3701 static int copy_model_tests (MODEL *targ, const MODEL *src)
3702 {
3703 int i, n = src->ntests;
3704
3705 if (n <= 0 || src->tests == NULL) {
3706 return 0;
3707 }
3708
3709 targ->tests = malloc(n * sizeof *targ->tests);
3710 if (targ->tests == NULL) return 1;
3711
3712 for (i=0; i<n; i++) {
3713 copy_test(&targ->tests[i], &src->tests[i]);
3714 }
3715
3716 return 0;
3717 }
3718
real_add_test_to_model(MODEL * pmod,ModelTest * test)3719 static int real_add_test_to_model (MODEL *pmod, ModelTest *test)
3720 {
3721 ModelTest *tests = NULL;
3722 int nt = pmod->ntests;
3723 int err = 0;
3724
3725 tests = realloc(pmod->tests, (nt + 1) * sizeof *tests);
3726
3727 if (tests == NULL) {
3728 err = E_ALLOC;
3729 } else {
3730 pmod->tests = tests;
3731 pmod->ntests += 1;
3732 copy_test(&pmod->tests[nt], test);
3733 }
3734
3735 return err;
3736 }
3737
attach_model_params_from_xml(xmlNodePtr node,xmlDocPtr doc,MODEL * pmod)3738 static int attach_model_params_from_xml (xmlNodePtr node, xmlDocPtr doc,
3739 MODEL *pmod)
3740 {
3741 char **S;
3742 int slop = (pmod->ci == GARCH);
3743 int np = 0;
3744 int err = 0;
3745
3746 S = gretl_xml_get_strings_array(node, doc, &np, slop, &err);
3747
3748 if (!err) {
3749 pmod->params = S;
3750 pmod->nparams = np;
3751 }
3752
3753 return err;
3754 }
3755
attach_model_tests_from_xml(MODEL * pmod,xmlNodePtr node)3756 int attach_model_tests_from_xml (MODEL *pmod, xmlNodePtr node)
3757 {
3758 ModelTest test;
3759 xmlNodePtr cur = node->xmlChildrenNode;
3760 int got, err = 0;
3761
3762 while (cur != NULL && !err) {
3763 gretl_test_init(&test, 0);
3764 got = 0;
3765 got += gretl_xml_get_prop_as_int(cur, "type", &test.type);
3766 got += gretl_xml_get_prop_as_uchar(cur, "teststat", &test.teststat);
3767 got += gretl_xml_get_prop_as_int(cur, "dfn", &test.dfn);
3768 got += gretl_xml_get_prop_as_double(cur, "dfd", &test.dfd);
3769 got += gretl_xml_get_prop_as_int(cur, "order", &test.order);
3770 got += gretl_xml_get_prop_as_double(cur, "value", &test.value);
3771 got += gretl_xml_get_prop_as_double(cur, "pvalue", &test.pvalue);
3772 got += gretl_xml_get_prop_as_double(cur, "crit", &test.crit);
3773 got += gretl_xml_get_prop_as_double(cur, "alpha", &test.alpha);
3774 if (got < 7) {
3775 err = E_DATA;
3776 } else {
3777 gretl_xml_get_prop_as_opt(cur, "opt", &test.opt);
3778 gretl_xml_get_prop_as_string(cur, "param", &test.param);
3779 err = real_add_test_to_model(pmod, &test);
3780 free(test.param); /* copied by real_add_test_to_model */
3781 }
3782 cur = cur->next;
3783 }
3784
3785 return err;
3786 }
3787
serialize_test(const ModelTest * src,PRN * prn)3788 static void serialize_test (const ModelTest *src, PRN *prn)
3789 {
3790 pprintf(prn, "<test type=\"%d\" ", src->type);
3791
3792 if (src->param != NULL && *src->param != '\0') {
3793 pprintf(prn, "param=\"%s\" ", src->param);
3794 }
3795
3796 pprintf(prn, "teststat=\"%d\" ", (int) src->teststat);
3797 pprintf(prn, "dfn=\"%d\" ", src->dfn);
3798 pprintf(prn, "dfd=\"%g\" ", src->dfd);
3799 pprintf(prn, "order=\"%d\" ", src->order);
3800 pprintf(prn, "value=\"%.15g\" ", src->value);
3801 pprintf(prn, "pvalue=\"%.15g\" ", src->pvalue);
3802
3803 if (!na(src->crit)) {
3804 pprintf(prn, "crit=\"%g\" ", src->crit);
3805 pprintf(prn, "alpha=\"%g\" ", src->alpha);
3806 }
3807
3808 if (src->opt != OPT_NONE) {
3809 pprintf(prn, "opt=\"%u\" ", (unsigned) src->opt);
3810 }
3811
3812 pputs(prn, "/>\n");
3813 }
3814
bundlize_test(const ModelTest * src)3815 static gretl_bundle *bundlize_test (const ModelTest *src)
3816 {
3817 gretl_bundle *b = gretl_bundle_new();
3818
3819 if (b == NULL) {
3820 return NULL;
3821 }
3822
3823 if (src->param != NULL && *src->param != '\0') {
3824 gretl_bundle_set_string(b, "param", src->param);
3825 }
3826
3827 if (src->dfn > 0) {
3828 gretl_bundle_set_scalar(b, "dfn", (double) src->dfn);
3829 }
3830 if (src->dfd > 0) {
3831 gretl_bundle_set_scalar(b, "dfd", src->dfd);
3832 }
3833 if (src->order > 0) {
3834 gretl_bundle_set_scalar(b, "order", src->order);
3835 }
3836 if (!na(src->value)) {
3837 gretl_bundle_set_scalar(b, "test", src->value);
3838 }
3839 if (!na(src->pvalue)) {
3840 gretl_bundle_set_scalar(b, "pvalue", src->pvalue);
3841 }
3842 if (!na(src->crit)) {
3843 gretl_bundle_set_scalar(b, "crit", src->crit);
3844 gretl_bundle_set_scalar(b, "alpha", src->alpha);
3845 }
3846
3847 return b;
3848 }
3849
serialize_model_tests(const MODEL * pmod,PRN * prn)3850 static int serialize_model_tests (const MODEL *pmod, PRN *prn)
3851 {
3852 int i, n = pmod->ntests;
3853
3854 if (n <= 0 || pmod->tests == NULL) {
3855 return 0;
3856 }
3857
3858 pprintf(prn, "<tests count=\"%d\">\n", pmod->ntests);
3859
3860 for (i=0; i<n; i++) {
3861 serialize_test(&pmod->tests[i], prn);
3862 }
3863
3864 pputs(prn, "</tests>\n");
3865
3866 return 0;
3867 }
3868
testvals_differ(double x,double y)3869 static int testvals_differ (double x, double y)
3870 {
3871 double eq_tol = 1.0e-10;
3872 double reldiff;
3873
3874 if (x == 0.0) {
3875 reldiff = fabs(y);
3876 } else if (y == 0.0) {
3877 reldiff = fabs(x);
3878 } else if (x > y) {
3879 reldiff = fabs((x - y) / y);
3880 } else {
3881 reldiff = fabs((y - x) / x);
3882 }
3883
3884 return reldiff > eq_tol;
3885 }
3886
test_params_differ(const char * p,const char * s)3887 static int test_params_differ (const char *p, const char *s)
3888 {
3889 int ret = 0;
3890
3891 if (s == NULL && p != NULL) {
3892 ret = 1;
3893 } else if (p == NULL && s != NULL) {
3894 ret = 1;
3895 } else if (p != NULL && s != NULL) {
3896 ret = strcmp(p, s);
3897 }
3898
3899 return ret;
3900 }
3901
model_tests_differ(ModelTest * mt1,ModelTest * mt2)3902 static int model_tests_differ (ModelTest *mt1, ModelTest *mt2)
3903 {
3904 int ret = 0;
3905
3906 if (mt1->type != mt2->type) {
3907 ret = 1;
3908 } else if (mt1->order != mt2->order) {
3909 ret = 1;
3910 } else if (mt1->teststat != mt2->teststat) {
3911 ret = 1;
3912 } else if (test_params_differ(mt1->param, mt2->param)) {
3913 ret = 1;
3914 } else if (testvals_differ(mt1->value, mt2->value)) {
3915 ret = 1;
3916 }
3917
3918 return ret;
3919 }
3920
3921 /**
3922 * model_test_free:
3923 * @test: object to free.
3924 */
3925
model_test_free(ModelTest * test)3926 void model_test_free (ModelTest *test)
3927 {
3928 free(test->param);
3929 free(test);
3930 }
3931
3932 /**
3933 * model_test_new:
3934 * @ttype: type of test to add.
3935 *
3936 * Returns: new #ModelTest pointer, or %NULL on failure.
3937 */
3938
model_test_new(ModelTestType ttype)3939 ModelTest *model_test_new (ModelTestType ttype)
3940 {
3941 ModelTest *test = malloc(sizeof *test);
3942
3943 if (test != NULL) {
3944 gretl_test_init(test, ttype);
3945 }
3946
3947 return test;
3948 }
3949
3950 /**
3951 * maybe_add_test_to_model:
3952 * @pmod: pointer to model.
3953 * @test: model test to be added.
3954 *
3955 * Adds a #ModelTest to @pmod, if the test in question has
3956 * not already been performed and recorded. Note that this
3957 * function takes care of freeing @test.
3958 *
3959 * Returns: 1 if the test was added, otherwise 0.
3960 */
3961
maybe_add_test_to_model(MODEL * pmod,ModelTest * test)3962 int maybe_add_test_to_model (MODEL *pmod, ModelTest *test)
3963 {
3964 int i, done = 0, add = 0;
3965
3966 if (test == NULL || test->teststat == GRETL_STAT_NONE) {
3967 return 0;
3968 }
3969
3970 for (i=0; i<pmod->ntests; i++) {
3971 if (!model_tests_differ(test, &pmod->tests[i])) {
3972 done = 1;
3973 }
3974 }
3975
3976 if (!done) {
3977 int n = pmod->ntests + 1;
3978 ModelTest *tests;
3979
3980 tests = realloc(pmod->tests, n * sizeof *tests);
3981
3982 if (tests != NULL) {
3983 pmod->tests = tests;
3984 copy_test(&pmod->tests[n-1], test);
3985 pmod->ntests += 1;
3986 add = 1;
3987 }
3988 }
3989
3990 free(test->param);
3991 free(test);
3992
3993 return add;
3994 }
3995
gretl_model_get_normality_test(const MODEL * pmod,PRN * prn)3996 int gretl_model_get_normality_test (const MODEL *pmod, PRN *prn)
3997 {
3998 ModelTest *test = NULL;
3999 int i, err = 0;
4000
4001 for (i=0; i<pmod->ntests; i++) {
4002 if (pmod->tests[i].type == GRETL_TEST_NORMAL) {
4003 test = &pmod->tests[i];
4004 break;
4005 }
4006 }
4007
4008 if (test == NULL) {
4009 err = E_BADSTAT;
4010 } else {
4011 record_test_result(test->value, test->pvalue);
4012 gretl_model_test_print(pmod, i, prn);
4013 }
4014
4015 return err;
4016 }
4017
model_test_set_teststat(ModelTest * test,unsigned char ts)4018 void model_test_set_teststat (ModelTest *test, unsigned char ts)
4019 {
4020 test->teststat = ts;
4021 }
4022
model_test_set_order(ModelTest * test,int order)4023 void model_test_set_order (ModelTest *test, int order)
4024 {
4025 test->order = order;
4026 }
4027
model_test_set_dfn(ModelTest * test,int df)4028 void model_test_set_dfn (ModelTest *test, int df)
4029 {
4030 test->dfn = df;
4031 }
4032
model_test_set_dfd(ModelTest * test,double df)4033 void model_test_set_dfd (ModelTest *test, double df)
4034 {
4035 test->dfd = df;
4036 }
4037
model_test_set_value(ModelTest * test,double val)4038 void model_test_set_value (ModelTest *test, double val)
4039 {
4040 test->value = val;
4041 }
4042
model_test_set_pvalue(ModelTest * test,double pval)4043 void model_test_set_pvalue (ModelTest *test, double pval)
4044 {
4045 test->pvalue = pval;
4046 }
4047
model_test_set_opt(ModelTest * test,gretlopt opt)4048 void model_test_set_opt (ModelTest *test, gretlopt opt)
4049 {
4050 test->opt = opt;
4051 }
4052
model_test_set_crit_and_alpha(ModelTest * test,double crit,double alpha)4053 void model_test_set_crit_and_alpha (ModelTest *test,
4054 double crit,
4055 double alpha)
4056 {
4057 test->crit = crit;
4058 test->alpha = alpha;
4059 }
4060
model_test_set_param(ModelTest * test,const char * s)4061 void model_test_set_param (ModelTest *test, const char *s)
4062 {
4063 test->param = gretl_strdup(s);
4064 }
4065
model_test_set_allocated_param(ModelTest * test,char * s)4066 void model_test_set_allocated_param (ModelTest *test, char *s)
4067 {
4068 test->param = s;
4069 }
4070
4071 struct test_strings {
4072 ModelTestType ID;
4073 const char *descrip;
4074 const char *H0;
4075 };
4076
4077 static struct test_strings tstrings[] = {
4078 { GRETL_TEST_ADD,
4079 N_("Test for addition of variables"),
4080 N_("parameters are zero for the variables") },
4081 { GRETL_TEST_ARCH,
4082 N_("Test for ARCH of order %s"),
4083 N_("no ARCH effect is present") },
4084 { GRETL_TEST_AUTOCORR,
4085 N_("LM test for autocorrelation up to order %s"),
4086 N_("no autocorrelation") },
4087 { GRETL_TEST_CHOW,
4088 N_("Chow test for structural break at observation %s"),
4089 N_("no structural break") },
4090 { GRETL_TEST_CHOWDUM,
4091 N_("Chow test for structural difference with respect to %s"),
4092 N_("no structural difference") },
4093 { GRETL_TEST_CUSUM,
4094 N_("CUSUM test for parameter stability"),
4095 N_("no change in parameters") },
4096 { GRETL_TEST_QLR,
4097 N_("QLR test for structural break"),
4098 N_("no structural break") },
4099 { GRETL_TEST_GROUPWISE,
4100 N_("Likelihood ratio test for groupwise heteroskedasticity"),
4101 N_("the units have a common error variance") },
4102 { GRETL_TEST_LOGS,
4103 N_("Non-linearity test (logs)"),
4104 N_("relationship is linear") },
4105 { GRETL_TEST_NORMAL,
4106 N_("Test for normality of residual"),
4107 N_("error is normally distributed") },
4108 { GRETL_TEST_OMIT,
4109 N_("Test for omission of variables"),
4110 N_("parameters are zero for the variables") },
4111 { GRETL_TEST_RESET,
4112 N_("RESET test for specification"),
4113 N_("specification is adequate") },
4114 { GRETL_TEST_SQUARES,
4115 N_("Non-linearity test (squares)"),
4116 N_("relationship is linear") },
4117 { GRETL_TEST_WHITES,
4118 N_("White's test for heteroskedasticity"),
4119 N_("heteroskedasticity not present") },
4120 { GRETL_TEST_BP,
4121 N_("Breusch-Pagan test for heteroskedasticity"),
4122 N_("heteroskedasticity not present") },
4123 { GRETL_TEST_SARGAN,
4124 N_("Sargan over-identification test"),
4125 N_("all instruments are valid") },
4126 { GRETL_TEST_IV_HAUSMAN,
4127 N_("Hausman test"),
4128 N_("OLS estimates are consistent") },
4129 { GRETL_TEST_PANEL_HAUSMAN,
4130 N_("Hausman test"),
4131 N_("GLS estimates are consistent") },
4132 { GRETL_TEST_PANEL_F,
4133 N_("Test for differing group intercepts"),
4134 N_("The groups have a common intercept") },
4135 { GRETL_TEST_PANEL_BP,
4136 N_("Breusch-Pagan test"),
4137 N_("Variance of the unit-specific error = 0") },
4138 { GRETL_TEST_PANEL_TIMEDUM,
4139 N_("Wald joint test on time dummies"),
4140 N_("No time effects") },
4141 { GRETL_TEST_HET_1,
4142 N_("Pesaran-Taylor test for heteroskedasticity"),
4143 N_("heteroskedasticity not present") },
4144 { GRETL_TEST_COMFAC,
4145 N_("Test of common factor restriction"),
4146 N_("restriction is acceptable") },
4147 { GRETL_TEST_INDEP,
4148 N_("Test of independence"),
4149 N_("rho = 0") },
4150 { GRETL_TEST_RE,
4151 N_("LR test for rho = 0"),
4152 NULL },
4153 { GRETL_TEST_WITHIN_F,
4154 N_("Joint test on named regressors"),
4155 NULL },
4156 { GRETL_TEST_RE_WALD,
4157 N_("Joint test on named regressors"),
4158 NULL },
4159 { GRETL_TEST_PANEL_WELCH,
4160 N_("Robust test for differing group intercepts"),
4161 N_("The groups have a common intercept") },
4162 { GRETL_TEST_XDEPEND,
4163 N_("Pesaran CD test for cross-sectional dependence"),
4164 N_("No cross-sectional dependence") },
4165 { GRETL_TEST_PANEL_AR,
4166 N_("Wooldridge test for autocorrelation in panel data"),
4167 N_("No first-order autocorrelation (rho = -0.5)") },
4168 { GRETL_TEST_MAX, NULL, NULL }
4169 };
4170
test_get_descrip(const ModelTest * test)4171 static const char *test_get_descrip (const ModelTest *test)
4172 {
4173 static const char *dfree =
4174 N_("Distribution free Wald test for heteroskedasticity");
4175 static const char *lr_indep =
4176 N_("Likelihood ratio test of independence");
4177
4178 if (test->type == GRETL_TEST_GROUPWISE &&
4179 test->teststat == GRETL_STAT_WALD_CHISQ) {
4180 return dfree;
4181 } else if (test->type == GRETL_TEST_INDEP &&
4182 test->teststat == GRETL_STAT_LR) {
4183 return lr_indep;
4184 } else {
4185 int i;
4186
4187 for (i=0; tstrings[i].ID < GRETL_TEST_MAX; i++) {
4188 if (test->type == tstrings[i].ID) {
4189 return tstrings[i].descrip;
4190 }
4191 }
4192 }
4193
4194 return NULL;
4195 }
4196
test_get_opt_string(const ModelTest * test)4197 static const char *test_get_opt_string (const ModelTest *test)
4198 {
4199 if (test->type == GRETL_TEST_RESET) {
4200 if (test->opt & OPT_R) {
4201 return N_("squares only");
4202 } else if (test->opt & OPT_C) {
4203 return N_("cubes only");
4204 }
4205 } else if (test->type == GRETL_TEST_WHITES) {
4206 if (test->opt & OPT_X) {
4207 return N_("squares only");
4208 }
4209 } else if (test->type == GRETL_TEST_BP) {
4210 if (test->opt & OPT_R) {
4211 return N_("robust variant");
4212 }
4213 }
4214
4215 return NULL;
4216 }
4217
print_test_opt(const ModelTest * test,PRN * prn)4218 static void print_test_opt (const ModelTest *test, PRN *prn)
4219 {
4220 const char *optstr = test_get_opt_string(test);
4221
4222 if (optstr != NULL) {
4223 pprintf(prn, " (%s)", _(optstr));
4224 }
4225 }
4226
gretl_test_print_heading(const ModelTest * test,PRN * prn)4227 static int gretl_test_print_heading (const ModelTest *test, PRN *prn)
4228 {
4229 const char *descrip, *param = NULL;
4230 char ordstr[16];
4231
4232 descrip = test_get_descrip(test);
4233 if (descrip == NULL) {
4234 return 1;
4235 }
4236
4237 if (test->order > 0) {
4238 sprintf(ordstr, "%d", test->order);
4239 param = ordstr;
4240 } else if (test->type == GRETL_TEST_CHOW ||
4241 test->type == GRETL_TEST_CHOWDUM) {
4242 param = test->param;
4243 }
4244
4245 if (param != NULL) {
4246 if (plain_format(prn)) {
4247 pprintf(prn, _(descrip), param);
4248 } else {
4249 pprintf(prn, _(descrip), param);
4250 }
4251 } else {
4252 if (plain_format(prn)) {
4253 pputs(prn, _(descrip));
4254 } else {
4255 pputs(prn, _(descrip));
4256 }
4257 }
4258
4259 if (test->opt) {
4260 print_test_opt(test, prn);
4261 }
4262
4263 return 0;
4264 }
4265
print_add_omit_varnames(const char * s,PRN * prn)4266 static int print_add_omit_varnames (const char *s, PRN *prn)
4267 {
4268 const char *endings[] = {
4269 "\n ",
4270 "\\\\\n \\qquad ",
4271 "\\par\n "
4272 };
4273 const char *ending = NULL;
4274
4275 if (s == NULL || *s == '\0') {
4276 return 1;
4277 }
4278
4279 if (plain_format(prn)) {
4280 ending = endings[0];
4281 } else if (tex_format(prn)) {
4282 ending = endings[1];
4283 } else if (rtf_format(prn)) {
4284 ending = endings[2];
4285 } else {
4286 return 1;
4287 }
4288
4289 pputs(prn, ending);
4290
4291 while (*s) {
4292 if (*s == ' ') {
4293 pputs(prn, ending);
4294 } else {
4295 pputc(prn, *s);
4296 }
4297 s++;
4298 }
4299
4300 return 0;
4301 }
4302
gretl_test_print_h_0(const ModelTest * test,int heading,PRN * prn)4303 static void gretl_test_print_h_0 (const ModelTest *test, int heading,
4304 PRN *prn)
4305 {
4306 const char *H0 = NULL;
4307 int i;
4308
4309 if (test->type == GRETL_TEST_PANEL_AR &&
4310 test->teststat == GRETL_STAT_STUDENT) {
4311 H0 = N_("No first-order autocorrelation (rho = 0)");
4312 } else {
4313 for (i=0; tstrings[i].ID < GRETL_TEST_MAX; i++) {
4314 if (test->type == tstrings[i].ID) {
4315 H0 = tstrings[i].H0;
4316 }
4317 }
4318 }
4319
4320 if (test->type == GRETL_TEST_INDEP &&
4321 test->teststat == GRETL_STAT_WALD_CHISQ &&
4322 H0 != NULL && !strcmp(H0, "rho = 0")) {
4323 /* adjust for QML biprobit */
4324 H0 = N_("atanh(rho) = 0");
4325 }
4326
4327 if (H0 == NULL && test->type != GRETL_TEST_WITHIN_F &&
4328 test->type != GRETL_TEST_RE_WALD) {
4329 return;
4330 }
4331
4332 if (heading) {
4333 pputs(prn, " -");
4334 if (tex_format(prn)) {
4335 pputc(prn, '-');
4336 }
4337 }
4338
4339 if (H0 == NULL) {
4340 return;
4341 }
4342
4343 if (plain_format(prn)) {
4344 pprintf(prn, "\n %s: ", _("Null hypothesis"));
4345 pputs(prn, _(H0));
4346 } else if (tex_format(prn)) {
4347 pprintf(prn, "\\\\\n\\quad %s: ", _("Null hypothesis"));
4348 if (!strcmp(H0, "rho = 0")) {
4349 pputs(prn, "$\\rho = 0$");
4350 } else {
4351 pputs(prn, _(H0));
4352 }
4353 } else if (rtf_format(prn)) {
4354 pprintf(prn, "\\par\n %s: ", _("Null hypothesis"));
4355 pputs(prn, _(H0));
4356 }
4357
4358 if (test->type == GRETL_TEST_ADD || test->type == GRETL_TEST_OMIT) {
4359 print_add_omit_varnames(test->param, prn);
4360 }
4361 }
4362
get_test_stat_string(const ModelTest * test,PRN * prn)4363 static gchar *get_test_stat_string (const ModelTest *test, PRN *prn)
4364 {
4365 int tex = tex_format(prn);
4366 gchar *ret = NULL;
4367
4368 switch (test->teststat) {
4369 case GRETL_STAT_LM:
4370 ret = g_strdup_printf("LM = %g", test->value);
4371 break;
4372 case GRETL_STAT_F:
4373 case GRETL_STAT_RESET:
4374 if (tex) {
4375 ret = g_strdup_printf("$F(%d, %g)$ = %g", test->dfn, test->dfd, test->value);
4376 } else {
4377 ret = g_strdup_printf("F(%d, %g) = %g", test->dfn, test->dfd, test->value);
4378 }
4379 break;
4380 case GRETL_STAT_SUP_WALD:
4381 if (tex) {
4382 ret = g_strdup_printf("max $\\chi^2(%d)$ = %g (%s)", test->dfn,
4383 test->value, test->param);
4384 } else {
4385 ret = g_strdup_printf(_("chi-square(%d) = %g at observation %s"),
4386 test->dfn, test->value, test->param);
4387 }
4388 break;
4389 case GRETL_STAT_LMF:
4390 ret = g_strdup_printf("LMF = %g", test->value);
4391 break;
4392 case GRETL_STAT_WF:
4393 if (tex) {
4394 ret = g_strdup_printf("Welch $F(%d, %.1f)$ = %g", test->dfn, test->dfd, test->value);
4395 } else {
4396 ret = g_strdup_printf("Welch F(%d, %.1f) = %g", test->dfn, test->dfd, test->value);
4397 }
4398 break;
4399 case GRETL_STAT_HARVEY_COLLIER:
4400 if (tex) {
4401 ret = g_strdup_printf("Harvey--Collier $t(%d)$ = %g", test->dfn, test->value);
4402 } else {
4403 ret = g_strdup_printf("Harvey-Collier t(%d) = %g", test->dfn, test->value);
4404 }
4405 break;
4406 case GRETL_STAT_NORMAL_CHISQ:
4407 if (tex) {
4408 ret = g_strdup_printf("$\\chi^2(2)$ = %g", test->value);
4409 } else {
4410 ret = g_strdup_printf("%s(2) = %g", _("Chi-square"), test->value);
4411 }
4412 break;
4413 case GRETL_STAT_LR:
4414 case GRETL_STAT_WALD_CHISQ:
4415 case GRETL_STAT_LB_CHISQ:
4416 if (tex) {
4417 if (na(test->value)) {
4418 ret = g_strdup_printf("$\\chi^2(%d)$ = NA (%s)", test->dfn, _("failed"));
4419 } else {
4420 ret = g_strdup_printf("$\\chi^2(%d)$ = %g", test->dfn, test->value);
4421 }
4422 } else {
4423 if (na(test->value)) {
4424 ret = g_strdup_printf("%s(%d) = NA (%s)", _("Chi-square"), test->dfn,
4425 _("failed"));
4426 } else {
4427 ret = g_strdup_printf("%s(%d) = %g", _("Chi-square"), test->dfn, test->value);
4428 }
4429 }
4430 break;
4431 case GRETL_STAT_Z:
4432 if (tex) {
4433 ret = g_strdup_printf("$z$ = %g", test->value);
4434 } else {
4435 ret = g_strdup_printf("z = %g", test->value);
4436 }
4437 break;
4438 case GRETL_STAT_STUDENT:
4439 if (tex) {
4440 ret = g_strdup_printf("$t$(%d) = %g", test->dfn, test->value);
4441 } else {
4442 ret = g_strdup_printf("t(%d) = %g", test->dfn, test->value);
4443 }
4444 break;
4445 default:
4446 break;
4447 }
4448
4449 return ret;
4450 }
4451
get_test_pval_string(const ModelTest * test,PRN * prn)4452 static gchar *get_test_pval_string (const ModelTest *test, PRN *prn)
4453 {
4454 int tex = tex_format(prn);
4455 gchar *ret = NULL;
4456
4457 switch (test->teststat) {
4458 case GRETL_STAT_LM:
4459 if (tex) {
4460 ret = g_strdup_printf("$P$($\\chi^2(%d) >$ %g) = %g",
4461 test->dfn, test->value, test->pvalue);
4462 } else {
4463 ret = g_strdup_printf("P(%s(%d) > %g) = %g", _("Chi-square"),
4464 test->dfn, test->value, test->pvalue);
4465 }
4466 break;
4467 case GRETL_STAT_F:
4468 case GRETL_STAT_RESET:
4469 if (tex) {
4470 ret = g_strdup_printf("$P$($F(%d, %g) >$ %g) = %g",
4471 test->dfn, test->dfd, test->value, test->pvalue);
4472 } else {
4473 ret = g_strdup_printf("P(F(%d, %g) > %g) = %g",
4474 test->dfn, test->dfd, test->value, test->pvalue);
4475 }
4476 break;
4477 case GRETL_STAT_WF:
4478 if (tex) {
4479 ret = g_strdup_printf("$P$($F(%d, %.1f) >$ %g) = %g",
4480 test->dfn, test->dfd, test->value, test->pvalue);
4481 } else {
4482 ret = g_strdup_printf("P(F(%d, %.1f) > %g) = %g",
4483 test->dfn, test->dfd, test->value, test->pvalue);
4484 }
4485 break;
4486 case GRETL_STAT_LMF:
4487 if (tex) {
4488 ret = g_strdup_printf("$P$($F(%d, %g) >$ %g) = %g",
4489 test->dfn, test->dfd, test->value, test->pvalue);
4490 } else {
4491 ret = g_strdup_printf("P(F(%d, %g) > %g) = %g",
4492 test->dfn, test->dfd, test->value, test->pvalue);
4493 }
4494 break;
4495 case GRETL_STAT_HARVEY_COLLIER:
4496 if (tex) {
4497 ret = g_strdup_printf("$P$($t_{%d} >$ %g) = %g",
4498 test->dfn, test->value, test->pvalue);
4499 } else {
4500 ret = g_strdup_printf("P(t(%d) > %g) = %g",
4501 test->dfn, test->value, test->pvalue);
4502 }
4503 break;
4504 case GRETL_STAT_STUDENT:
4505 if (tex) {
4506 ret = g_strdup_printf("$P$($|t| >$ %g) = %g",
4507 fabs(test->value), test->pvalue);
4508 } else {
4509 ret = g_strdup_printf("P(|t| > %g) = %g",
4510 fabs(test->value), test->pvalue);
4511 }
4512 break;
4513 case GRETL_STAT_NORMAL_CHISQ:
4514 case GRETL_STAT_LR:
4515 case GRETL_STAT_WALD_CHISQ:
4516 case GRETL_STAT_SUP_WALD:
4517 case GRETL_STAT_Z:
4518 if (na(test->value)) {
4519 ; /* leave @ret NULL */
4520 } else if (na(test->pvalue)) {
4521 ret = g_strdup("NA");
4522 } else {
4523 ret = g_strdup_printf("%g", test->pvalue);
4524 }
4525 break;
4526 default:
4527 break;
4528 }
4529
4530 return ret;
4531 }
4532
csv_print_test(const ModelTest * src,PRN * prn)4533 static void csv_print_test (const ModelTest *src, PRN *prn)
4534 {
4535 const char *descrip = NULL;
4536 char c = prn_delim(prn);
4537
4538 descrip = test_get_descrip(src);
4539
4540 if (descrip != NULL) {
4541 const char *optstr = test_get_opt_string(src);
4542
4543 if (optstr != NULL) {
4544 pprintf(prn, "\"%s, %s\"\n", descrip, optstr);
4545 } else {
4546 pprintf(prn, "\"%s\"\n", descrip);
4547 }
4548 }
4549
4550 if (src->param != NULL && *src->param != '\0') {
4551 pprintf(prn, "\"%s\"%c\"%s\"\n", _("parameter"), c, src->param);
4552 }
4553
4554 if (src->dfn > 0 && src->dfd > 0) {
4555 pprintf(prn, "\"%s\"%c%d\n", _("dfn"), c, src->dfn);
4556 pprintf(prn, "\"%s\"%c%g\n", _("dfd"), c, src->dfd);
4557 } else if (src->dfn > 0) {
4558 pprintf(prn, "\"%s\"%c%d\n", _("df"), c, src->dfn);
4559 }
4560 if (src->order) {
4561 pprintf(prn, "\"%s\"%c%d\n", _("lag order"), c, src->order);
4562 }
4563 pprintf(prn, "\"%s\"%c%.15g\n", _("test statistic"), c, src->value);
4564 if (!na(src->pvalue)) {
4565 pprintf(prn, "\"%s\"%c%.15g\n", _("p-value"), c, src->pvalue);
4566 }
4567
4568 if (!na(src->crit)) {
4569 double a = 100 * src->alpha;
4570 gchar *buf;
4571
4572 buf = g_strdup_printf(_("%g percent critical value"), a);
4573 pprintf(prn, "\"%s\"%c%g\n", buf, c, src->crit);
4574 g_free(buf);
4575 }
4576
4577 pputc(prn, '\n');
4578 }
4579
4580 #define asy_test(t) (t == GRETL_STAT_WALD_CHISQ || \
4581 t == GRETL_STAT_Z)
4582
4583 #define asy_pval(test) (test->teststat == GRETL_STAT_SUP_WALD || \
4584 (test->opt & OPT_A))
4585
4586 #define boot_pval(test) (test->opt & OPT_B)
4587
gretl_model_test_print_direct(const ModelTest * test,int heading,PRN * prn)4588 void gretl_model_test_print_direct (const ModelTest *test, int heading, PRN *prn)
4589 {
4590 const char *tstr;
4591 gchar *buf = NULL;
4592
4593 if (rtf_format(prn)) {
4594 pputs(prn, "\\par \\ql ");
4595 }
4596
4597 if (heading) {
4598 gretl_test_print_heading(test, prn);
4599 }
4600
4601 gretl_test_print_h_0(test, heading, prn);
4602
4603 tstr = asy_test(test->teststat)? N_("Asymptotic test statistic") :
4604 N_("Test statistic");
4605
4606 buf = get_test_stat_string(test, prn);
4607
4608 if (plain_format(prn)) {
4609 pprintf(prn, "\n %s: %s\n", _(tstr), buf);
4610 } else if (tex_format(prn)) {
4611 pprintf(prn, "\\\\\n\\quad %s: %s\\\\\n", _(tstr), buf);
4612 } else if (rtf_format(prn)) {
4613 pprintf(prn, "\\par\n %s: %s\\par\n", _(tstr), buf);
4614 }
4615
4616 g_free(buf);
4617 buf = get_test_pval_string(test, prn);
4618
4619 if (buf != NULL) {
4620 const char *pvstr =
4621 asy_pval(test) ? N_("with asymptotic p-value") :
4622 boot_pval(test) ? N_("with bootstrap p-value") :
4623 N_("with p-value");
4624
4625 if (plain_format(prn)) {
4626 pprintf(prn, " %s = %s\n\n", _(pvstr), buf);
4627 } else if (tex_format(prn)) {
4628 pprintf(prn, "\\quad %s = %s\\\\\n", _(pvstr), buf);
4629 } else if (rtf_format(prn)) {
4630 pprintf(prn, " %s = %s\\par\n\n", _(pvstr), buf);
4631 }
4632 } else if (!na(test->crit) && !na(test->alpha)) {
4633 double a = test->alpha * 100.0;
4634
4635 if (plain_format(prn)) {
4636 buf = g_strdup_printf(_("%g percent critical value"), a);
4637 pprintf(prn, " (%s = %.2f)\n\n", buf, test->crit);
4638 } else if (tex_format(prn)) {
4639 buf = g_strdup_printf(_("%g percent critical value"), a);
4640 pprintf(prn, "\\quad (%s = %.2f)\\\\\n", buf, test->crit);
4641 } else if (rtf_format(prn)) {
4642 buf = g_strdup_printf(_("%g percent critical value"), a);
4643 pprintf(prn, " (%s = %.2f)\\par\n\n", buf, test->crit);
4644 }
4645 } else {
4646 pputc(prn, '\n');
4647 }
4648
4649 g_free(buf);
4650 }
4651
gretl_model_test_print(const MODEL * pmod,int i,PRN * prn)4652 void gretl_model_test_print (const MODEL *pmod, int i, PRN *prn)
4653 {
4654 if (i >= 0 && i < pmod->ntests) {
4655 if (csv_format(prn)) {
4656 csv_print_test(&pmod->tests[i], prn);
4657 } else {
4658 gretl_model_test_print_direct(&pmod->tests[i], 1, prn);
4659 }
4660 }
4661 }
4662
gretl_model_print_last_test(const MODEL * pmod,PRN * prn)4663 void gretl_model_print_last_test (const MODEL *pmod, PRN *prn)
4664 {
4665 gretl_model_test_print(pmod, pmod->ntests - 1, prn);
4666 }
4667
copy_ar_info(const ARINFO * src)4668 static ARINFO *copy_ar_info (const ARINFO *src)
4669 {
4670 ARINFO *targ;
4671
4672 targ = malloc(sizeof *targ);
4673 if (targ == NULL) return NULL;
4674
4675 if (src->arlist != NULL) {
4676 int m = src->arlist[0];
4677
4678 if ((targ->rho = copyvec(src->rho, m)) == NULL) {
4679 free(targ);
4680 return NULL;
4681 }
4682
4683 if ((targ->sderr = copyvec(src->sderr, m)) == NULL) {
4684 free(targ->rho);
4685 free(targ);
4686 return NULL;
4687 }
4688
4689 targ->arlist = gretl_list_copy(src->arlist);
4690 if (targ->arlist == NULL) {
4691 free(targ->rho);
4692 free(targ->sderr);
4693 free(targ);
4694 return NULL;
4695 }
4696 }
4697
4698 return targ;
4699 }
4700
copy_model_params(MODEL * targ,const MODEL * src)4701 static int copy_model_params (MODEL *targ, const MODEL *src)
4702 {
4703 int err = 0;
4704
4705 if (src->nparams > 0) {
4706 targ->params = strings_array_dup(src->params, src->nparams);
4707
4708 if (targ->params == NULL) {
4709 err = E_ALLOC;
4710 } else {
4711 targ->nparams = src->nparams;
4712 }
4713 }
4714
4715 return err;
4716 }
4717
serialize_coeff_sep(model_data_item * item,PRN * prn)4718 static void serialize_coeff_sep (model_data_item *item, PRN *prn)
4719 {
4720 CoeffSep *cs = (CoeffSep *) item->ptr;
4721
4722 pprintf(prn, " pos=\"%d\"", cs->pos);
4723 if (*cs->str != '\0') {
4724 pputs(prn, " string=\"");
4725 gretl_xml_put_string(cs->str, prn);
4726 pputc(prn, '"');
4727 }
4728 pputs(prn, "/>\n");
4729 }
4730
serialize_vcv_info(model_data_item * item,PRN * prn)4731 static void serialize_vcv_info (model_data_item *item, PRN *prn)
4732 {
4733 VCVInfo *vi = (VCVInfo *) item->ptr;
4734
4735 pprintf(prn, " vmaj=\"%d\"", vi->vmaj);
4736 pprintf(prn, " vmin=\"%d\"", vi->vmin);
4737 if (vi->order > 0) {
4738 pprintf(prn, " order=\"%d\"", vi->order);
4739 }
4740 if (vi->flags > 0) {
4741 pprintf(prn, " flags=\"%d\"", vi->flags);
4742 }
4743 if (vi->bw > 0 && !na(vi->bw)) {
4744 pprintf(prn, " bw=\"%.12g\"", vi->bw);
4745 }
4746 pputs(prn, "/>\n");
4747 }
4748
4749 /* FIXME updating and placement of these functions */
4750
4751 struct type_mapper {
4752 int type;
4753 const char *name;
4754 const char *compat;
4755 };
4756
4757 static struct type_mapper mapper[] = {
4758 { GRETL_TYPE_INT, "int", "1" },
4759 { GRETL_TYPE_LIST, "list", "2" },
4760 { GRETL_TYPE_DOUBLE, "double", "3" },
4761 { GRETL_TYPE_INT_ARRAY, "intarray", "4" },
4762 { GRETL_TYPE_DOUBLE_ARRAY, "doublearray", "5" },
4763 { GRETL_TYPE_STRING, "string", "6" },
4764 { GRETL_TYPE_CMPLX_ARRAY, "cmplxarray", "8" },
4765 { GRETL_TYPE_STRUCT, "struct" , "9" },
4766 { GRETL_TYPE_MATRIX, "matrix" , "10" },
4767 { GRETL_TYPE_NONE, "none", "0" },
4768 };
4769
type_from_type_string(const char * s)4770 static int type_from_type_string (const char *s)
4771 {
4772 int i;
4773
4774 if (isdigit(*s)) {
4775 /* backward compatibility */
4776 for (i=0; mapper[i].type != GRETL_TYPE_NONE; i++) {
4777 if (*s == mapper[i].compat[0]) {
4778 return mapper[i].type;
4779 }
4780 }
4781 if (*s == '7') {
4782 /* note: was "chararray" */
4783 return GRETL_TYPE_STRING;
4784 }
4785 } else {
4786 for (i=0; mapper[i].type != GRETL_TYPE_NONE; i++) {
4787 if (!strcmp(s, mapper[i].name)) {
4788 return mapper[i].type;
4789 }
4790 }
4791 if (!strcmp(s, "chararray")) {
4792 return GRETL_TYPE_STRING;
4793 }
4794 }
4795
4796 return GRETL_TYPE_NONE;
4797 }
4798
gretl_type_name(int t)4799 static const char *gretl_type_name (int t)
4800 {
4801 int i;
4802
4803 for (i=0; mapper[i].type != GRETL_TYPE_NONE; i++) {
4804 if (t == mapper[i].type) {
4805 return mapper[i].name;
4806 }
4807 }
4808
4809 return "unknown";
4810 }
4811
serialize_model_data_items(const MODEL * pmod,PRN * prn)4812 static void serialize_model_data_items (const MODEL *pmod, PRN *prn)
4813 {
4814 model_data_item *item;
4815 int i, j, nelem;
4816
4817 pprintf(prn, "<data-items count=\"%d\">\n", pmod->n_data_items);
4818
4819 for (i=0; i<pmod->n_data_items; i++) {
4820 item = pmod->data_items[i];
4821 nelem = 0;
4822
4823 pprintf(prn, "<data-item key=\"%s\" type=\"%s\"",
4824 item->key, gretl_type_name(item->type));
4825
4826 if (!strcmp(item->key, "coeffsep")) {
4827 serialize_coeff_sep(item, prn);
4828 continue;
4829 } else if (!strcmp(item->key, "vcv_info")) {
4830 serialize_vcv_info(item, prn);
4831 continue;
4832 }
4833
4834 if (item->type == GRETL_TYPE_INT_ARRAY) {
4835 nelem = item->size / sizeof(int);
4836 } else if (item->type == GRETL_TYPE_DOUBLE_ARRAY) {
4837 nelem = item->size / sizeof(double);
4838 } else if (item->type == GRETL_TYPE_CMPLX_ARRAY) {
4839 nelem = item->size / sizeof(cmplx);
4840 }
4841
4842 if (nelem > 0) {
4843 pprintf(prn, " count=\"%d\">\n", nelem);
4844 } else if (item->type == GRETL_TYPE_STRING) {
4845 pputc(prn, '>');
4846 } else {
4847 pputs(prn, ">\n");
4848 }
4849
4850 if (item->type == GRETL_TYPE_INT) {
4851 pprintf(prn, "%d", *(int *) item->ptr);
4852 } else if (item->type == GRETL_TYPE_DOUBLE) {
4853 double x = *(double *) item->ptr;
4854
4855 if (na(x)) {
4856 pputs(prn, "NA");
4857 } else {
4858 pprintf(prn, "%.15g", x);
4859 }
4860 } else if (item->type == GRETL_TYPE_INT_ARRAY) {
4861 int *vals = (int *) item->ptr;
4862
4863 for (j=0; j<nelem; j++) {
4864 pprintf(prn, "%d ", vals[j]);
4865 }
4866 } else if (item->type == GRETL_TYPE_DOUBLE_ARRAY) {
4867 double *vals = (double *) item->ptr;
4868
4869 for (j=0; j<nelem; j++) {
4870 if (na(vals[j])) {
4871 pputs(prn, "NA ");
4872 } else {
4873 pprintf(prn, "%.15g ", vals[j]);
4874 }
4875 }
4876 } else if (item->type == GRETL_TYPE_CMPLX_ARRAY) {
4877 cmplx *vals = (cmplx *) item->ptr;
4878 double x;
4879
4880 for (j=0; j<nelem; j++) {
4881 x = vals[j].r;
4882 if (na(x)) {
4883 pputs(prn, "NA ");
4884 } else {
4885 pprintf(prn, "%.15g ", x);
4886 }
4887 x = vals[j].i;
4888 if (na(x)) {
4889 pputs(prn, "NA ");
4890 } else {
4891 pprintf(prn, "%.15g ", x);
4892 }
4893 }
4894 } else if (item->type == GRETL_TYPE_LIST) {
4895 int *list = (int *) item->ptr;
4896
4897 for (j=0; j<=list[0]; j++) {
4898 pprintf(prn, "%d ", list[j]);
4899 }
4900 } else if (item->type == GRETL_TYPE_STRING) {
4901 gretl_xml_put_string((const char *) item->ptr, prn);
4902 } else if (item->type == GRETL_TYPE_MATRIX) {
4903 gretl_matrix *m = (gretl_matrix *) item->ptr;
4904
4905 gretl_matrix_serialize(m, NULL, prn);
4906 } else if (item->type == GRETL_TYPE_ARRAY) {
4907 gretl_array *A = (gretl_array *) item->ptr;
4908
4909 gretl_array_serialize(A, prn);
4910 } else {
4911 ; /* no-op: not handled */
4912 }
4913
4914 pputs(prn, "</data-item>\n");
4915 }
4916
4917 pputs(prn, "</data-items>\n");
4918 }
4919
copy_model_data_items(MODEL * targ,const MODEL * src)4920 static int copy_model_data_items (MODEL *targ, const MODEL *src)
4921 {
4922 int i, n = src->n_data_items;
4923 int err = 0;
4924
4925 targ->data_items = malloc(n * sizeof *targ->data_items);
4926 if (targ->data_items == NULL) {
4927 return 1;
4928 }
4929
4930 for (i=0; i<n; i++) {
4931 targ->data_items[i] = NULL;
4932 }
4933
4934 for (i=0; i<n; i++) {
4935 targ->data_items[i] = replicate_data_item(src->data_items[i]);
4936 if (targ->data_items[i] == NULL) {
4937 err = 1;
4938 break;
4939 }
4940 }
4941
4942 if (err) {
4943 for (i=0; i<n; i++) {
4944 free(targ->data_items[i]);
4945 }
4946 free(targ->data_items);
4947 targ->data_items = NULL;
4948 targ->n_data_items = 0;
4949 } else {
4950 targ->n_data_items = n;
4951 }
4952
4953 return err;
4954 }
4955
4956 /* ensure we get a valid identifier for the bundle key */
4957
item_key_to_bundle_key(char * targ,const char * src)4958 static char *item_key_to_bundle_key (char *targ,
4959 const char *src)
4960 {
4961 strcpy(targ, src);
4962 return gretl_charsub(targ, '-', '_');
4963 }
4964
my_matrix_from_list(model_data_item * item)4965 static gretl_matrix *my_matrix_from_list (model_data_item *item)
4966 {
4967 int *list = item->ptr;
4968 gretl_matrix *m;
4969 int i;
4970
4971 m = gretl_matrix_alloc(1, list[0]);
4972
4973 if (m != NULL) {
4974 for (i=0; i<list[0]; i++) {
4975 m->val[i] = list[i+1];
4976 }
4977 }
4978
4979 return m;
4980 }
4981
matrix_from_cmplx(model_data_item * item)4982 static gretl_matrix *matrix_from_cmplx (model_data_item *item)
4983 {
4984 int nr = item->size / sizeof(cmplx);
4985 cmplx *c = item->ptr;
4986 gretl_matrix *m;
4987 int i;
4988
4989 m = gretl_matrix_alloc(nr, 2);
4990
4991 if (m != NULL) {
4992 for (i=0; i<nr; i++) {
4993 gretl_matrix_set(m, i, 0, c[i].r);
4994 gretl_matrix_set(m, i, 1, c[i].i);
4995 }
4996 }
4997
4998 return m;
4999 }
5000
bundlize_model_data_items(const MODEL * pmod,gretl_bundle * b)5001 int bundlize_model_data_items (const MODEL *pmod, gretl_bundle *b)
5002 {
5003 model_data_item *item;
5004 char bkey[VNAMELEN];
5005 gretl_matrix *m;
5006 double xval;
5007 int i, ival;
5008 int err = 0;
5009
5010 for (i=0; i<pmod->n_data_items && !err; i++) {
5011 item = pmod->data_items[i];
5012 item_key_to_bundle_key(bkey, item->key);
5013 if (gretl_bundle_has_key(b, bkey)) {
5014 continue; /* item already present */
5015 } else if (item->type == GRETL_TYPE_INT) {
5016 ival = *(int *) item->ptr;
5017 err = gretl_bundle_set_scalar(b, bkey, ival);
5018 } else if (item->type == GRETL_TYPE_DOUBLE) {
5019 xval = *(double *) item->ptr;
5020 err = gretl_bundle_set_scalar(b, bkey, xval);
5021 } else if (item->type == GRETL_TYPE_MATRIX) {
5022 err = gretl_bundle_set_matrix(b, bkey, item->ptr);
5023 } else if (item->type == GRETL_TYPE_LIST) {
5024 if (pmod->ci == MIDASREG && !strcmp(bkey, "seplist")) {
5025 ; /* internal use only, skip it */
5026 } else if (pmod->ci == ARMA && !strcmp(bkey, "ainfo")) {
5027 m = my_matrix_from_list(item);
5028 if (m != NULL) {
5029 err = gretl_bundle_donate_data(b, bkey, m,
5030 GRETL_TYPE_MATRIX,
5031 0);
5032 }
5033 } else {
5034 err = gretl_bundle_set_list(b, bkey, item->ptr);
5035 }
5036 } else if (item->type == GRETL_TYPE_ARRAY) {
5037 gretl_bundle_set_data(b, bkey, item->ptr, item->type, 0);
5038 } else if (item->type == GRETL_TYPE_CMPLX_ARRAY) {
5039 m = matrix_from_cmplx(item);
5040 if (m != NULL) {
5041 err = gretl_bundle_donate_data(b, bkey, m,
5042 GRETL_TYPE_MATRIX,
5043 0);
5044 }
5045 }
5046 }
5047
5048 if (!err && pmod->tests != NULL) {
5049 const ModelTest *test;
5050 const char *tkey;
5051 gretl_bundle *bt;
5052
5053 for (i=0; i<pmod->ntests && !err; i++) {
5054 test = &pmod->tests[i];
5055 tkey = test_type_key(test->type);
5056 if (tkey != NULL && !gretl_bundle_has_key(b, tkey)) {
5057 bt = bundlize_test(test);
5058 if (bt != NULL) {
5059 err = gretl_bundle_donate_data(b, tkey, bt,
5060 GRETL_TYPE_BUNDLE, 0);
5061 }
5062 }
5063 }
5064 }
5065
5066 if (!err) {
5067 /* estimation start and end, 1-based */
5068 gretl_bundle_set_int(b, "t1", pmod->t1 + 1);
5069 gretl_bundle_set_int(b, "t2", pmod->t2 + 1);
5070 /* contemporaneous sample range, 1-based */
5071 gretl_bundle_set_int(b, "sample_t1", pmod->smpl.t1 + 1);
5072 gretl_bundle_set_int(b, "sample_t2", pmod->smpl.t2 + 1);
5073 }
5074
5075 if (pmod->ci == PANEL && (pmod->opt & OPT_F) &&
5076 !gretl_bundle_has_key(b, "within_rsq")) {
5077 /* fixed effects: add within R-squared */
5078 gretl_bundle_set_scalar(b, "within_rsq", pmod->adjrsq);
5079 }
5080
5081 if (!na(pmod->dw) && !gretl_bundle_has_key(b, "dw")) {
5082 /* add Durbin-Watson statistic if available */
5083 gretl_bundle_set_scalar(b, "dw", pmod->dw);
5084 }
5085
5086 return err;
5087 }
5088
display_model_data_items(const MODEL * pmod)5089 void display_model_data_items (const MODEL *pmod)
5090 {
5091 model_data_item *item;
5092 int i, n = pmod->n_data_items;
5093
5094 fprintf(stderr, "model has %d data items\n", n);
5095
5096 for (i=0; i<n; i++) {
5097 item = pmod->data_items[i];
5098 fprintf(stderr, "%d '%s': ", i, item->key);
5099 if (item->type == GRETL_TYPE_INT) {
5100 fprintf(stderr, "%d\n", *(int *) item->ptr);
5101 } else if (item->type == GRETL_TYPE_DOUBLE) {
5102 fprintf(stderr, "%.15g\n", *(double *) item->ptr);
5103 } else {
5104 fprintf(stderr, "%p\n", (void *) item->ptr);
5105 }
5106 }
5107 }
5108
copy_model(MODEL * targ,MODEL * src)5109 static int copy_model (MODEL *targ, MODEL *src)
5110 {
5111 int k = src->ncoeff;
5112 int m = k * (k + 1) / 2;
5113 int err = 0;
5114
5115 /* monolithic copy of structure */
5116 *targ = *src;
5117
5118 /* temporarily zero various count members just in case
5119 copying of allocated info fails */
5120 targ->ntests = 0;
5121 targ->nparams = 0;
5122 targ->n_data_items = 0;
5123
5124 /* now work on pointer members */
5125 gretl_model_init_pointers(targ);
5126
5127 if (src->coeff != NULL &&
5128 (targ->coeff = copyvec(src->coeff, src->ncoeff)) == NULL) {
5129 return 1;
5130 }
5131
5132 if (src->sderr != NULL &&
5133 (targ->sderr = copyvec(src->sderr, src->ncoeff)) == NULL) {
5134 return 1;
5135 }
5136
5137 if (src->uhat != NULL &&
5138 (targ->uhat = copyvec(src->uhat, src->full_n)) == NULL) {
5139 return 1;
5140 }
5141
5142 if (src->yhat != NULL &&
5143 (targ->yhat = copyvec(src->yhat, src->full_n)) == NULL) {
5144 return 1;
5145 }
5146
5147 if (src->submask != NULL &&
5148 (targ->submask = copy_subsample_mask(src->submask, &err)) == NULL) {
5149 return 1;
5150 }
5151
5152 if (src->missmask != NULL &&
5153 (targ->missmask = gretl_strdup(src->missmask)) == NULL) {
5154 return 1;
5155 }
5156
5157 if (src->xpx != NULL && (targ->xpx = copyvec(src->xpx, m)) == NULL) {
5158 return 1;
5159 }
5160
5161 if (src->vcv != NULL && (targ->vcv = copyvec(src->vcv, m)) == NULL) {
5162 return 1;
5163 }
5164
5165 if (src->arinfo != NULL &&
5166 (targ->arinfo = copy_ar_info(src->arinfo)) == NULL) {
5167 return 1;
5168 }
5169
5170 if (src->ntests > 0 && src->tests != NULL) {
5171 copy_model_tests(targ, src);
5172 if (targ->tests == NULL) {
5173 return 1;
5174 }
5175 targ->ntests = src->ntests;
5176 }
5177
5178 if (src->name != NULL) {
5179 targ->name = gretl_strdup(src->name);
5180 if (targ->name == NULL) {
5181 return 1;
5182 }
5183 }
5184
5185 if (src->depvar != NULL) {
5186 targ->depvar = gretl_strdup(src->depvar);
5187 if (targ->depvar == NULL) {
5188 return 1;
5189 }
5190 }
5191
5192 if (src->nparams > 0 && src->params != NULL) {
5193 err = copy_model_params(targ, src);
5194 if (err) {
5195 return err;
5196 }
5197 }
5198
5199 if (src->n_data_items > 0) {
5200 copy_model_data_items(targ, src);
5201 if (targ->data_items == NULL) {
5202 return 1;
5203 }
5204 targ->n_data_items = src->n_data_items;
5205 }
5206
5207 if (src->list != NULL &&
5208 (targ->list = gretl_list_copy(src->list)) == NULL) {
5209 return 1;
5210 }
5211
5212 if (gretl_is_between_model(src) && src->dataset != NULL) {
5213 /* special for "between" panel model: transfer the
5214 group-means dataset to @targ
5215 */
5216 targ->dataset = src->dataset;
5217 src->dataset = NULL;
5218 }
5219
5220 return 0;
5221 }
5222
gretl_model_serialize(const MODEL * pmod,SavedObjectFlags flags,PRN * prn)5223 int gretl_model_serialize (const MODEL *pmod, SavedObjectFlags flags,
5224 PRN *prn)
5225 {
5226 char *xmlname = NULL;
5227 int k = pmod->ncoeff;
5228 int m = k * (k + 1) / 2;
5229 int err = 0;
5230
5231 if (pmod->name == NULL || *pmod->name == '\0') {
5232 xmlname = gretl_strdup("none");
5233 } else {
5234 xmlname = gretl_xml_encode(pmod->name);
5235 }
5236
5237 pprintf(prn, "<gretl-model ID=\"%d\" name=\"%s\" saveflags=\"%d\" ",
5238 pmod->ID, xmlname, (int) flags);
5239 free(xmlname);
5240
5241 if (pmod->depvar != NULL) {
5242 pprintf(prn, "depvar=\"%s\" ", pmod->depvar);
5243 }
5244
5245 pprintf(prn, "t1=\"%d\" t2=\"%d\" nobs=\"%d\" ",
5246 pmod->t1, pmod->t2, pmod->nobs);
5247 pprintf(prn, "full_n=\"%d\" ncoeff=\"%d\" dfn=\"%d\" dfd=\"%d\" ",
5248 pmod->full_n, pmod->ncoeff, pmod->dfn, pmod->dfd);
5249 pprintf(prn, "ifc=\"%d\" ci=\"%s\" opt=\"%u\" nwt=\"%d\" aux=\"%d\" ",
5250 pmod->ifc, gretl_command_word(pmod->ci), pmod->opt,
5251 pmod->nwt, pmod->aux);
5252
5253 gretl_push_c_numeric_locale();
5254
5255 gretl_xml_put_double("ess", pmod->ess, prn);
5256 gretl_xml_put_double("tss", pmod->tss, prn);
5257 gretl_xml_put_double("sigma", pmod->sigma, prn);
5258 gretl_xml_put_double("rsq", pmod->rsq, prn);
5259 gretl_xml_put_double("adjrsq", pmod->adjrsq, prn);
5260 gretl_xml_put_double("fstt", pmod->fstt, prn);
5261 gretl_xml_put_double("chi-square", pmod->chisq, prn);
5262 gretl_xml_put_double("lnL", pmod->lnL, prn);
5263 gretl_xml_put_double("ybar", pmod->ybar, prn);
5264 gretl_xml_put_double("sdy", pmod->sdy, prn);
5265
5266 gretl_xml_put_double("crit0", pmod->criterion[0], prn);
5267 gretl_xml_put_double("crit1", pmod->criterion[1], prn);
5268 gretl_xml_put_double("crit2", pmod->criterion[2], prn);
5269
5270 gretl_xml_put_double("dw", pmod->dw, prn);
5271 gretl_xml_put_double("rho", pmod->rho, prn);
5272
5273 pputs(prn, ">\n");
5274
5275 if (pmod->smpl.rseed > 0) {
5276 pprintf(prn, "<sample t1=\"%d\" t2=\"%d\" rseed=\"%u\"/>\n",
5277 pmod->smpl.t1, pmod->smpl.t2, pmod->smpl.rseed);
5278 } else {
5279 pprintf(prn, "<sample t1=\"%d\" t2=\"%d\"/>\n",
5280 pmod->smpl.t1, pmod->smpl.t2);
5281 }
5282
5283 gretl_xml_put_double_array("coeff", pmod->coeff, k, prn);
5284 gretl_xml_put_double_array("sderr", pmod->sderr, k, prn);
5285
5286 if (pmod->uhat != NULL) {
5287 gretl_xml_put_double_array("uhat", pmod->uhat, pmod->full_n, prn);
5288 }
5289
5290 if (pmod->yhat != NULL) {
5291 gretl_xml_put_double_array("yhat", pmod->yhat, pmod->full_n, prn);
5292 }
5293
5294 if (pmod->submask != NULL) {
5295 write_model_submask(pmod, prn);
5296 }
5297
5298 if (pmod->missmask != NULL) {
5299 pputs(prn, "<missmask>");
5300 pputs(prn, pmod->missmask);
5301 pputs(prn, "</missmask>\n");
5302 }
5303
5304 if (pmod->xpx != NULL) {
5305 gretl_xml_put_double_array("xpx", pmod->xpx, m, prn);
5306 }
5307
5308 if (pmod->vcv != NULL) {
5309 gretl_xml_put_double_array("vcv", pmod->vcv, m, prn);
5310 }
5311
5312 if (pmod->arinfo != NULL && pmod->arinfo->arlist != NULL) {
5313 int r = pmod->arinfo->arlist[0];
5314
5315 pputs(prn, "<arinfo>\n");
5316 gretl_xml_put_tagged_list("arlist", pmod->arinfo->arlist, prn);
5317 gretl_xml_put_double_array("rho", pmod->arinfo->rho, r, prn);
5318 gretl_xml_put_double_array("sderr", pmod->arinfo->sderr, r, prn);
5319 pputs(prn, "</arinfo>\n");
5320 }
5321
5322 if (pmod->ntests > 0 && pmod->tests != NULL) {
5323 serialize_model_tests(pmod, prn);
5324 }
5325
5326 if (pmod->nparams > 0 && pmod->params != NULL) {
5327 gretl_xml_put_strings_array("params",
5328 (const char **) pmod->params,
5329 pmod->nparams, prn);
5330 }
5331
5332 if (pmod->list != NULL) {
5333 gretl_xml_put_tagged_list("list", pmod->list, prn);
5334 }
5335
5336 if (pmod->n_data_items > 0) {
5337 serialize_model_data_items(pmod, prn);
5338 }
5339
5340 pputs(prn, "</gretl-model>\n");
5341
5342 /* note: the DATASET element of pmod is not
5343 handled here */
5344
5345 gretl_pop_c_numeric_locale();
5346
5347 return err;
5348 }
5349
5350 /* next block: functions for reconstituting a model from
5351 its XML representation */
5352
model_submask_from_xml(xmlNodePtr node,xmlDocPtr doc,MODEL * pmod)5353 static int model_submask_from_xml (xmlNodePtr node, xmlDocPtr doc,
5354 MODEL *pmod)
5355 {
5356 char *mask;
5357 int err;
5358
5359 err = gretl_xml_get_submask(node, doc, &mask);
5360 if (!err) {
5361 pmod->submask = mask;
5362 }
5363
5364 return err;
5365 }
5366
5367 static int
retrieve_model_coeff_separator(xmlNodePtr cur,MODEL * pmod)5368 retrieve_model_coeff_separator (xmlNodePtr cur, MODEL *pmod)
5369 {
5370 CoeffSep *cs = malloc(sizeof *cs);
5371 char *tmp = NULL;
5372 int err = 0;
5373
5374 if (cs != NULL) {
5375 cs->str[0] = '\0';
5376 gretl_xml_get_prop_as_int(cur, "pos", &cs->pos);
5377 gretl_xml_get_prop_as_string(cur, "string", &tmp);
5378 if (tmp != NULL) {
5379 strcpy(cs->str, tmp);
5380 free(tmp);
5381 }
5382 err = gretl_model_set_data(pmod, "coeffsep", cs,
5383 GRETL_TYPE_STRUCT,
5384 sizeof *cs);
5385 }
5386
5387 return err;
5388 }
5389
5390 static int
retrieve_model_vcv_info(xmlNodePtr cur,MODEL * pmod)5391 retrieve_model_vcv_info (xmlNodePtr cur, MODEL *pmod)
5392 {
5393 VCVInfo *vi = vcv_info_new();
5394 int err = 0;
5395
5396 if (vi != NULL) {
5397 int ival;
5398 double x;
5399
5400 gretl_xml_get_prop_as_int(cur, "vmaj", &vi->vmaj);
5401 gretl_xml_get_prop_as_int(cur, "vmin", &vi->vmin);
5402 if (gretl_xml_get_prop_as_int(cur, "order", &ival)) {
5403 vi->order = ival;
5404 }
5405 if (gretl_xml_get_prop_as_int(cur, "flags", &ival)) {
5406 vi->flags = ival;
5407 }
5408 if (gretl_xml_get_prop_as_double(cur, "bw", &x)) {
5409 vi->bw = x;
5410 }
5411 err = gretl_model_set_data(pmod, "vcv_info", vi,
5412 GRETL_TYPE_STRUCT,
5413 sizeof *vi);
5414 }
5415
5416 return err;
5417 }
5418
5419 /* backward compatibility: transcribe old-style "gretl_set_int"
5420 settings to bits in pmod->opt, where applicable
5421 */
5422
gretl_model_set_int_compat(MODEL * pmod,const char * key,int ival)5423 static int gretl_model_set_int_compat (MODEL *pmod,
5424 const char *key,
5425 int ival)
5426 {
5427 gretlopt opt = pmod->opt;
5428 int err = 0;
5429
5430 if (ival == 0) {
5431 return 0;
5432 }
5433
5434 if (pmod->ci == PANEL) {
5435 if (!strcmp(key, "between")) {
5436 pmod->opt |= OPT_B;
5437 } else if (!strcmp(key, "fixed-effects")) {
5438 pmod->opt |= OPT_F;
5439 } else if (!strcmp(key, "random-effects")) {
5440 pmod->opt |= OPT_U;
5441 } else if (!strcmp(key, "unit-weights")) {
5442 pmod->opt |= OPT_H;
5443 }
5444 } else if (pmod->ci == AR1) {
5445 if (!strcmp(key, "hilu")) {
5446 pmod->opt |= OPT_H;
5447 } else if (!strcmp(key, "pwe")) {
5448 pmod->opt |= OPT_P;
5449 }
5450 } else if (pmod->ci == GMM) {
5451 if (!strcmp(key, "iterated")) {
5452 pmod->opt |= OPT_I;
5453 } else if (!strcmp(key, "two-step")) {
5454 pmod->opt |= OPT_T;
5455 }
5456 } else if (pmod->ci == HECKIT) {
5457 if (!strcmp(key, "two-step")) {
5458 pmod->opt |= OPT_T;
5459 }
5460 } else if (pmod->ci == LOGIT || pmod->ci == PROBIT) {
5461 if (!strcmp(key, "show-pvals")) {
5462 pmod->opt |= OPT_P;
5463 }
5464 }
5465
5466 if (!strcmp(key, "robust")) {
5467 pmod->opt |= OPT_R;
5468 } else if (!strcmp(key, "time-dummies")) {
5469 pmod->opt |= OPT_D;
5470 } else if (!strcmp(key, "no-df-corr")) {
5471 pmod->opt |= OPT_N;
5472 }
5473
5474 if (pmod->opt == opt) {
5475 /* not handled above */
5476 err = gretl_model_set_int(pmod, key, ival);
5477 }
5478
5479 return err;
5480 }
5481
5482 /* backward compat: convert ad hoc settings to VCVInfo */
5483
compat_compose_vcv_info(MODEL * pmod)5484 static void compat_compose_vcv_info (MODEL *pmod)
5485 {
5486 VCVInfo vi = { 0, 0, 0, 0, NADBL };
5487 int ival;
5488
5489 if (gretl_model_get_int(pmod, "hc")) {
5490 vi.vmaj = VCV_HC;
5491 vi.vmin = gretl_model_get_int(pmod, "hc_version");
5492 } else if ((ival = gretl_model_get_int(pmod, "ml_vcv"))) {
5493 vi.vmaj = VCV_ML;
5494 vi.vmin = ival;
5495 } else if (gretl_model_get_int(pmod, "panel_hac")) {
5496 vi.vmaj = VCV_PANEL;
5497 vi.vmin = PANEL_HAC;
5498 } else if (gretl_model_get_int(pmod, "panel_bk")) {
5499 vi.vmaj = VCV_PANEL;
5500 vi.vmin = PANEL_BK;
5501 } else if (gretl_model_get_int(pmod, "using_hac") ||
5502 gretl_model_get_int(pmod, "hac_kernel") ||
5503 gretl_model_get_int(pmod, "hac_lag")) {
5504 vi.vmaj = VCV_HAC;
5505 vi.vmin = gretl_model_get_int(pmod, "hac_kernel");
5506 vi.order = gretl_model_get_int(pmod, "hac_lag");
5507 vi.flags = gretl_model_get_int(pmod, "hac_prewhiten");
5508 if (vi.vmin == KERNEL_QS) {
5509 vi.bw = gretl_model_get_double(pmod, "qs_bandwidth");
5510 }
5511 } else if (pmod->ci == LAD && gretl_model_get_int(pmod, "rq")) {
5512 double a = gretl_model_get_double(pmod, "rq_alpha");
5513
5514 if (na(a)) {
5515 /* doing VCV, not confidence intervals */
5516 vi.vmaj = VCV_RQ;
5517 vi.vmin = (gretl_model_get_int(pmod, "rq_nid"))?
5518 RQ_NID : RQ_ASY;
5519 }
5520 }
5521
5522 if (vi.vmaj > 0) {
5523 VCVInfo *pvi = vcv_info_new();
5524
5525 if (pvi != NULL) {
5526 *pvi = vi;
5527 gretl_model_set_data(pmod, "vcv_info", pvi,
5528 GRETL_TYPE_STRUCT,
5529 sizeof *pvi);
5530 }
5531 }
5532 }
5533
5534 /* backward compat: convert old list separator */
5535
5536 #define old_listsep(v) (v == 999 || v == 9999)
5537
maybe_convert_listsep(int * list,const DATASET * dset)5538 static void maybe_convert_listsep (int *list, const DATASET *dset)
5539 {
5540 int i;
5541
5542 for (i=1; i<=list[0]; i++) {
5543 if (old_listsep(list[i]) && list[i] >= dset->v) {
5544 list[i] = LISTSEP;
5545 }
5546 }
5547 }
5548
data_list_to_matrix(const int * list)5549 static gretl_matrix *data_list_to_matrix (const int *list)
5550 {
5551 gretl_matrix *y = gretl_column_vector_alloc(list[0]);
5552 int i;
5553
5554 if (y != NULL) {
5555 for (i=0; i<y->rows; i++) {
5556 y->val[i] = list[i+1];
5557 }
5558 }
5559
5560 return y;
5561 }
5562
5563 #define XDEBUG 1
5564
model_data_items_from_xml(xmlNodePtr node,xmlDocPtr doc,MODEL * pmod,int fixopt)5565 static int model_data_items_from_xml (xmlNodePtr node, xmlDocPtr doc,
5566 MODEL *pmod, int fixopt)
5567 {
5568 xmlNodePtr cur;
5569 int n_items;
5570 int got_vcv = 0;
5571 int err = 0;
5572
5573 if (!gretl_xml_get_prop_as_int(node, "count", &n_items)) {
5574 return 1;
5575 }
5576
5577 cur = node->xmlChildrenNode;
5578
5579 while (cur != NULL && !err) {
5580 char *key = NULL;
5581 char *typestr = NULL;
5582 int ignore_err = 0;
5583 int t, nelem = 0;
5584
5585 if (!gretl_xml_get_prop_as_string(cur, "type", &typestr) ||
5586 !gretl_xml_get_prop_as_string(cur, "key", &key)) {
5587 err = 1;
5588 break;
5589 }
5590
5591 #if XDEBUG
5592 fprintf(stderr, "data_item: type='%s', key='%s'\n", typestr, key);
5593 #endif
5594
5595 t = type_from_type_string(typestr);
5596
5597 if (!strcmp(key, "roots")) {
5598 /* could be a backward-compat problem */
5599 ignore_err = 1;
5600 }
5601
5602 if (!strcmp(key, "coeffsep") || !strcmp(key, "10")) {
5603 /* special, with backward compatibility */
5604 err = retrieve_model_coeff_separator(cur, pmod);
5605 } else if (!strcmp(key, "vcv_info")) {
5606 err = retrieve_model_vcv_info(cur, pmod);
5607 if (!err) {
5608 got_vcv = 1;
5609 }
5610 } else if (t == GRETL_TYPE_INT) {
5611 int ival;
5612
5613 if (!gretl_xml_node_get_int(cur, doc, &ival)) {
5614 err = 1;
5615 } else if (fixopt) {
5616 err = gretl_model_set_int_compat(pmod, key, ival);
5617 } else {
5618 err = gretl_model_set_int(pmod, key, ival);
5619 }
5620 } else if (t == GRETL_TYPE_DOUBLE) {
5621 double xval;
5622
5623 if (!gretl_xml_node_get_double(cur, doc, &xval)) {
5624 err = 1;
5625 } else {
5626 err = gretl_model_set_double(pmod, key, xval);
5627 }
5628 } else if (t == GRETL_TYPE_LIST) {
5629 if (!strcmp(key, "xlist")) {
5630 /* ad hoc (for forecasting): will be recreated if need be */
5631 ;
5632 } else if (!strcmp(key, "yvals")) {
5633 /* backward compatibility for mnlogit models */
5634 int *list = gretl_xml_get_list(cur, doc, &err);
5635
5636 if (!err && list != NULL) {
5637 gretl_matrix *yv = data_list_to_matrix(list);
5638
5639 if (yv != NULL) {
5640 err = gretl_model_set_matrix_as_data(pmod, key, yv);
5641 }
5642 }
5643 free(list);
5644 } else {
5645 int *list = gretl_xml_get_list(cur, doc, &err);
5646
5647 if (!err && list != NULL) {
5648 err = gretl_model_set_list_as_data(pmod, key, list);
5649 }
5650 }
5651 } else if (t == GRETL_TYPE_STRING) {
5652 char *s;
5653 int gots;
5654
5655 if (!strcmp(key, "pmask") || !strcmp(key, "qmask")) {
5656 /* these masks must not have leading white space */
5657 gots = gretl_xml_node_get_trimmed_string(cur, doc, &s);
5658 } else {
5659 gots = gretl_xml_node_get_string(cur, doc, &s);
5660 }
5661 if (!gots) {
5662 err = 1;
5663 } else {
5664 err = gretl_model_set_string_as_data(pmod, key, s);
5665 }
5666 } else if (t == GRETL_TYPE_INT_ARRAY) {
5667 int *ivals = gretl_xml_get_int_array(cur, doc, &nelem, &err);
5668
5669 if (nelem > 0) {
5670 err = gretl_model_set_data(pmod, key, ivals, t,
5671 nelem * sizeof *ivals);
5672 }
5673 } else if (t == GRETL_TYPE_DOUBLE_ARRAY) {
5674 double *xvals = gretl_xml_get_double_array(cur, doc, &nelem, &err);
5675
5676 if (nelem > 0) {
5677 err = gretl_model_set_data(pmod, key, xvals, t,
5678 nelem * sizeof *xvals);
5679 }
5680 } else if (t == GRETL_TYPE_CMPLX_ARRAY) {
5681 cmplx *cvals = gretl_xml_get_cmplx_array(cur, doc, &nelem, &err);
5682
5683 if (nelem > 0) {
5684 err = gretl_model_set_data(pmod, key, cvals, t,
5685 nelem * sizeof *cvals);
5686 }
5687 } else if (t == GRETL_TYPE_MATRIX) {
5688 xmlNodePtr child = cur->xmlChildrenNode;
5689 gretl_matrix *m;
5690
5691 if (child == NULL) {
5692 err = E_DATA;
5693 } else {
5694 m = gretl_xml_get_matrix(child, doc, &err);
5695 if (!err) {
5696 err = gretl_model_set_matrix_as_data(pmod, key, m);
5697 }
5698 }
5699 }
5700
5701 if (err && ignore_err) {
5702 fprintf(stderr, "ignoring error reloading %s\n", key);
5703 err = 0;
5704 }
5705
5706 xmlFree(key);
5707 xmlFree(typestr);
5708
5709 cur = cur->next;
5710 }
5711
5712 if (!err && !got_vcv) {
5713 compat_compose_vcv_info(pmod);
5714 }
5715
5716 #if XDEBUG
5717 fprintf(stderr, "data items from XML, returning %d\n", err);
5718 #endif
5719
5720 return err;
5721 }
5722
arinfo_from_xml(xmlNodePtr node,xmlDocPtr doc,MODEL * pmod)5723 static int arinfo_from_xml (xmlNodePtr node, xmlDocPtr doc,
5724 MODEL *pmod)
5725 {
5726 xmlNodePtr cur;
5727 int n, err = 0;
5728
5729 if (gretl_model_add_arinfo(pmod, 0)) {
5730 return 1;
5731 }
5732
5733 cur = node->xmlChildrenNode;
5734
5735 while (cur != NULL && !err) {
5736 if (!xmlStrcmp(cur->name, (XUC) "arlist")) {
5737 pmod->arinfo->arlist = gretl_xml_get_list(cur, doc, &err);
5738 } else if (!xmlStrcmp(cur->name, (XUC) "rho")) {
5739 pmod->arinfo->rho = gretl_xml_get_double_array(cur, doc, &n, &err);
5740 } else if (!xmlStrcmp(cur->name, (XUC) "sderr")) {
5741 pmod->arinfo->sderr = gretl_xml_get_double_array(cur, doc, &n, &err);
5742 }
5743 cur = cur->next;
5744 }
5745
5746 return err;
5747 }
5748
5749 /**
5750 * gretl_model_from_XML:
5751 * @node: XML node from which to read.
5752 * @doc: pointer to XML document.
5753 * @dset: dataset information.
5754 * @err: location to receive error code.
5755 *
5756 * Reads info on a gretl model from the given XML @node
5757 * and @doc, and reconstitutes the model in memory.
5758 *
5759 * Returns: allocated model, or %NULL on failure.
5760 */
5761
gretl_model_from_XML(xmlNodePtr node,xmlDocPtr doc,const DATASET * dset,int * err)5762 MODEL *gretl_model_from_XML (xmlNodePtr node, xmlDocPtr doc,
5763 const DATASET *dset,
5764 int *err)
5765 {
5766 MODEL *pmod;
5767 char *buf = NULL;
5768 char *modname = NULL;
5769 xmlNodePtr cur;
5770 int fixopt = 0;
5771 int n, got = 0;
5772
5773 pmod = gretl_model_new();
5774 if (pmod == NULL) {
5775 *err = E_ALLOC;
5776 return NULL;
5777 }
5778
5779 got += gretl_xml_get_prop_as_int(node, "ID", &pmod->ID);
5780 got += gretl_xml_get_prop_as_int(node, "t1", &pmod->t1);
5781 got += gretl_xml_get_prop_as_int(node, "t2", &pmod->t2);
5782 got += gretl_xml_get_prop_as_int(node, "nobs", &pmod->nobs);
5783 got += gretl_xml_get_prop_as_int(node, "full_n", &pmod->full_n);
5784 got += gretl_xml_get_prop_as_int(node, "ncoeff", &pmod->ncoeff);
5785 got += gretl_xml_get_prop_as_int(node, "dfn", &pmod->dfn);
5786 got += gretl_xml_get_prop_as_int(node, "dfd", &pmod->dfd);
5787 got += gretl_xml_get_prop_as_int(node, "ifc", &pmod->ifc);
5788
5789 got += gretl_xml_get_prop_as_string(node, "ci", &buf);
5790
5791 if (got < 10) {
5792 fprintf(stderr, "gretl_model_from_XML: got(1) = %d (expected 12)\n", got);
5793 if (buf != NULL) {
5794 free(buf);
5795 }
5796 *err = E_DATA;
5797 goto bailout;
5798 }
5799
5800 gretl_xml_get_prop_as_int(node, "nwt", &pmod->nwt);
5801 gretl_xml_get_prop_as_int(node, "aux", &pmod->aux);
5802
5803 if (gretl_xml_get_prop_as_int(node, "opt", &n)) {
5804 pmod->opt = (unsigned) n;
5805 } else {
5806 /* backward compat may be needed */
5807 fixopt = 1;
5808 }
5809
5810 if (isdigit(*buf)) {
5811 /* backward compatibility: the model command is given
5812 as a number, not a string */
5813 pmod->ci = atoi(buf);
5814 } else {
5815 pmod->ci = gretl_command_number(buf);
5816 }
5817
5818 free(buf);
5819
5820 gretl_push_c_numeric_locale();
5821
5822 gretl_xml_get_prop_as_double(node, "ess", &pmod->ess);
5823 gretl_xml_get_prop_as_double(node, "tss", &pmod->tss);
5824 gretl_xml_get_prop_as_double(node, "sigma", &pmod->sigma);
5825 gretl_xml_get_prop_as_double(node, "rsq", &pmod->rsq);
5826 gretl_xml_get_prop_as_double(node, "adjrsq", &pmod->adjrsq);
5827 gretl_xml_get_prop_as_double(node, "fstt", &pmod->fstt);
5828 gretl_xml_get_prop_as_double(node, "chi-square", &pmod->chisq);
5829 gretl_xml_get_prop_as_double(node, "lnL", &pmod->lnL);
5830 gretl_xml_get_prop_as_double(node, "ybar", &pmod->ybar);
5831 gretl_xml_get_prop_as_double(node, "sdy", &pmod->sdy);
5832
5833 gretl_xml_get_prop_as_double(node, "crit0", &pmod->criterion[0]);
5834 gretl_xml_get_prop_as_double(node, "crit1", &pmod->criterion[1]);
5835 gretl_xml_get_prop_as_double(node, "crit2", &pmod->criterion[2]);
5836
5837 gretl_xml_get_prop_as_double(node, "dw", &pmod->dw);
5838 gretl_xml_get_prop_as_double(node, "rho", &pmod->rho);
5839
5840 gretl_xml_get_prop_as_string(node, "name", &modname);
5841 gretl_xml_get_prop_as_string(node, "depvar", &pmod->depvar);
5842
5843 cur = node->xmlChildrenNode;
5844
5845 #if XDEBUG
5846 fprintf(stderr, "reading model child nodes...\n");
5847 #endif
5848
5849 while (cur != NULL) {
5850 if (!xmlStrcmp(cur->name, (XUC) "sample")) {
5851 gretl_xml_get_prop_as_int(cur, "t1", &pmod->smpl.t1);
5852 gretl_xml_get_prop_as_int(cur, "t2", &pmod->smpl.t2);
5853 gretl_xml_get_prop_as_unsigned_int(cur, "rseed", &pmod->smpl.rseed);
5854 } else if (!xmlStrcmp(cur->name, (XUC) "coeff")) {
5855 pmod->coeff = gretl_xml_get_double_array(cur, doc, &n, err);
5856 } else if (!xmlStrcmp(cur->name, (XUC) "sderr")) {
5857 pmod->sderr = gretl_xml_get_double_array(cur, doc, &n, err);
5858 } else if (!xmlStrcmp(cur->name, (XUC) "uhat")) {
5859 pmod->uhat = gretl_xml_get_double_array(cur, doc, &n, err);
5860 } else if (!xmlStrcmp(cur->name, (XUC) "yhat")) {
5861 pmod->yhat = gretl_xml_get_double_array(cur, doc, &n, err);
5862 } else if (!xmlStrcmp(cur->name, (XUC) "xpx")) {
5863 /* note: allow this one to fail silently */
5864 pmod->xpx = gretl_xml_get_double_array(cur, doc, &n, NULL);
5865 } else if (!xmlStrcmp(cur->name, (XUC) "vcv")) {
5866 pmod->vcv = gretl_xml_get_double_array(cur, doc, &n, err);
5867 } else if (!xmlStrcmp(cur->name, (XUC) "list")) {
5868 pmod->list = gretl_xml_get_list(cur, doc, err);
5869 } else if (!xmlStrcmp(cur->name, (XUC) "tests")) {
5870 *err = attach_model_tests_from_xml(pmod, cur);
5871 } else if (!xmlStrcmp(cur->name, (XUC) "params")) {
5872 *err = attach_model_params_from_xml(cur, doc, pmod);
5873 } else if (!xmlStrcmp(cur->name, (XUC) "submask")) {
5874 fprintf(stderr, "getting submask...\n");
5875 *err = model_submask_from_xml(cur, doc, pmod);
5876 fprintf(stderr, "got it\n");
5877 } else if (!xmlStrcmp(cur->name, (XUC) "missmask")) {
5878 if (!gretl_xml_node_get_string(cur, doc, &pmod->missmask)) {
5879 *err = 1;
5880 }
5881 } else if (!xmlStrcmp(cur->name, (XUC) "arinfo")) {
5882 *err = arinfo_from_xml(cur, doc, pmod);
5883 } else if (!xmlStrcmp(cur->name, (XUC) "data-items")) {
5884 *err = model_data_items_from_xml(cur, doc, pmod, fixopt);
5885 }
5886
5887 if (*err) {
5888 fprintf(stderr, "gretl_model_from_XML: block 3: err = %d reading '%s'\n",
5889 *err, cur->name);
5890 break;
5891 }
5892
5893 cur = cur->next;
5894 }
5895
5896 gretl_pop_c_numeric_locale();
5897
5898 /* backward compatibility */
5899 if (!*err && pmod->nparams > pmod->ncoeff && pmod->depvar == NULL) {
5900 int i, np = pmod->nparams - 1;
5901
5902 pmod->depvar = pmod->params[0];
5903 for (i=0; i<np; i++) {
5904 pmod->params[i] = pmod->params[i+1];
5905 }
5906 pmod->params[np] = NULL;
5907 pmod->nparams = np;
5908 }
5909
5910 if (!*err && pmod->opt != 0) {
5911 /* compat for some options whose flags have been changed */
5912 if (pmod->opt & OPT_E) {
5913 if (pmod->ci == GARCH) {
5914 /* standardized residuals */
5915 pmod->opt &= ~OPT_E;
5916 pmod->opt |= OPT_Z;
5917 }
5918 }
5919 if (pmod->opt & OPT_W) {
5920 if (pmod->ci == PANEL || pmod->ci == IVREG) {
5921 /* weights */
5922 pmod->opt &= ~OPT_W;
5923 pmod->opt |= OPT_H;
5924 } else if (pmod->ci == DURATION) {
5925 /* Weibull */
5926 pmod->opt &= ~OPT_W;
5927 pmod->opt |= OPT_B;
5928 }
5929 }
5930 }
5931
5932 if (!*err && pmod->list != NULL) {
5933 maybe_convert_listsep(pmod->list, dset);
5934 }
5935
5936 if (modname != NULL) {
5937 gretl_model_set_name(pmod, modname);
5938 free(modname);
5939 }
5940
5941 bailout:
5942
5943 if (*err) {
5944 if (pmod != NULL) {
5945 gretl_model_free(pmod);
5946 pmod = NULL;
5947 }
5948 }
5949
5950 return pmod;
5951 }
5952
5953 /**
5954 * gretl_model_copy:
5955 * @pmod: pointer to #MODEL to copy.
5956 *
5957 * Does a deep copy of @pmod: allocates a new #MODEL pointer
5958 * which has its own allocated copies of all the pointer
5959 * members of @pmod. The only feature of @pmod that is
5960 * not duplicated is the reference count, which is set
5961 * to zero in the copy.
5962 *
5963 * Returns: the copied model, or %NULL on failure.
5964 */
5965
gretl_model_copy(MODEL * pmod)5966 MODEL *gretl_model_copy (MODEL *pmod)
5967 {
5968 MODEL *new = malloc(sizeof *new);
5969
5970 #if MDEBUG
5971 fprintf(stderr, "gretl_model_copy: copying %p, allocated at %p\n",
5972 pmod, new);
5973 #endif
5974
5975 if (new != NULL) {
5976 int err = copy_model(new, pmod);
5977
5978 if (err) {
5979 clear_model(new);
5980 free(new);
5981 new = NULL;
5982 } else {
5983 /* nota bene */
5984 new->refcount = 0;
5985 }
5986 }
5987
5988 return new;
5989 }
5990
5991 /**
5992 * swap_models:
5993 * @targ: pointer to target #MODEL.
5994 * @src: pointer to source #MODEL.
5995 *
5996 * Swaps the content of the two model pointers.
5997 */
5998
swap_models(MODEL * targ,MODEL * src)5999 void swap_models (MODEL *targ, MODEL *src)
6000 {
6001 MODEL tmp = *targ;
6002
6003 #if MDEBUG
6004 fprintf(stderr, "swap_models: %p <-> %p\n", targ, src);
6005 #endif
6006
6007 *targ = *src;
6008 *src = tmp;
6009 }
6010
6011 #define normality_test(ci, opt) (ci == MODTEST && (opt & OPT_N))
6012
6013 /**
6014 * command_ok_for_model:
6015 * @test_ci: index of command to be tested.
6016 * @opt: option for command to be tested.
6017 * @pmod: the gretl model in question.
6018 *
6019 * Check to see if the model-related command in question is
6020 * meaningful and acceptable in the context of the estimator
6021 * associated with @pmod. Note, though, that this function may
6022 * give a "false positive": to be quite sure, we may need to
6023 * know more about the model (e.g. specific options used).
6024 * See also model_test_ok().
6025 *
6026 * Returns: 1 if the command seems OK, otherwise 0.
6027 */
6028
command_ok_for_model(int test_ci,gretlopt opt,const MODEL * pmod)6029 int command_ok_for_model (int test_ci, gretlopt opt,
6030 const MODEL *pmod)
6031 {
6032 int between = 0;
6033 int regular_panel = 0;
6034 int mci = pmod->ci;
6035 int ok = 1;
6036
6037 if (mci == MIDASREG) {
6038 /* treat as a case of NLS */
6039 mci = NLS;
6040 }
6041
6042 if (mci == NLS && test_ci == FCAST) {
6043 return 1;
6044 }
6045
6046 if (mci == BIPROBIT) {
6047 /* can we support anything else? */
6048 return test_ci == RESTRICT;
6049 }
6050
6051 if (test_ci == BKW) {
6052 /* most models should be OK */
6053 return pmod->ncoeff > 1 && (pmod->vcv != NULL || pmod->xpx != NULL);
6054 }
6055
6056 if (NONLIST_MODEL(mci)) {
6057 return (test_ci == RESTRICT ||
6058 test_ci == TABPRINT ||
6059 (mci != MLE && normality_test(test_ci, opt)));
6060 }
6061
6062 between = gretl_is_between_model(pmod);
6063 regular_panel = gretl_is_regular_panel_model(pmod);
6064
6065 switch (test_ci) {
6066
6067 case ADD:
6068 if (mci == ARMA || mci == GARCH ||
6069 mci == HECKIT || mci == INTREG) {
6070 ok = 0;
6071 } else if (mci == PANEL && (pmod->opt & OPT_B)) {
6072 ok = 0;
6073 } else if (opt & OPT_L) {
6074 /* --lm variant: OLS only? */
6075 ok = (mci == OLS);
6076 }
6077 break;
6078
6079 case OMIT:
6080 if (mci == ARMA || mci == GARCH || mci == INTREG ||
6081 mci == DPANEL) {
6082 ok = 0;
6083 } else if (between) {
6084 ok = 0;
6085 }
6086 break;
6087
6088 case VIF:
6089 if (mci == IVREG || mci == ARMA || mci == GARCH ||
6090 mci == PANEL || mci == DPANEL) {
6091 ok = 0;
6092 }
6093 break;
6094
6095 case EQNPRINT:
6096 if (mci == ARMA || mci == DPANEL ||
6097 mci == HECKIT || mci == INTREG) {
6098 ok = 0;
6099 }
6100 break;
6101
6102 case MODTEST:
6103 if (opt & OPT_H) {
6104 /* ARCH */
6105 ok = (mci != ARCH && mci != GARCH);
6106 } else if (opt & OPT_C) {
6107 /* common factor restriction */
6108 ok = (mci == AR1);
6109 } else if (opt & OPT_D) {
6110 /* x-sectional dependence */
6111 ok = !between;
6112 } else if (opt & OPT_N) {
6113 /* normality */
6114 if (mci == LOGIT || mci == HECKIT || mci == DURATION) {
6115 /* POISSON, NEGBIN? */
6116 ok = 0;
6117 }
6118 } else if (mci != OLS) {
6119 if (mci == IVREG && (opt & (OPT_A | OPT_W))) {
6120 /* Autocorr. and H'sked. supported for IVREG */
6121 ok = 1;
6122 } else if (mci == ARMA && (opt & OPT_A)) {
6123 ok = 1;
6124 } else if (mci == PANEL && (opt & OPT_P)) {
6125 ok = 1;
6126 } else if (regular_panel && (opt & OPT_A)) {
6127 ok = 1;
6128 } else {
6129 ok = 0;
6130 }
6131 }
6132 break;
6133
6134 case CHOW:
6135 case CUSUM:
6136 case QLRTEST:
6137 case LEVERAGE:
6138 case RESET:
6139 case PANSPEC:
6140 /* OLS-only tests */
6141 ok = (mci == OLS);
6142 break;
6143
6144 case RESTRICT:
6145 if (mci == LAD || mci == QUANTREG) {
6146 ok = 0;
6147 }
6148 break;
6149
6150 default:
6151 break;
6152 }
6153
6154 return ok;
6155 }
6156
6157 /**
6158 * model_test_ok:
6159 * @ci: index of a model test command.
6160 * @opt: option associated with test command, if any.
6161 * @pmod: the model to be tested.
6162 * @dset: dataset information.
6163 *
6164 * A more rigorous version of command_ok_for_model(). Use
6165 * this function if the extra information is available.
6166 *
6167 * Returns: 1 if the test command @ci (with possible option
6168 * @opt) is acceptable in the context of the model @pmod, and
6169 * the dataset described by @dset, otherwise 0.
6170 */
6171
model_test_ok(int ci,gretlopt opt,const MODEL * pmod,const DATASET * dset)6172 int model_test_ok (int ci, gretlopt opt, const MODEL *pmod,
6173 const DATASET *dset)
6174 {
6175 int ok = command_ok_for_model(ci, opt, pmod);
6176
6177 /* for now we'll treat MIDASREG as a case of NLS */
6178 if (ci == MIDASREG) {
6179 ci = NLS;
6180 }
6181
6182 if (ok && pmod->missmask != NULL) {
6183 /* can't do these with embedded missing obs? */
6184 if (gretl_is_regular_panel_model(pmod)) {
6185 ; /* OK? */
6186 } else if (ci == CUSUM || ci == BDS ||
6187 (ci == MODTEST && (opt & (OPT_A | OPT_H)))) {
6188 ok = 0;
6189 }
6190 }
6191
6192 if (ok && pmod->ncoeff == 1) {
6193 if (ci == COEFFSUM) {
6194 ok = 0;
6195 } else if (pmod->ifc && ci == MODTEST) {
6196 /* const only: rule out squares, logs, h'sked */
6197 if (opt & (OPT_W | OPT_B | OPT_S | OPT_L)) {
6198 ok = 0;
6199 }
6200 }
6201 }
6202
6203 if (ci == MODTEST && (opt & OPT_A) &&
6204 gretl_is_regular_panel_model(pmod)) {
6205 return 1;
6206 }
6207
6208 if (ok && !dataset_is_time_series(dset)) {
6209 /* time-series-only tests */
6210 if (ci == CUSUM || ci == QLRTEST || ci == BDS ||
6211 (ci == MODTEST && (opt & (OPT_H | OPT_A)))) {
6212 ok = 0;
6213 }
6214 }
6215
6216 if (ok && !dataset_is_panel(dset)) {
6217 /* panel-only tests */
6218 if (ci == PANSPEC) {
6219 ok = 0;
6220 } else if (ci == MODTEST && (opt & (OPT_P | OPT_D))) {
6221 ok = 0;
6222 }
6223 }
6224
6225 if (ok && pmod->ncoeff - pmod->ifc <= 1 && ci == VIF) {
6226 /* needs at least two independent vars */
6227 ok = 0;
6228 }
6229
6230 if (ok && ci == MODTEST && (opt & OPT_C)) {
6231 /* common factor test */
6232 if (pmod->opt & OPT_P) {
6233 /* ?? check what this means! */
6234 ok = 0;
6235 }
6236 }
6237
6238 if (ok && ci == RESTRICT) {
6239 if (gretl_model_get_int(pmod, "null-model")) {
6240 ok = 0;
6241 }
6242 }
6243
6244 return ok;
6245 }
6246
6247 static int gretl_model_count = 0;
6248
get_model_count(void)6249 int get_model_count (void)
6250 {
6251 return gretl_model_count;
6252 }
6253
set_model_count(int c)6254 void set_model_count (int c)
6255 {
6256 gretl_model_count = c;
6257 }
6258
model_count_plus(void)6259 int model_count_plus (void)
6260 {
6261 return ++gretl_model_count;
6262 }
6263
model_count_minus(MODEL * pmod)6264 void model_count_minus (MODEL *pmod)
6265 {
6266 if (pmod == NULL) {
6267 --gretl_model_count;
6268 } else if (pmod->ID > 0) {
6269 gretl_model_count = pmod->ID - 1;
6270 pmod->ID = 0;
6271 }
6272 }
6273
set_model_id(MODEL * pmod,gretlopt opt)6274 void set_model_id (MODEL *pmod, gretlopt opt)
6275 {
6276 /* always record model estimation time */
6277 pmod->esttime = gretl_monotonic_time();
6278
6279 /* Experimental: limit setting of sequential model
6280 number to "visible" models estimated in main
6281 script or session (not in functions).
6282 */
6283 if (opt & OPT_A) {
6284 /* An auxiliary model? Likely, but OPT_A has
6285 special meaning for a few estimators
6286 */
6287 if (pmod->ci != DPANEL && pmod->ci != GARCH &&
6288 pmod->ci != ARMA) {
6289 return;
6290 }
6291 }
6292
6293 if ((opt & OPT_Q) && !(opt & OPT_W)) {
6294 /* --quiet and not --window, so "invisible" */
6295 return;
6296 } else if (gretl_function_depth() > 0) {
6297 /* model was estimated inside a function */
6298 return;
6299 }
6300
6301 if (pmod->errcode == 0) {
6302 pmod->ID = ++gretl_model_count;
6303 }
6304 }
6305
model_list_to_string(int * list,char * buf)6306 void model_list_to_string (int *list, char *buf)
6307 {
6308 char numstr[8];
6309 int i;
6310
6311 for (i=1; i<=list[0]; i++) {
6312 if (list[i] == LISTSEP) {
6313 strcat(buf, "; ");
6314 } else {
6315 sprintf(numstr, "%d ", list[i]);
6316 strcat(buf, numstr);
6317 }
6318 }
6319 }
6320
highest_numbered_var_in_model(const MODEL * pmod,const DATASET * dset)6321 int highest_numbered_var_in_model (const MODEL *pmod,
6322 const DATASET *dset)
6323 {
6324 int i, v, vmax = 0;
6325 int gotsep = 0;
6326
6327 if (pmod->ci == MLE || pmod->ci == GMM || pmod->list == NULL) {
6328 return 0;
6329 }
6330
6331 for (i=1; i<=pmod->list[0]; i++) {
6332 v = pmod->list[i];
6333 if (v == LISTSEP) {
6334 gotsep = 1;
6335 continue;
6336 }
6337 if (v >= dset->v) {
6338 /* temporary variables, already gone? */
6339 continue;
6340 }
6341 if ((pmod->ci == ARMA || pmod->ci == GARCH) && !gotsep) {
6342 /* real vars start after LISTSEP */
6343 continue;
6344 }
6345 #if 0
6346 fprintf(stderr, "highest numbered... checking var %d\n", v);
6347 #endif
6348 if (v > vmax) {
6349 vmax = v;
6350 }
6351 if (pmod->ci == NLS) {
6352 /* only the dependent var can be tested */
6353 break;
6354 }
6355 }
6356
6357 /* clustered standard errors? */
6358 v = gretl_model_get_cluster_var(pmod);
6359 if (v > vmax) vmax = v;
6360
6361 /* auxiliary variables for some model types */
6362
6363 if (pmod->ci == WLS) {
6364 if (pmod->nwt > vmax) vmax = pmod->nwt;
6365 } else if (pmod->ci == INTREG) {
6366 v = gretl_model_get_int(pmod, "lovar");
6367 if (v > vmax) vmax = v;
6368 v = gretl_model_get_int(pmod, "hivar");
6369 if (v > vmax) vmax = v;
6370 } else if (COUNT_MODEL(pmod->ci)) {
6371 v = gretl_model_get_int(pmod, "offset_var");
6372 if (v > vmax) vmax = v;
6373 } else if (pmod->ci == DURATION) {
6374 v = gretl_model_get_int(pmod, "cens_var");
6375 if (v > vmax) vmax = v;
6376 }
6377
6378 return vmax;
6379 }
6380
mle_criteria(MODEL * pmod,int addk)6381 int mle_criteria (MODEL *pmod, int addk)
6382 {
6383 int err = 0;
6384
6385 if (na(pmod->lnL)) {
6386 pmod->criterion[C_AIC] = NADBL;
6387 pmod->criterion[C_BIC] = NADBL;
6388 pmod->criterion[C_HQC] = NADBL;
6389 err = 1;
6390 } else {
6391 int k = pmod->ncoeff + addk;
6392 int n = pmod->nobs;
6393
6394 pmod->criterion[C_AIC] = -2.0 * pmod->lnL + 2.0 * k;
6395 pmod->criterion[C_BIC] = -2.0 * pmod->lnL + k * log(n);
6396 pmod->criterion[C_HQC] = -2.0 * pmod->lnL + 2 * k * log(log(n));
6397 }
6398
6399 return err;
6400 }
6401
model_use_zscore(const MODEL * pmod)6402 int model_use_zscore (const MODEL *pmod)
6403 {
6404 if (pmod == NULL) {
6405 /* modprint */
6406 return 1;
6407 } else if (gretl_model_get_int(pmod, "dfcorr")) {
6408 /* override ASYMPTOTIC_MODEL if need be */
6409 return 0;
6410 } else if (pmod->ci == OLS && (pmod->opt & OPT_N)) {
6411 /* OLS with --no-df-corr */
6412 return 1;
6413 } else if (ASYMPTOTIC_MODEL(pmod->ci)) {
6414 return 1;
6415 } else if (pmod->ci == PANEL && (pmod->opt & OPT_U)) {
6416 return 1;
6417 } else if ((pmod->opt & OPT_R) && libset_get_bool(ROBUST_Z)) {
6418 return 1;
6419 } else {
6420 return 0;
6421 }
6422 }
6423
coeff_pval(int ci,double x,int df)6424 double coeff_pval (int ci, double x, int df)
6425 {
6426 double p = NADBL;
6427
6428 if (!na(x)) {
6429 if (df == 0 || ASYMPTOTIC_MODEL(ci)) {
6430 p = normal_pvalue_2(x);
6431 } else {
6432 p = student_pvalue_2(df, x);
6433 }
6434 }
6435
6436 return p;
6437 }
6438
model_coeff_pval(const MODEL * pmod,double x)6439 double model_coeff_pval (const MODEL *pmod, double x)
6440 {
6441 double p = NADBL;
6442
6443 if (!na(x)) {
6444 if (model_use_zscore(pmod)) {
6445 p = normal_pvalue_2(x);
6446 } else if (pmod->ci == AR) {
6447 /* backward compat (?) */
6448 int k = pmod->arinfo->arlist[0];
6449 int dfd = pmod->dfd;
6450
6451 if (k > 1) {
6452 dfd += pmod->ncoeff - k;
6453 }
6454 p = student_pvalue_2(dfd, x);
6455 } else {
6456 p = student_pvalue_2(pmod->dfd, x);
6457 }
6458 }
6459
6460 return p;
6461 }
6462
6463 /**
6464 * gretl_model_allocate_param_names:
6465 * @pmod: pointer to target model.
6466 * @k: number of strings to allocate.
6467 *
6468 * Allocate an array of @k strings to hold the names given to
6469 * the associated coefficients, in a model where these strings
6470 * are not simply given by the names of the independent variables.
6471 *
6472 * Returns: 0 on success, non-zero on error.
6473 */
6474
gretl_model_allocate_param_names(MODEL * pmod,int k)6475 int gretl_model_allocate_param_names (MODEL *pmod, int k)
6476 {
6477 pmod->params = strings_array_new_with_length(k, PNAMELEN);
6478
6479 if (pmod->params == NULL) {
6480 pmod->errcode = E_ALLOC;
6481 }
6482
6483 if (!pmod->errcode) {
6484 pmod->nparams = k;
6485 }
6486
6487 return pmod->errcode;
6488 }
6489
6490 /**
6491 * gretl_model_set_param_name:
6492 * @pmod: pointer to target model.
6493 * @i: 0-based index of value to set.
6494 * @name: string value to set.
6495 *
6496 * Returns: 0 on success, non-zero on error.
6497 */
6498
gretl_model_set_param_name(MODEL * pmod,int i,const char * name)6499 int gretl_model_set_param_name (MODEL *pmod, int i, const char *name)
6500 {
6501 if (pmod->params == NULL || i < 0 || i >= pmod->nparams ||
6502 name == NULL) {
6503 return E_DATA;
6504 } else {
6505 pmod->params[i][0] = '\0';
6506 if (strlen(name) >= PNAMELEN) {
6507 strncat(pmod->params[i], name, PNAMELEN - 2);
6508 strcat(pmod->params[i], "~");
6509 } else {
6510 strncat(pmod->params[i], name, PNAMELEN - 1);
6511 }
6512 return 0;
6513 }
6514 }
6515
count_coeffs(int p,const char * pmask,int q,const char * qmask,int P,int Q,int r,int ifc)6516 static int count_coeffs (int p, const char *pmask,
6517 int q, const char *qmask,
6518 int P, int Q, int r, int ifc)
6519 {
6520 int i, n = P + Q + r + ifc;
6521
6522 for (i=0; i<p; i++) {
6523 if (arma_included(pmask, i)) {
6524 n++;
6525 }
6526 }
6527
6528 for (i=0; i<q; i++) {
6529 if (arma_included(qmask, i)) {
6530 n++;
6531 }
6532 }
6533
6534 return n;
6535 }
6536
6537 /**
6538 * gretl_model_add_arma_varnames:
6539 * @pmod: pointer to target model.
6540 * @dset: dataset information.
6541 * @yno: ID number of dependent variable.
6542 * @p: non-seasonal (max) AR order.
6543 * @q: non-seasonal (max) MA order.
6544 * @pmask: for masking out specific AR lags.
6545 * @qmask: for masking out specific MA lags.
6546 * @P: seasonal AR order.
6547 * @Q: seasonal MA order.
6548 * @r: number of exogenous regressors (excluding the constant).
6549 *
6550 * Composes a set of names to be given to the regressors in an
6551 * ARMA model.
6552 *
6553 * Returns: 0 on success, non-zero on error.
6554 */
6555
gretl_model_add_arma_varnames(MODEL * pmod,const DATASET * dset,int yno,int p,int q,const char * pmask,const char * qmask,int P,int Q,int r)6556 int gretl_model_add_arma_varnames (MODEL *pmod, const DATASET *dset,
6557 int yno, int p, int q,
6558 const char *pmask, const char *qmask,
6559 int P, int Q,
6560 int r)
6561 {
6562 int nullmod = 0;
6563 int nc, xstart;
6564 int i, j;
6565
6566 nc = count_coeffs(p, pmask, q, qmask, P, Q, r, pmod->ifc);
6567
6568 if (pmod->depvar != NULL) {
6569 free(pmod->depvar);
6570 }
6571
6572 pmod->depvar = gretl_strdup(dset->varname[yno]);
6573 if (pmod->depvar == NULL) {
6574 pmod->errcode = E_ALLOC;
6575 return 1;
6576 }
6577
6578 if (pmod->nparams > 0 && pmod->params != NULL) {
6579 for (i=0; i<pmod->nparams; i++) {
6580 free(pmod->params[i]);
6581 }
6582 free(pmod->params);
6583 }
6584
6585 nullmod = gretl_model_get_int(pmod, "null-model");
6586
6587 if (nc == 0 && nullmod) {
6588 /* special case of null model */
6589 nc = 1;
6590 }
6591
6592 pmod->params = strings_array_new_with_length(nc, VNAMELEN);
6593 if (pmod->params == NULL) {
6594 free(pmod->depvar);
6595 pmod->depvar = NULL;
6596 pmod->errcode = E_ALLOC;
6597 return 1;
6598 }
6599
6600 pmod->nparams = nc;
6601
6602 if (pmod->ifc || nullmod) {
6603 strcpy(pmod->params[0], dset->varname[0]);
6604 j = 1;
6605 } else {
6606 j = 0;
6607 }
6608
6609 for (i=0; i<p; i++) {
6610 if (arma_included(pmask, i)) {
6611 sprintf(pmod->params[j++], "phi_%d", i + 1);
6612 }
6613 }
6614
6615 for (i=0; i<P; i++) {
6616 sprintf(pmod->params[j++], "Phi_%d", i + 1);
6617 }
6618
6619 for (i=0; i<q; i++) {
6620 if (arma_included(qmask, i)) {
6621 sprintf(pmod->params[j++], "theta_%d", i + 1);
6622 }
6623 }
6624
6625 for (i=0; i<Q; i++) {
6626 sprintf(pmod->params[j++], "Theta_%d", i + 1);
6627 }
6628
6629 xstart = arma_depvar_pos(pmod) + 1;
6630
6631 for (i=0; i<r; i++) {
6632 strcpy(pmod->params[j++], dset->varname[pmod->list[xstart+i]]);
6633 }
6634
6635 return 0;
6636 }
6637
6638 /**
6639 * gretl_model_add_panel_varnames:
6640 * @pmod: pointer to target model.
6641 * @dset: dataset information.
6642 * @ulist: list of index numbers of cross-sectional
6643 * units included in the model.
6644 *
6645 * Composes a set of names to be given to the regressors in an
6646 * panel model.
6647 *
6648 * Returns: 0 on success, non-zero on error.
6649 */
6650
gretl_model_add_panel_varnames(MODEL * pmod,const DATASET * dset,const int * ulist)6651 int gretl_model_add_panel_varnames (MODEL *pmod, const DATASET *dset,
6652 const int *ulist)
6653 {
6654 int np = pmod->ncoeff;
6655 int i, j, v;
6656
6657 pmod->depvar = gretl_strdup(dset->varname[pmod->list[1]]);
6658 if (pmod->depvar == NULL) {
6659 pmod->errcode = E_ALLOC;
6660 return 1;
6661 }
6662
6663 pmod->params = strings_array_new_with_length(np, VNAMELEN);
6664 if (pmod->params == NULL) {
6665 pmod->errcode = E_ALLOC;
6666 return 1;
6667 }
6668
6669 pmod->nparams = np;
6670
6671 j = 1;
6672 for (i=0; i<np; i++) {
6673 v = pmod->list[i+2];
6674 if (v < dset->v) {
6675 strcpy(pmod->params[i], dset->varname[v]);
6676 } else if (ulist != NULL) {
6677 sprintf(pmod->params[i], "ahat_%d", ulist[j++]);
6678 } else {
6679 sprintf(pmod->params[i], "ahat_%d", j++);
6680 }
6681 }
6682
6683 return 0;
6684 }
6685
6686 /**
6687 * gretl_model_add_allocated_varnames:
6688 * @pmod: pointer to target model.
6689 * @vnames: array of names of independent variables.
6690 *
6691 * Attaches an allocated set of variable names to be used
6692 * when printing model results, for use in special cases
6693 * where we can't just reference names from the list of
6694 * regressors attached to the model. The number of strings
6695 * must match the number of coefficients, given by the
6696 * %ncoeff member of @pmod.
6697 *
6698 * Note that @pmod "takes charge" of the array @vnames:
6699 * this will be freed when the model is destroyed.
6700 */
6701
gretl_model_add_allocated_varnames(MODEL * pmod,char ** vnames)6702 void gretl_model_add_allocated_varnames (MODEL *pmod, char **vnames)
6703 {
6704 pmod->nparams = pmod->ncoeff;
6705 pmod->params = vnames;
6706 }
6707
6708 /**
6709 * gretl_model_add_y_median:
6710 * @pmod: pointer to target model.
6711 * @y: array containing the dependent variable.
6712 *
6713 * Calculates the median of @y using the valid observations
6714 * with the model's sample range and attaches the median
6715 * to the model as data under the key %ymedian.
6716 *
6717 * Returns: 0 on success or error code on error.
6718 */
6719
gretl_model_add_y_median(MODEL * pmod,const double * y)6720 int gretl_model_add_y_median (MODEL *pmod, const double *y)
6721 {
6722 int T = pmod->t2 - pmod->t1 + 1;
6723 double *sy, m;
6724 int t, n, ok, n2p;
6725
6726 sy = malloc(T * sizeof *sy);
6727
6728 if (sy == NULL) {
6729 return E_ALLOC;
6730 }
6731
6732 n = 0;
6733 for (t=pmod->t1; t<=pmod->t2; t++) {
6734 if (pmod->uhat != NULL) {
6735 ok = !na(pmod->uhat[t]);
6736 } else {
6737 ok = !model_missing(pmod, t);
6738 }
6739 if (ok) {
6740 sy[n++] = y[t];
6741 }
6742 }
6743
6744 if (n == 0) {
6745 free(sy);
6746 return E_DATA;
6747 }
6748
6749 qsort(sy, n, sizeof *sy, gretl_compare_doubles);
6750
6751 n2p = (T = n / 2) + 1;
6752 m = (n % 2)? sy[n2p - 1] : 0.5 * (sy[T - 1] + sy[n2p - 1]);
6753
6754 gretl_model_set_double(pmod, "ymedian", m);
6755
6756 free(sy);
6757
6758 return 0;
6759 }
6760
gretl_model_add_normality_test(MODEL * pmod,double X2)6761 int gretl_model_add_normality_test (MODEL *pmod, double X2)
6762 {
6763 ModelTest *test = model_test_new(GRETL_TEST_NORMAL);
6764 int err = 0;
6765
6766 if (test != NULL) {
6767 model_test_set_teststat(test, GRETL_STAT_NORMAL_CHISQ);
6768 model_test_set_dfn(test, 2);
6769 model_test_set_value(test, X2);
6770 model_test_set_pvalue(test, chisq_cdf_comp(2, X2));
6771 maybe_add_test_to_model(pmod, test);
6772 } else {
6773 err = E_ALLOC;
6774 }
6775
6776 return err;
6777 }
6778
gretl_model_get_test(MODEL * pmod,ModelTestType ttype)6779 ModelTest *gretl_model_get_test (MODEL *pmod, ModelTestType ttype)
6780 {
6781 int i;
6782
6783 for (i=0; i<pmod->ntests; i++) {
6784 if (pmod->tests[i].type == ttype) {
6785 return &pmod->tests[i];
6786 }
6787 }
6788
6789 return NULL;
6790 }
6791
xneq(double x,double y)6792 static int xneq (double x, double y)
6793 {
6794 double reldiff;
6795
6796 if (x == 0.0) {
6797 reldiff = fabs(y);
6798 } else if (y == 0.0) {
6799 reldiff = fabs(x);
6800 } else if (x > y) {
6801 reldiff = fabs((x - y) / y);
6802 } else {
6803 reldiff = fabs((y - x) / x);
6804 }
6805
6806 return reldiff > 1.5e-12;
6807 }
6808
6809 /* try to tell if an OLS model with two independent variables is
6810 actually a quadratic model (x_2 = x_1^2). */
6811
model_is_quadratic(const MODEL * pmod,const int * xlist,const DATASET * dset)6812 static int model_is_quadratic (const MODEL *pmod, const int *xlist,
6813 const DATASET *dset)
6814 {
6815 const double *x1 = dset->Z[xlist[2]];
6816 const double *x2 = dset->Z[xlist[3]];
6817 int t;
6818
6819 for (t=pmod->t1; t<=pmod->t2; t++) {
6820 if (!na(x1[t]) && xneq(x2[t], x1[t] * x1[t])) {
6821 return 0;
6822 }
6823 }
6824
6825 return 1;
6826 }
6827
6828 /* we'll be a bit conservative here so as not to make
6829 fools of ourselves */
6830
6831 #define fitted_formula_ok(c) (c == OLS || c == WLS || \
6832 c == HSK || c == IVREG || \
6833 c == LAD || c == LOGISTIC)
6834
6835 /**
6836 * gretl_model_get_fitted_formula:
6837 * @pmod: pointer to target model.
6838 * @xvar: ID number of variable that _may_ be "x" in the model.
6839 * @dset: dataset information.
6840 *
6841 * If @pmod is a simple linear, quadratic or logistic model,
6842 * and if @xvar is in fact the "x" variable from the model,
6843 * returns a string representing the formula for generating the
6844 * fitted values as a function of x. This formula may be used
6845 * in the context of a fitted versus actual plot.
6846 *
6847 * Returns: formula for fitted values, or %NULL if this is
6848 * not available.
6849 */
6850
gretl_model_get_fitted_formula(const MODEL * pmod,int xvar,const DATASET * dset)6851 char *gretl_model_get_fitted_formula (const MODEL *pmod, int xvar,
6852 const DATASET *dset)
6853 {
6854 const DATASET *mset;
6855 int *xlist = NULL;
6856 char *ret = NULL;
6857
6858 if (xvar == 0 || pmod->ncoeff > 3 || !fitted_formula_ok(pmod->ci)) {
6859 return NULL;
6860 }
6861
6862 xlist = gretl_model_get_x_list(pmod);
6863 if (xlist == NULL) {
6864 return NULL;
6865 }
6866
6867 if (pmod->dataset != NULL) {
6868 mset = pmod->dataset;
6869 } else {
6870 mset = dset;
6871 }
6872
6873 gretl_push_c_numeric_locale();
6874
6875 if (pmod->ci == LOGISTIC) {
6876 if (pmod->ifc && pmod->ncoeff == 2 && xvar == pmod->list[3]) {
6877 double lmax = gretl_model_get_double(pmod, "lmax");
6878
6879 if (!na(lmax)) {
6880 ret = gretl_strdup_printf("yformula: %g/(1.0+exp(-(%g%s%g*x)))",
6881 lmax, pmod->coeff[0],
6882 (pmod->coeff[1] >= 0)? "+" : "",
6883 pmod->coeff[1]);
6884 }
6885 }
6886 } else if (!pmod->ifc && pmod->ncoeff == 1 && xvar == xlist[1]) {
6887 ret = gretl_strdup_printf("yformula: %g*x", pmod->coeff[0]);
6888 } else if (pmod->ifc && pmod->ncoeff == 2 && xvar == xlist[2]) {
6889 ret = gretl_strdup_printf("yformula: %g%s%g*x", pmod->coeff[0],
6890 (pmod->coeff[1] >= 0)? "+" : "",
6891 pmod->coeff[1]);
6892 } else if (pmod->ifc && pmod->ncoeff == 3 && xvar == xlist[2]) {
6893 if (model_is_quadratic(pmod, xlist, mset)) {
6894 ret = gretl_strdup_printf("yformula: %g%s%g*x%s%g*x**2", pmod->coeff[0],
6895 (pmod->coeff[1] >= 0)? "+" : "",
6896 pmod->coeff[1],
6897 (pmod->coeff[2] >= 0)? "+" : "",
6898 pmod->coeff[2]);
6899 }
6900 }
6901
6902 gretl_pop_c_numeric_locale();
6903
6904 free(xlist);
6905
6906 return ret;
6907 }
6908
6909 /**
6910 * gretl_model_set_name:
6911 * @pmod: pointer to target model.
6912 * @name: the name to give the model.
6913 *
6914 * Sets the name of the given model; this is used in
6915 * printing the model and in displaying it in the
6916 * gretl GUI. Note that a model's name must be no more
6917 * than #MAXSAVENAME bytes in length, including the
6918 * terminating NUL byte; @name is truncated if it is
6919 * too long.
6920 */
6921
gretl_model_set_name(MODEL * pmod,const char * name)6922 void gretl_model_set_name (MODEL *pmod, const char *name)
6923 {
6924 if (name == pmod->name) {
6925 return;
6926 }
6927
6928 if (pmod->name == NULL) {
6929 pmod->name = calloc(1, MAXSAVENAME);
6930 }
6931
6932 if (pmod->name != NULL) {
6933 *pmod->name = '\0';
6934 if (strlen(name) >= MAXSAVENAME) {
6935 char *tmp = g_strdup(name);
6936
6937 strcpy(pmod->name, gretl_utf8_truncate(tmp, MAXSAVENAME-1));
6938 g_free(tmp);
6939 } else {
6940 strcpy(pmod->name, name);
6941 }
6942 }
6943 }
6944
6945 /**
6946 * gretl_model_get_name:
6947 * @pmod: pointer to gretl model.
6948 *
6949 * Returns: the name that has been set for @pmod, if any.
6950 * Note that the value returned may be %NULL.
6951 */
6952
gretl_model_get_name(const MODEL * pmod)6953 const char *gretl_model_get_name (const MODEL *pmod)
6954 {
6955 if (pmod != NULL) {
6956 return pmod->name;
6957 }
6958
6959 return NULL;
6960 }
6961
6962 /**
6963 * gretl_model_get_scalar:
6964 * @pmod: pointer to target model.
6965 * @idx: index for the scalar value that is wanted.
6966 * @dset: dataset struct.
6967 * @err: location to receive error code (required).
6968 *
6969 * Retrieves a specified scalar statistic from @pmod:
6970 * @idx must be less than %M_SCALAR_MAX.
6971 *
6972 * Returns: the requested statistic, or #NADBL on failure,
6973 * in which case @err will contain a non-zero error code.
6974 */
6975
gretl_model_get_scalar(MODEL * pmod,ModelDataIndex idx,DATASET * dset,int * err)6976 double gretl_model_get_scalar (MODEL *pmod, ModelDataIndex idx,
6977 DATASET *dset, int *err)
6978 {
6979 double x = NADBL;
6980
6981 if (pmod == NULL) {
6982 *err = E_BADSTAT;
6983 return x;
6984 }
6985
6986 if (idx == M_GMMCRIT && pmod->ci != GMM) {
6987 *err = E_BADSTAT;
6988 return x;
6989 }
6990
6991 switch (idx) {
6992 case M_ESS:
6993 case M_GMMCRIT:
6994 x = pmod->ess;
6995 break;
6996 case M_RSQ:
6997 x = pmod->rsq;
6998 break;
6999 case M_LNL:
7000 x = pmod->lnL;
7001 break;
7002 case M_AIC:
7003 x = pmod->criterion[C_AIC];
7004 break;
7005 case M_BIC:
7006 x = pmod->criterion[C_BIC];
7007 break;
7008 case M_HQC:
7009 x = pmod->criterion[C_HQC];
7010 break;
7011 case M_SIGMA:
7012 x = pmod->sigma;
7013 break;
7014 case M_TRSQ:
7015 if (!na(pmod->rsq)) {
7016 x = pmod->nobs * pmod->rsq;
7017 }
7018 break;
7019 case M_DF:
7020 x = (double) pmod->dfd;
7021 break;
7022 case M_NCOEFF:
7023 x = (double) pmod->ncoeff; /* is this always available? */
7024 break;
7025 case M_T:
7026 x = (double) pmod->nobs;
7027 break;
7028 case M_DW:
7029 x = pmod->dw;
7030 break;
7031 case M_DWPVAL:
7032 x = get_DW_pvalue_for_model(pmod, dset, err);
7033 break;
7034 case M_FSTT:
7035 x = pmod->fstt;
7036 break;
7037 case M_CHISQ:
7038 x = pmod->chisq;
7039 break;
7040 default:
7041 break;
7042 }
7043
7044 if (idx != M_DWPVAL && na(x)) {
7045 *err = E_BADSTAT;
7046 }
7047
7048 return x;
7049 }
7050
7051 /**
7052 * gretl_model_get_series:
7053 * @x: series to fill (must be of length dset->n).
7054 * @pmod: pointer to target model.
7055 * @dset: dataset information.
7056 * @idx: index for the series that is wanted.
7057 *
7058 * Retrieves from @pmod a copy of a specified series (for
7059 * example, regression residuals); @idx must be greater than
7060 * %M_SCALAR_MAX and less than %M_SERIES_MAX.
7061 *
7062 * Returns: 0 on success, non-zero code on error.
7063 */
7064
gretl_model_get_series(double * x,MODEL * pmod,const DATASET * dset,ModelDataIndex idx)7065 int gretl_model_get_series (double *x, MODEL *pmod,
7066 const DATASET *dset,
7067 ModelDataIndex idx)
7068 {
7069 const double *src = NULL;
7070 int t;
7071
7072 if (pmod->t2 - pmod->t1 + 1 > dset->n ||
7073 model_sample_problem(pmod, dset)) {
7074 gretl_errmsg_set(
7075 (idx == M_UHAT)?
7076 _("Can't retrieve uhat: data set has changed") :
7077 (idx == M_YHAT)?
7078 _("Can't retrieve yhat: data set has changed") :
7079 (idx == M_H)?
7080 _("Can't retrieve ht: data set has changed") :
7081 _("Can't retrieve series: data set has changed"));
7082 return E_BADSTAT;
7083 }
7084
7085 if (pmod->ci == BIPROBIT && (idx == M_UHAT || idx == M_YHAT)) {
7086 return E_BADSTAT;
7087 }
7088
7089 if (idx == M_UHAT) {
7090 src = pmod->uhat;
7091 } else if (idx == M_YHAT) {
7092 src = pmod->yhat;
7093 } else if (idx == M_LLT) {
7094 src = gretl_model_get_data(pmod, "llt");
7095 } else if (idx == M_AHAT) {
7096 src = gretl_model_get_data(pmod, "ahat");
7097 } else if (idx == M_H) {
7098 src = gretl_model_get_data(pmod, "garch_h");
7099 }
7100
7101 if (src == NULL && idx != M_SAMPLE) {
7102 return E_BADSTAT;
7103 }
7104
7105 /* allow for internal "just testing" usage */
7106 if (x == NULL) {
7107 return 0;
7108 }
7109
7110 if (idx == M_SAMPLE) {
7111 for (t=0; t<dset->n; t++) {
7112 if (t < pmod->t1 || t > pmod->t2) {
7113 x[t] = 0.0;
7114 } else {
7115 x[t] = model_missing(pmod, t)? 0 : 1;
7116 }
7117 }
7118 } else {
7119 for (t=0; t<dset->n; t++) {
7120 if (t < pmod->t1 || t > pmod->t2) {
7121 x[t] = NADBL;
7122 } else {
7123 x[t] = src[t];
7124 }
7125 }
7126 }
7127
7128 return 0;
7129 }
7130
7131 static gretl_matrix *
model_get_estvec(const MODEL * pmod,int idx,int * err)7132 model_get_estvec (const MODEL *pmod, int idx, int *err)
7133 {
7134 const double *src;
7135 gretl_matrix *v = NULL;
7136 int i;
7137
7138 if (gretl_model_get_data(pmod, "rq_tauvec")) {
7139 *err = E_BADSTAT;
7140 return NULL;
7141 }
7142
7143 src = (idx == M_COEFF)? pmod->coeff : pmod->sderr;
7144
7145 if (src == NULL) {
7146 *err = E_BADSTAT;
7147 return NULL;
7148 }
7149
7150 v = gretl_column_vector_alloc(pmod->ncoeff);
7151 if (v == NULL) {
7152 *err = E_ALLOC;
7153 } else {
7154 for (i=0; i<pmod->ncoeff; i++) {
7155 gretl_vector_set(v, i, src[i]);
7156 }
7157 }
7158
7159 return v;
7160 }
7161
7162 static gretl_matrix *
model_get_hatvec(const MODEL * pmod,int idx,int * err)7163 model_get_hatvec (const MODEL *pmod, int idx, int *err)
7164 {
7165 gretl_matrix *v = NULL;
7166 const double *src = NULL;
7167 int t;
7168
7169 if (idx == M_UHAT) {
7170 src = pmod->uhat;
7171 } else if (idx == M_YHAT) {
7172 src = pmod->yhat;
7173 } else if (idx == M_LLT) {
7174 src = gretl_model_get_data(pmod, "llt");
7175 }
7176
7177 if (src == NULL) {
7178 *err = E_BADSTAT;
7179 return NULL;
7180 }
7181
7182 /* FIXME: we reject NAs when creating a gretl_matrix -- but should
7183 we just substitute NaNs?
7184 */
7185 for (t=pmod->t1; t<=pmod->t2; t++) {
7186 if (na(src[t])) {
7187 *err = E_MISSDATA;
7188 break;
7189 }
7190 }
7191
7192 if (!*err) {
7193 v = gretl_column_vector_alloc(pmod->t2 - pmod->t1 + 1);
7194 if (v == NULL) {
7195 *err = E_ALLOC;
7196 } else {
7197 for (t=pmod->t1; t<=pmod->t2; t++) {
7198 gretl_vector_set(v, t - pmod->t1, src[t]);
7199 }
7200 }
7201 }
7202
7203 return v;
7204 }
7205
7206 static gretl_matrix *
model_get_special_vec(const MODEL * pmod,ModelDataIndex idx,int * err)7207 model_get_special_vec (const MODEL *pmod, ModelDataIndex idx, int *err)
7208 {
7209 gretl_matrix *v = NULL;
7210 double *mdata = NULL;
7211 int t;
7212
7213 if (idx == M_AHAT) {
7214 mdata = gretl_model_get_data(pmod, "ahat");
7215 } else if (idx == M_H) {
7216 mdata = gretl_model_get_data(pmod, "garch_h");
7217 }
7218
7219 if (mdata == NULL) {
7220 fprintf(stderr, "model_get_special_vec: mdata is NULL\n");
7221 *err = E_BADSTAT;
7222 }
7223
7224 v = gretl_column_vector_alloc(pmod->t2 - pmod->t1 + 1);
7225
7226 if (v == NULL) {
7227 *err = E_ALLOC;
7228 } else {
7229 gretl_matrix_set_t1(v, pmod->t1);
7230 gretl_matrix_set_t2(v, pmod->t2);
7231 for (t=pmod->t1; t<=pmod->t2; t++) {
7232 /* FIXME: is this indexation always right? */
7233 gretl_vector_set(v, t - pmod->t1, mdata[t]);
7234 }
7235 }
7236
7237 return v;
7238 }
7239
model_get_rhovec(const MODEL * pmod,int * err)7240 static gretl_matrix *model_get_rhovec (const MODEL *pmod, int *err)
7241 {
7242 gretl_matrix *r = NULL;
7243
7244 if (pmod->ci == AR1) {
7245 double x = gretl_model_get_double(pmod, "rho_gls");
7246
7247 r = gretl_matrix_from_scalar(x);
7248 } else if (pmod->ci != AR) {
7249 r = gretl_matrix_from_scalar(pmod->rho);
7250 } else if (pmod->arinfo == NULL ||
7251 pmod->arinfo->arlist == NULL ||
7252 pmod->arinfo->rho == NULL) {
7253 *err = E_INVARG;
7254 } else {
7255 int i, l0 = pmod->arinfo->arlist[0];
7256 int lmax = pmod->arinfo->arlist[l0];
7257
7258 r = gretl_vector_alloc(lmax);
7259 if (r != NULL) {
7260 gretl_matrix_zero(r);
7261 for (i=1; i<=l0; i++) {
7262 gretl_vector_set(r, pmod->arinfo->arlist[i] - 1,
7263 pmod->arinfo->rho[i-1]);
7264 }
7265 }
7266 }
7267
7268 return r;
7269 }
7270
7271 static gretl_matrix *
model_get_intervals_matrix(const MODEL * pmod,int * err)7272 model_get_intervals_matrix (const MODEL *pmod, int *err)
7273 {
7274 gretl_matrix *m, *ci;
7275
7276 ci = gretl_model_get_data(pmod, "coeff_intervals");
7277 if (ci == NULL) {
7278 *err = E_BADSTAT;
7279 return NULL;
7280 }
7281
7282 m = gretl_matrix_copy(ci);
7283
7284 if (m == NULL) {
7285 *err = E_ALLOC;
7286 }
7287
7288 return m;
7289 }
7290
model_get_special_test(const MODEL * pmod,int type,int * err)7291 static gretl_matrix *model_get_special_test (const MODEL *pmod,
7292 int type, int *err)
7293 {
7294 gretl_matrix *r = NULL;
7295 int i, found = 0;
7296
7297 for (i=0; i<pmod->ntests; i++) {
7298 found = (pmod->tests[i].type == type);
7299 if (found) {
7300 r = gretl_vector_alloc(3);
7301 if (r == NULL) {
7302 *err = E_ALLOC;
7303 return NULL;
7304 }
7305 r->val[0] = pmod->tests[i].value;
7306 r->val[1] = pmod->tests[i].dfn;
7307 r->val[2] = pmod->tests[i].pvalue;
7308 break;
7309 }
7310 }
7311
7312 if (!found) {
7313 *err = E_BADSTAT;
7314 }
7315
7316 return r;
7317 }
7318
7319 /**
7320 * gretl_model_get_matrix:
7321 * @pmod: pointer to target model.
7322 * @idx: index for the matrix that is wanted.
7323 * @err: location to receive error code (required).
7324 *
7325 * Retrieves from @pmod a copy of a specified matrix (for
7326 * example, regression residuals); @idx must be greater than
7327 * %M_SERIES_MAX and less than %M_MATRIX_MAX.
7328 *
7329 * Returns: the allocated matrix, or %NULL on failure,
7330 * in which case @err will contain a non-zero error code.
7331 */
7332
gretl_model_get_matrix(MODEL * pmod,ModelDataIndex idx,int * err)7333 gretl_matrix *gretl_model_get_matrix (MODEL *pmod, ModelDataIndex idx,
7334 int *err)
7335 {
7336 gretl_matrix *M = NULL;
7337
7338 if (pmod == NULL) {
7339 *err = E_BADSTAT;
7340 return NULL;
7341 }
7342
7343 if ((idx == M_UHAT || idx == M_YHAT) &&
7344 (pmod->ci == BIPROBIT || gretl_is_between_model(pmod))) {
7345 /* special: uhat and yhat are matrices */
7346 const char *utag = pmod->ci == BIPROBIT ? "bp_uhat" : "uhat";
7347 const char *ytag = pmod->ci == BIPROBIT ? "bp_yhat" : "yhat";
7348 gretl_matrix *P;
7349
7350 if (idx == M_UHAT) {
7351 P = gretl_model_get_data(pmod, utag);
7352 } else {
7353 P = gretl_model_get_data(pmod, ytag);
7354 }
7355
7356 if (P == NULL) {
7357 *err = E_BADSTAT;
7358 } else {
7359 M = gretl_matrix_copy(P);
7360 if (M == NULL) {
7361 *err = E_ALLOC;
7362 }
7363 }
7364 return M;
7365 }
7366
7367 switch (idx) {
7368 case M_UHAT:
7369 case M_YHAT:
7370 case M_LLT:
7371 M = model_get_hatvec(pmod, idx, err);
7372 break;
7373 case M_COEFF:
7374 case M_SE:
7375 M = model_get_estvec(pmod, idx, err);
7376 break;
7377 case M_COEFF_CI:
7378 M = model_get_intervals_matrix(pmod, err);
7379 break;
7380 case M_VCV:
7381 M = gretl_vcv_matrix_from_model(pmod, NULL, err);
7382 break;
7383 case M_AHAT:
7384 if (gretl_model_get_data(pmod, "ahat") == NULL) {
7385 *err = E_BADSTAT;
7386 } else {
7387 M = model_get_special_vec(pmod, M_AHAT, err);
7388 }
7389 break;
7390 case M_H:
7391 case M_SIGMA:
7392 if (gretl_model_get_data(pmod, "garch_h") == NULL) {
7393 *err = E_BADSTAT;
7394 } else {
7395 M = model_get_special_vec(pmod, M_H, err);
7396 }
7397 break;
7398 case M_RHO:
7399 M = model_get_rhovec(pmod, err);
7400 break;
7401 case M_HAUSMAN:
7402 if (pmod->ci == IVREG) {
7403 M = model_get_special_test(pmod, GRETL_TEST_IV_HAUSMAN, err);
7404 } else if (pmod->ci == PANEL) {
7405 M = model_get_special_test(pmod, GRETL_TEST_PANEL_HAUSMAN, err);
7406 } else {
7407 *err = E_BADSTAT;
7408 }
7409 break;
7410 case M_SARGAN:
7411 M = model_get_special_test(pmod, GRETL_TEST_SARGAN, err);
7412 break;
7413 case M_EHAT:
7414 M = gretl_model_get_data(pmod, "ehat");
7415 if (M == NULL) {
7416 *err = E_BADSTAT;
7417 } else {
7418 M = gretl_matrix_copy(M);
7419 if (M == NULL) {
7420 *err = E_ALLOC;
7421 }
7422 }
7423 break;
7424 default:
7425 *err = E_BADSTAT;
7426 break;
7427 }
7428
7429 if (M == NULL && !*err) {
7430 *err = E_ALLOC;
7431 }
7432
7433 return M;
7434 }
7435
get_vcv_element(MODEL * pmod,const char * s,const DATASET * dset)7436 static double get_vcv_element (MODEL *pmod, const char *s,
7437 const DATASET *dset)
7438 {
7439 char v1str[VNAMELEN], v2str[VNAMELEN];
7440 char fmt[16];
7441 int v1 = 0, v2 = 0;
7442 int i, j, k, gotit;
7443 double ret = NADBL;
7444
7445 if (pmod == NULL || pmod->list == NULL) {
7446 return NADBL;
7447 }
7448
7449 sprintf(fmt, "%%%d[^,],%%%ds", VNAMELEN-1, VNAMELEN-1);
7450
7451 if (sscanf(s, fmt, v1str, v2str) != 2) {
7452 return NADBL;
7453 }
7454
7455 v1 = gretl_model_get_param_number(pmod, dset, v1str);
7456 v2 = gretl_model_get_param_number(pmod, dset, v2str);
7457
7458 if (v1 < 0 || v2 < 0) {
7459 return NADBL;
7460 }
7461
7462 /* make model vcv matrix if need be */
7463 if (pmod->vcv == NULL && makevcv(pmod, pmod->sigma)) {
7464 return NADBL;
7465 }
7466
7467 /* now find the right entry */
7468 if (v1 > v2) {
7469 k = v1;
7470 v1 = v2;
7471 v2 = k;
7472 }
7473
7474 gotit = 0;
7475 k = 0;
7476 for (i=0; i<pmod->ncoeff && !gotit; i++) {
7477 for (j=0; j<pmod->ncoeff; j++) {
7478 if (j < i) {
7479 continue;
7480 }
7481 if (i == v1 && j == v2) {
7482 ret = pmod->vcv[k];
7483 gotit = 1;
7484 break;
7485 }
7486 k++;
7487 }
7488 }
7489
7490 return ret;
7491 }
7492
7493 /* retrieve a specific element from one of the arrays of data
7494 on a model */
7495
7496 double
gretl_model_get_data_element(MODEL * pmod,int idx,const char * s,const DATASET * dset,int * err)7497 gretl_model_get_data_element (MODEL *pmod, int idx, const char *s,
7498 const DATASET *dset, int *err)
7499 {
7500 GretlObjType type;
7501 double x = NADBL;
7502 int vi = 0;
7503
7504 if (pmod == NULL) {
7505 pmod = get_genr_model(&type);
7506 if (pmod == NULL || type != GRETL_OBJ_EQN) {
7507 pmod = get_last_model(&type);
7508 if (pmod == NULL || type != GRETL_OBJ_EQN) {
7509 *err = E_INVARG;
7510 return x;
7511 }
7512 }
7513 }
7514
7515 if (idx == M_RHO) {
7516 int k = atoi(s);
7517
7518 if (k == 1 && pmod->ci == AR1) {
7519 x = gretl_model_get_double(pmod, "rho_gls");
7520 } else if (k == 1 && pmod->ci != AR) {
7521 x = pmod->rho;
7522 } else if (pmod->arinfo == NULL ||
7523 pmod->arinfo->arlist == NULL ||
7524 pmod->arinfo->rho == NULL) {
7525 *err = E_INVARG;
7526 } else {
7527 vi = in_gretl_list(pmod->arinfo->arlist, k);
7528 if (vi > 0) {
7529 x = pmod->arinfo->rho[vi-1];
7530 } else {
7531 *err = E_INVARG;
7532 }
7533 }
7534 } else if (idx == M_VCV) {
7535 x = get_vcv_element(pmod, s, dset);
7536 if (na(x)) {
7537 *err = E_INVARG;
7538 }
7539 } else if (idx == M_COEFF || idx == M_SE) {
7540 vi = gretl_model_get_param_number(pmod, dset, s);
7541 if (vi < 0) {
7542 *err = E_INVARG;
7543 } else {
7544 if (idx == M_COEFF && pmod->coeff != NULL) {
7545 x = pmod->coeff[vi];
7546 } else if (idx == M_SE && pmod->sderr != NULL) {
7547 x = pmod->sderr[vi];
7548 } else {
7549 *err = E_INVARG;
7550 }
7551 }
7552 }
7553
7554 return x;
7555 }
7556
exact_fit_check(const MODEL * pmod,PRN * prn)7557 int exact_fit_check (const MODEL *pmod, PRN *prn)
7558 {
7559 if (pmod->rsq == 1.0) {
7560 pputs(prn, _("The model exhibits an exact linear fit"));
7561 pputc(prn, '\n');
7562 return 1;
7563 }
7564
7565 return 0;
7566 }
7567
maybe_suppress_time_dummies(MODEL * pmod,int ndum)7568 void maybe_suppress_time_dummies (MODEL *pmod, int ndum)
7569 {
7570 const char *s = get_optval_string(pmod->ci, OPT_D);
7571
7572 if (s != NULL && !strcmp(s, "noprint")) {
7573 gretl_model_set_int(pmod, "skipdums", ndum);
7574 }
7575 }
7576