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 #include "../cephes/mconf.h"
20
21 #include "libgretl.h"
22 #include "version.h"
23 #include "bhhh_max.h"
24 #include "libset.h"
25 #include "arma_priv.h"
26
27 #ifdef WIN32
28 # include "gretl_win32.h"
29 # include <io.h>
30 #endif
31
32 #ifndef WIN32
33
tramo_x12a_spawn(const char * workdir,const char * fmt,...)34 static int tramo_x12a_spawn (const char *workdir, const char *fmt, ...)
35 {
36 va_list ap;
37 int i, nargs;
38 int ok;
39 int status = 0, ret = 0;
40 GError *error = NULL;
41 gchar **argv = NULL;
42 gchar *sout = NULL, *serr = NULL;
43 char *s;
44
45 argv = malloc(2 * sizeof *argv);
46 if (argv == NULL) return 1;
47 argv[0] = g_strdup(fmt);
48 argv[1] = NULL;
49 i = nargs = 1;
50
51 va_start(ap, fmt);
52 while ((s = va_arg(ap, char *))) {
53 i++;
54 argv = realloc(argv, (i+1) * sizeof *argv);
55 if (argv == NULL) {
56 status = 1;
57 break;
58 }
59 argv[i-1] = g_strdup(s);
60 argv[i] = NULL;
61 }
62 va_end(ap);
63
64 if (status == 1) return 1;
65
66 nargs = i;
67
68 ok = g_spawn_sync(workdir,
69 argv,
70 NULL,
71 G_SPAWN_SEARCH_PATH,
72 NULL,
73 NULL,
74 &sout,
75 &serr,
76 &status,
77 &error);
78
79 if (!ok) {
80 fprintf(stderr, "spawn: '%s'\n", error->message);
81 g_error_free(error);
82 ret = 1;
83 } else if (serr && *serr) {
84 fprintf(stderr, "stderr: '%s'\n", serr);
85 ret = 1;
86 } else if (status != 0) {
87 fprintf(stderr, "status=%d: stdout: '%s'\n", status, sout);
88 ret = 1;
89 }
90
91 if (serr != NULL) g_free(serr);
92 if (sout != NULL) g_free(sout);
93
94 if (ret != 0) fputc(' ', stderr);
95 for (i=0; i<nargs; i++) {
96 if (ret != 0) fprintf(stderr, "%s ", argv[i]);
97 free(argv[i]);
98 }
99 free(argv);
100 if (ret != 0) fputc('\n', stderr);
101
102 return ret;
103 }
104
105 #endif /* !WIN32 */
106
107 #define X12A_CODE
108 #include "arma_common.c"
109 #undef X12A_CODE
110
add_unique_output_file(MODEL * pmod,const char * path)111 static int add_unique_output_file (MODEL *pmod, const char *path)
112 {
113 char outname[FILENAME_MAX];
114 char unique[FILENAME_MAX];
115 char line[256];
116 FILE *f0, *f1;
117 int err;
118
119 sprintf(outname, "%s.out", path);
120 f0 = gretl_fopen(outname, "r");
121 if (f0 == NULL) {
122 return E_FOPEN;
123 }
124
125 err = gretl_path_compose(unique, FILENAME_MAX, outname, ".XXXXXX");
126 if (err) {
127 fclose(f0);
128 return err;
129 }
130
131 f1 = gretl_mktemp(unique, "w");
132 if (f1 == NULL) {
133 fclose(f0);
134 return E_FOPEN;
135 }
136
137 while (fgets(line, sizeof line, f0)) {
138 fputs(line, f1);
139 }
140
141 fclose(f0);
142 fclose(f1);
143 gretl_remove(outname);
144
145 gretl_model_set_data(pmod, "x12a_output", g_strdup(unique),
146 GRETL_TYPE_STRING, strlen(unique) + 1);
147
148 return 0;
149 }
150
print_iterations(const char * path,PRN * prn)151 static int print_iterations (const char *path, PRN *prn)
152 {
153 char fname[MAXLEN];
154 FILE *fp;
155 char line[129];
156 int print = 0;
157 int err;
158
159 err = gretl_path_compose(fname, MAXLEN, path, ".out");
160 if (err) {
161 return err;
162 }
163
164 fp = gretl_fopen(fname, "r");
165 if (fp == NULL) {
166 fprintf(stderr, "Couldn't read from '%s'\n", fname);
167 err = 1;
168 } else {
169 while (fgets(line, sizeof line, fp)) {
170 if (!strncmp(line, " MODEL EST", 10)) print = 1;
171 if (print) pputs(prn, line);
172 if (!strncmp(line, " Estimatio", 10)) break;
173 }
174 fclose(fp);
175 }
176
177 return err;
178 }
179
180 #define non_yearly_frequency(f) (f != 1 && f != 4 && f != 12)
181
x12_date_to_n(const char * s,const DATASET * dset)182 static int x12_date_to_n (const char *s, const DATASET *dset)
183 {
184 char date[12] = {0};
185 int len = strlen(s);
186
187 if (non_yearly_frequency(dset->pd)) {
188 int t, maj = 0, min = 0;
189 char fmt[16];
190
191 sprintf(fmt, "%%%dd%%2d", len - 2);
192 if (sscanf(s, fmt, &maj, &min) == 2) {
193 t = (maj - 1) * dset->pd + min - 1;
194 } else {
195 t = -1;
196 }
197
198 return t;
199 }
200
201 if (dset->pd > 1) {
202 if (len <= 4) {
203 gchar *tmp = g_strndup(s, len - 2);
204
205 sprintf(date, "%s:%s", tmp, s + len - 2);
206 g_free(tmp);
207 } else {
208 strncat(date, s, 4);
209 strcat(date, ":");
210 strncat(date, s + 4, 4);
211 }
212 } else {
213 strncat(date, s, 4);
214 }
215
216 return dateton(date, dset);
217 }
218
219 /* Parse the statistics from the X12ARIMA output file foo.lks */
220
get_ll_stats(const char * fname,MODEL * pmod)221 static int get_ll_stats (const char *fname, MODEL *pmod)
222 {
223 FILE *fp;
224 char line[80], statname[12];
225 int nefobs = 0;
226 double x;
227
228 fp = gretl_fopen(fname, "r");
229 if (fp == NULL) {
230 fprintf(stderr, "Couldn't read from '%s'\n", fname);
231 return E_FOPEN;
232 }
233
234 pmod->sigma = NADBL;
235
236 gretl_push_c_numeric_locale();
237
238 while (fgets(line, sizeof line, fp)) {
239 if (sscanf(line, "%11s %lf", statname, &x) == 2) {
240 if (!strcmp(statname, "nefobs")) nefobs = (int) x;
241 else if (!strcmp(statname, "var")) pmod->sigma = sqrt(x);
242 else if (!strcmp(statname, "lnlkhd")) pmod->lnL = x;
243 else if (!strcmp(statname, "aic")) pmod->criterion[C_AIC] = x;
244 else if (!strcmp(statname, "bic")) pmod->criterion[C_BIC] = x;
245 else if (!strcmp(statname, "hnquin")) pmod->criterion[C_HQC] = x;
246 }
247 }
248
249 gretl_pop_c_numeric_locale();
250
251 fclose(fp);
252
253 if (nefobs > 0) {
254 pmod->nobs = nefobs;
255 pmod->t1 = pmod->t2 - nefobs + 1;
256 }
257
258 return 0;
259 }
260
261 /* Parse the roots information from the X12ARIMA output file foo.rts */
262
get_roots(const char * fname,MODEL * pmod,arma_info * ainfo)263 static int get_roots (const char *fname, MODEL *pmod,
264 arma_info *ainfo)
265 {
266 FILE *fp;
267 char line[132];
268 int i, nr, err = 0;
269 cmplx *roots;
270
271 nr = ainfo->p + ainfo->q + ainfo->P + ainfo->Q;
272 if (nr == 0) {
273 return 0;
274 }
275
276 roots = malloc(nr * sizeof *roots);
277 if (roots == NULL) {
278 return E_ALLOC;
279 }
280
281 fp = gretl_fopen(fname, "r");
282 if (fp == NULL) {
283 fprintf(stderr, "Couldn't read from '%s'\n", fname);
284 free(roots);
285 return E_FOPEN;
286 }
287
288 gretl_push_c_numeric_locale();
289
290 i = 0;
291 while (fgets(line, sizeof line, fp) && i < nr) {
292 double re, im;
293
294 if (!strncmp(line, "AR", 2) || !strncmp(line, "MA", 2)) {
295 if (sscanf(line, "%*s %*s %*s %lf %lf", &re, &im) == 2) {
296 roots[i].r = re;
297 roots[i].i = im;
298 i++;
299 }
300 }
301 }
302
303 gretl_pop_c_numeric_locale();
304
305 fclose(fp);
306
307 if (i != nr) {
308 fprintf(stderr, "Error reading '%s'\n", fname);
309 free(roots);
310 roots = NULL;
311 err = E_DATA;
312 }
313
314 if (roots != NULL) {
315 gretl_model_set_data(pmod, "roots", roots, GRETL_TYPE_CMPLX_ARRAY,
316 nr * sizeof *roots);
317 }
318
319 return err;
320 }
321
322 /* Note: X12ARIMA does not give the full covariance matrix: it
323 doesn't calculate covariances between the ARMA terms and any
324 regression terms. But the ARMA covariance matrix is found in
325 the *.acm file and if there are at least two regression terms
326 their covariance matrix is written to an *.rcm file.
327 */
328
get_x12a_vcv(char * fname,const char * path,MODEL * pmod,arma_info * ainfo)329 static int get_x12a_vcv (char *fname, const char *path,
330 MODEL *pmod, arma_info *ainfo)
331 {
332 FILE *fpa = NULL;
333 FILE *fpr = NULL;
334 char *p, *q, line[1024];
335 gretl_matrix *V = NULL;
336 double x;
337 int nc, narma, nr, apos;
338 int i, j, k, li;
339 int err = 0;
340
341 nc = pmod->ncoeff;
342 if (nc == 0) {
343 return 0;
344 }
345
346 narma = ainfo->np + ainfo->P + ainfo->nq + ainfo->Q;
347 nr = nc - narma;
348
349 if (narma >= 2) {
350 gretl_path_compose(fname, MAXLEN, path, ".acm");
351 fpa = gretl_fopen(fname, "r");
352 if (fpa == NULL) {
353 return E_FOPEN;
354 }
355 }
356
357 if (nr >= 2) {
358 gretl_path_compose(fname, MAXLEN, path, ".rcm");
359 fpr = gretl_fopen(fname, "r");
360 if (fpr == NULL) {
361 err = E_FOPEN;
362 goto bailout;
363 }
364 }
365
366 V = gretl_zero_matrix_new(nc, nc);
367 if (V == NULL) {
368 err = E_ALLOC;
369 goto bailout;
370 }
371
372 apos = ainfo->ifc;
373
374 gretl_push_c_numeric_locale();
375
376 /* fill in the ARMA terms first */
377 if (narma >= 2) {
378 int mapos = apos + ainfo->np + ainfo->P;
379
380 i = apos;
381 li = 0;
382 while (fgets(line, sizeof line, fpa)) {
383 /* skip the first two lines */
384 if (li > 1) {
385 p = line + strcspn(line, "+-");
386 for (j=0; j<narma; j++) {
387 x = strtod(p, &q);
388 k = j + apos;
389 if ((i < mapos && k >= mapos) ||
390 (i >= mapos && k < mapos)) {
391 /* flip sign of AR/MA covariances */
392 x = -x;
393 }
394 gretl_matrix_set(V, i, k, x);
395 p = q;
396 }
397 i++;
398 }
399 li++;
400 }
401 } else if (narma > 0) {
402 /* there must be a single AR or MA term */
403 k = ainfo->ifc;
404 x = pmod->sderr[k];
405 gretl_matrix_set(V, k, k, x*x);
406 }
407
408 /* then the regression terms, if any */
409 if (nr >= 2) {
410 int xpos = apos + narma;
411
412 i = ainfo->ifc ? 0 : xpos;
413 li = 0;
414 while (fgets(line, sizeof line, fpr)) {
415 /* skip the first two lines */
416 if (li > 1) {
417 p = line + strcspn(line, "+-");
418 k = ainfo->ifc ? 0 : xpos;
419 for (j=0; j<nr; j++) {
420 x = strtod(p, &q);
421 gretl_matrix_set(V, i, k, x);
422 if (k == 0) {
423 k = xpos;
424 } else {
425 k++;
426 }
427 p = q;
428 }
429 if (i == 0) {
430 i = xpos;
431 } else {
432 i++;
433 }
434 }
435 li++;
436 }
437 } else if (ainfo->ifc) {
438 /* just an intercept */
439 x = pmod->sderr[0];
440 gretl_matrix_set(V, 0, 0, x*x);
441 } else if (ainfo->nexo > 0) {
442 /* single exogenous variable */
443 k = nc - 1;
444 x = pmod->sderr[k];
445 gretl_matrix_set(V, k, k, x*x);
446 }
447
448 gretl_pop_c_numeric_locale();
449
450 err = gretl_model_write_vcv(pmod, V);
451 gretl_matrix_free(V);
452
453 bailout:
454
455 if (fpa != NULL) {
456 fclose(fpa);
457 }
458 if (fpr != NULL) {
459 fclose(fpr);
460 }
461
462 if (err) {
463 fprintf(stderr, "x12a: failed to retrieve covariance matrix\n");
464 }
465
466 return err;
467 }
468
469 /* Below: parse the coefficient estimates and standard errors from
470 the X-12-ARIMA output file foo.est
471 */
472
473 static int
get_estimates(const char * fname,MODEL * pmod,arma_info * ainfo)474 get_estimates (const char *fname, MODEL *pmod, arma_info *ainfo)
475 {
476 double *ar_coeff = pmod->coeff + ainfo->ifc;
477 double *ma_coeff = ar_coeff + ainfo->np + ainfo->P;
478 double *x_coeff = ma_coeff + ainfo->nq + ainfo->Q;
479 double *ar_sderr = pmod->sderr + ainfo->ifc;
480 double *ma_sderr = ar_sderr + ainfo->np + ainfo->P;
481 double *x_sderr = ma_sderr + ainfo->nq + ainfo->Q;
482 FILE *fp;
483 char line[128], word[16];
484 double b, se;
485 int i, j, k;
486 int err = 0;
487
488 fp = gretl_fopen(fname, "r");
489 if (fp == NULL) {
490 fprintf(stderr, "Couldn't read from '%s'\n", fname);
491 return E_FOPEN;
492 }
493
494 for (i=0; i<ainfo->nc; i++) {
495 pmod->coeff[i] = pmod->sderr[i] = NADBL;
496 }
497
498 gretl_push_c_numeric_locale();
499
500 i = j = k = 0;
501
502 while (fgets(line, sizeof line, fp) && i < ainfo->nc) {
503 if (sscanf(line, "%15s", word) == 1) {
504 if (!strcmp(word, "Constant")) {
505 if (sscanf(line, "%*s %*s %lf %lf", &b, &se) == 2) {
506 pmod->coeff[0] = b;
507 pmod->sderr[0] = se;
508 }
509 } else if (!strcmp(word, "User-defined")) {
510 if (sscanf(line, "%*s %*s %lf %lf", &b, &se) == 2) {
511 x_coeff[i] = b;
512 x_sderr[i] = se;
513 i++;
514 }
515 } else if (!strcmp(word, "AR")) {
516 if (sscanf(line, "%*s %*s %*s %*s %lf %lf", &b, &se) == 2) {
517 ar_coeff[j] = b;
518 ar_sderr[j] = se;
519 j++;
520 }
521 } else if (!strcmp(word, "MA")) {
522 if (sscanf(line, "%*s %*s %*s %*s %lf %lf", &b, &se) == 2) {
523 ma_coeff[k] = -b; /* MA sign conventions */
524 ma_sderr[k] = se;
525 k++;
526 }
527 }
528 }
529 }
530
531 gretl_pop_c_numeric_locale();
532
533 fclose(fp);
534
535 for (i=0; i<ainfo->nc; i++) {
536 if (na(pmod->coeff[i]) || na(pmod->sderr[i])) {
537 fprintf(stderr, "Error reading '%s'\n", fname);
538 err = E_DATA;
539 break;
540 }
541 }
542
543 return err;
544 }
545
546 /* On the first pass below we read the in-sample residuals and
547 write to @n_pre the count of pre-sample estimated innovations,
548 if any. On the second pass (if applicable) we pick up the
549 pre-sample values then stop.
550 */
551
read_rsd(double * uh,FILE * fp,const DATASET * dset,int * n_pre)552 static int read_rsd (double *uh, FILE *fp,
553 const DATASET *dset,
554 int *n_pre)
555 {
556 char line[64], date[9];
557 int read_pre = (n_pre == NULL);
558 int t, k, start = 0;
559 int nobs = 0;
560 double x;
561
562 k = 0;
563 while (fgets(line, sizeof line, fp)) {
564 if (*line == '-') {
565 start = 1;
566 continue;
567 }
568 if (start && sscanf(line, "%s %lf", date, &x) == 2) {
569 t = x12_date_to_n(date, dset);
570 if (t < 0) {
571 /* pre-sample */
572 if (read_pre) {
573 /* pass 2: writing pre-sample values */
574 uh[k++] = x;
575 } else {
576 /* pass 1: testing for pre-sample values */
577 *n_pre += 1;
578 }
579 } else if (read_pre) {
580 break;
581 }
582 if (t >= 0 && t < dset->n) {
583 uh[t] = x;
584 nobs++;
585 }
586 }
587 }
588
589 return nobs;
590 }
591
592 /* Parse the residuals from the x12a output file foo.rsd */
593
594 static int
get_uhat(const char * fname,MODEL * pmod,const DATASET * dset)595 get_uhat (const char *fname, MODEL *pmod, const DATASET *dset)
596 {
597 gretl_matrix *pre_innov = NULL;
598 FILE *fp;
599 int nobs = 0;
600 int n_pre = 0;
601 int err = 0;
602
603 fp = gretl_fopen(fname, "r");
604 if (fp == NULL) {
605 fprintf(stderr, "Couldn't read from '%s'\n", fname);
606 return E_FOPEN;
607 }
608
609 gretl_push_c_numeric_locale();
610
611 nobs = read_rsd(pmod->uhat, fp, dset, &n_pre);
612
613 if (nobs > 0 && n_pre > 0) {
614 pre_innov = gretl_matrix_alloc(n_pre, 1);
615 if (pre_innov != NULL) {
616 rewind(fp);
617 read_rsd(pre_innov->val, fp, dset, NULL);
618 gretl_model_set_data(pmod, "pre_innov", pre_innov,
619 GRETL_TYPE_MATRIX, 0);
620 }
621 }
622
623 gretl_pop_c_numeric_locale();
624 fclose(fp);
625
626 if (nobs == 0) {
627 fprintf(stderr, "Error reading '%s'\n", fname);
628 err = E_DATA;
629 }
630
631 return err;
632 }
633
634 static void
populate_x12a_arma_model(MODEL * pmod,const char * path,const DATASET * dset,arma_info * ainfo)635 populate_x12a_arma_model (MODEL *pmod, const char *path,
636 const DATASET *dset,
637 arma_info *ainfo)
638 {
639 char fname[MAXLEN];
640 int err;
641
642 pmod->t1 = ainfo->t1;
643 pmod->t2 = ainfo->t2;
644 pmod->ncoeff = ainfo->nc;
645 pmod->full_n = dset->n;
646
647 err = gretl_model_allocate_storage(pmod);
648 if (err) {
649 pmod->errcode = err;
650 return;
651 }
652
653 err = gretl_path_compose(fname, MAXLEN, path, ".est");
654 if (!err) {
655 err = get_estimates(fname, pmod, ainfo);
656 }
657
658 if (!err) {
659 gretl_path_compose(fname, MAXLEN, path, ".rsd");
660 err = get_uhat(fname, pmod, dset);
661 }
662
663 if (!err) {
664 gretl_path_compose(fname, MAXLEN, path, ".lks");
665 err = get_ll_stats(fname, pmod);
666 }
667
668 if (!err) {
669 gretl_path_compose(fname, MAXLEN, path, ".rts");
670 err = get_roots(fname, pmod, ainfo);
671 }
672
673 if (!err) {
674 /* if this fails we won't count it as a show-stopper */
675 get_x12a_vcv(fname, path, pmod, ainfo);
676 }
677
678 if (err) {
679 fprintf(stderr, "problem reading X-12-ARIMA model info\n");
680 pmod->errcode = err;
681 } else {
682 write_arma_model_stats(pmod, ainfo, dset);
683 }
684 }
685
686 static void
output_series_to_spc(const int * list,const DATASET * dset,int t1,int t2,FILE * fp)687 output_series_to_spc (const int *list, const DATASET *dset,
688 int t1, int t2, FILE *fp)
689 {
690 int i, t, done = 0;
691
692 for (t=t1; t<=t2 && !done; t++) {
693 for (i=1; i<=list[0] && !done; i++) {
694 if (na(dset->Z[list[i]][t])) {
695 fputs(" missingcode=-99999\n", fp);
696 done = 1;
697 }
698 }
699 }
700
701 fputs(" data=(\n", fp);
702
703 for (t=t1; t<=t2; t++) {
704 for (i=1; i<=list[0]; i++) {
705 if (na(dset->Z[list[i]][t])) {
706 fputs("-99999 ", fp);
707 } else {
708 fprintf(fp, "%.15g ", dset->Z[list[i]][t]);
709 }
710 }
711 fputc('\n', fp);
712 }
713
714 fputs(" )\n", fp);
715 }
716
arma_info_get_x_list(arma_info * ainfo)717 static int *arma_info_get_x_list (arma_info *ainfo)
718 {
719 int *xlist = NULL;
720 int start = arma_list_y_position(ainfo);
721 int i;
722
723 xlist = gretl_list_new(ainfo->nexo);
724
725 if (xlist != NULL) {
726 for (i=1; i<=xlist[0]; i++) {
727 xlist[i] = ainfo->alist[i + start];
728 }
729 }
730
731 return xlist;
732 }
733
734 static void
make_x12a_date_string(int t,const DATASET * dset,char * str)735 make_x12a_date_string (int t, const DATASET *dset, char *str)
736 {
737 if (non_yearly_frequency(dset->pd)) {
738 int maj = t / dset->pd + 1;
739 int min = t % dset->pd + 1;
740
741 sprintf(str, "%d.%d", maj, min);
742 } else {
743 ntolabel(str, t, dset);
744 gretl_charsub(str, ':', '.');
745 }
746 }
747
x12_pdq_string(arma_info * ainfo,FILE * fp)748 static void x12_pdq_string (arma_info *ainfo, FILE *fp)
749 {
750 fputc('(', fp);
751
752 if (ainfo->pmask == NULL) {
753 fprintf(fp, "%d", ainfo->p);
754 } else {
755 int i;
756
757 fputc('[', fp);
758 for (i=0; i<ainfo->p; i++) {
759 if (AR_included(ainfo, i)) {
760 fprintf(fp, "%d", i+1);
761 if (i < ainfo->p - 1) {
762 fputc(' ', fp);
763 }
764 }
765 }
766 fputc(']', fp);
767 }
768
769 fprintf(fp, " %d ", ainfo->d);
770
771 if (ainfo->qmask == NULL) {
772 fprintf(fp, "%d", ainfo->q);
773 } else {
774 int i;
775
776 fputc('[', fp);
777 for (i=0; i<ainfo->q; i++) {
778 if (MA_included(ainfo, i)) {
779 fprintf(fp, "%d", i+1);
780 if (i < ainfo->q - 1) {
781 fputc(' ', fp);
782 }
783 }
784 }
785 fputc(']', fp);
786 }
787
788 fputc(')', fp);
789 }
790
write_arma_spc_file(const char * fname,const DATASET * dset,arma_info * ainfo,int pdmax,gretlopt opt)791 static int write_arma_spc_file (const char *fname,
792 const DATASET *dset,
793 arma_info *ainfo, int pdmax,
794 gretlopt opt)
795 {
796 int maxobs = pdmax * 50;
797 int maxfcast = pdmax * 5;
798 int ylist[2];
799 int *xlist = NULL;
800 FILE *fp;
801 char datestr[12];
802 int nfcast = 0;
803 int t1 = ainfo->t1;
804 int tmax;
805 int i;
806
807 if (ainfo->nexo > 0) {
808 xlist = arma_info_get_x_list(ainfo);
809 if (xlist == NULL) {
810 return E_ALLOC;
811 }
812 }
813
814 fp = gretl_fopen(fname, "w");
815 if (fp == NULL) {
816 fprintf(stderr, "Couldn't write to '%s'\n", fname);
817 return 1;
818 }
819
820 gretl_push_c_numeric_locale();
821
822 fprintf(fp, "series{\n title=\"%s\"\n period=%d\n",
823 dset->varname[ainfo->yno], dset->pd);
824
825 t1 -= ainfo->maxlag;
826
827 make_x12a_date_string(t1, dset, datestr);
828 fprintf(fp, " start=%s\n", datestr);
829
830 ylist[0] = 1;
831 ylist[1] = ainfo->yno;
832
833 tmax = ainfo->t2;
834
835 if ((opt & OPT_F) && ainfo->t2 < dset->n - 1) {
836 int nobs;
837
838 tmax = dset->n - 1;
839 nobs = tmax - ainfo->t1 + 1; /* t1? */
840 if (nobs > maxobs) {
841 tmax -= nobs - maxobs;
842 }
843 nfcast = tmax - ainfo->t2;
844 if (nfcast > maxfcast) {
845 tmax -= nfcast - maxfcast;
846 nfcast -= nfcast - maxfcast;
847 }
848 #if 0
849 fprintf(stderr, "x12a: doing forecast: nfcast = %d\n", nfcast);
850 #endif
851 }
852
853 output_series_to_spc(ylist, dset, t1, tmax, fp);
854
855 if (tmax > ainfo->t2) {
856 make_x12a_date_string(ainfo->t2, dset, datestr);
857 fprintf(fp, " span = ( , %s)\n", datestr);
858 }
859
860 fputs("}\n", fp);
861
862 /* regression specification */
863 fputs("Regression {\n", fp);
864 if (ainfo->ifc) {
865 fputs(" variables = (const)\n", fp);
866 }
867 if (ainfo->nexo > 0) {
868 fputs(" user = ( ", fp);
869 for (i=1; i<=xlist[0]; i++) {
870 fprintf(fp, "%s ", dset->varname[xlist[i]]);
871 }
872 fputs(")\n", fp);
873 output_series_to_spc(xlist, dset, t1, tmax, fp);
874 free(xlist);
875 }
876 fputs("}\n", fp);
877
878 /* arima specification */
879 fputs("arima {\n model = ", fp);
880 x12_pdq_string(ainfo, fp);
881 if (ainfo->P > 0 || ainfo->D > 0 || ainfo->Q > 0) {
882 fprintf(fp, "(%d %d %d)\n}\n", ainfo->P, ainfo->D, ainfo->Q);
883 } else {
884 fputs("\n}\n", fp);
885 }
886
887 fputs("estimate {\n", fp);
888
889 /* could enable here: "tol = XX, maxiter = NN"
890 the default is tol = 1.0e-5, maxiter = 200
891 */
892
893 if (opt & OPT_V) {
894 fputs(" print = (acm itr lkf lks mdl est rts rcm)\n", fp);
895 } else {
896 fputs(" print = (acm lkf lks mdl est rts rcm)\n", fp);
897 }
898
899 fputs(" save = (rsd est lks acm rts rcm)\n", fp);
900
901 if (opt & OPT_C) {
902 fputs(" exact = none\n", fp);
903 }
904
905 fputs("}\n", fp);
906
907 if (nfcast > 0) {
908 fputs("forecast {\n save = (ftr)\n", fp);
909 fprintf(fp, " maxlead = %d\n}\n", nfcast);
910 }
911
912 gretl_pop_c_numeric_locale();
913
914 fclose(fp);
915
916 return 0;
917 }
918
print_x12a_error_file(const char * fname,PRN * prn)919 static void print_x12a_error_file (const char *fname, PRN *prn)
920 {
921 FILE *fp = gretl_fopen(fname, "r");
922
923 if (fp != NULL) {
924 char line[512];
925
926 while (fgets(line, sizeof line, fp)) {
927 pputs(prn, line);
928 }
929 fclose(fp);
930 }
931 }
932
delete_old_files(const char * path)933 static void delete_old_files (const char *path)
934 {
935 const char *exts[] = {
936 "acm", "itr", "lkf", "lks", "mdl", "est", "rts",
937 "rcm", "rsd", "ftr", "err", "log", "out", NULL
938 };
939 char old[MAXLEN];
940 int i, n = strlen(path);
941
942 for (i=0; exts[i] != NULL; i++) {
943 *old = '\0';
944 strncat(old, path, n - 3);
945 strcat(old, exts[i]);
946 gretl_remove(old);
947 }
948 }
949
x12a_maybe_allow_missvals(arma_info * ainfo)950 static void x12a_maybe_allow_missvals (arma_info *ainfo)
951 {
952 if (arma_exact_ml(ainfo)) {
953 ainfo->pflags |= ARMA_NAOK;
954 }
955 }
956
arma_x12_model(const int * list,const int * pqspec,const DATASET * dset,int pdmax,gretlopt opt,PRN * prn)957 MODEL arma_x12_model (const int *list, const int *pqspec,
958 const DATASET *dset, int pdmax,
959 gretlopt opt, PRN *prn)
960 {
961 int verbose = (opt & OPT_V);
962 const char *prog = gretl_x12_arima();
963 const char *workdir = gretl_x12_arima_dir();
964 char yname[VNAMELEN], path[MAXLEN];
965 MODEL armod;
966 arma_info ainfo_s, *ainfo;
967 int missv = 0, misst = 0;
968 #ifdef WIN32
969 char *cmd;
970 #endif
971 int err = 0;
972
973 #if 0
974 if (dset->t2 < dset->n - 1) {
975 /* FIXME this is temporary (OPT_F -> generate forecast) */
976 opt |= OPT_F;
977 }
978 #endif
979
980 ainfo = &ainfo_s;
981 arma_info_init(ainfo, opt | OPT_X, pqspec, dset);
982 if (opt & OPT_V) {
983 ainfo->prn = prn;
984 }
985 gretl_model_init(&armod, dset);
986
987 ainfo->alist = gretl_list_copy(list);
988 if (ainfo->alist == NULL) {
989 err = E_ALLOC;
990 }
991
992 if (!err) {
993 err = arma_check_list(ainfo, dset, opt);
994 }
995
996 if (err) {
997 armod.errcode = err;
998 goto bailout;
999 }
1000
1001 /* calculate maximum lag */
1002 calc_max_lag(ainfo);
1003
1004 x12a_maybe_allow_missvals(ainfo);
1005
1006 /* adjust sample range if need be */
1007 armod.errcode = arma_adjust_sample(ainfo, dset, &missv, &misst);
1008 if (armod.errcode) {
1009 goto bailout;
1010 } else if (missv > 0) {
1011 set_arma_missvals(ainfo);
1012 }
1013
1014 ainfo->y = (double *) dset->Z[ainfo->yno]; /* it's really const */
1015 strcpy(yname, dset->varname[ainfo->yno]);
1016
1017 /* write out an .spc file */
1018 sprintf(path, "%s%c%s.spc", workdir, SLASH, yname);
1019 write_arma_spc_file(path, dset, ainfo, pdmax, opt);
1020
1021 /* remove any files from an old run, in case of error */
1022 delete_old_files(path);
1023
1024 /* run the program */
1025 #ifdef WIN32
1026 cmd = g_strdup_printf("\"%s\" %s -r -p -q", prog, yname);
1027 err = win_run_sync(cmd, workdir);
1028 g_free(cmd);
1029 #else
1030 err = tramo_x12a_spawn(workdir, prog, yname, "-r", "-p", "-q", "-n", NULL);
1031 #endif
1032
1033 if (!err) {
1034 sprintf(path, "%s%c%s", workdir, SLASH, yname);
1035 populate_x12a_arma_model(&armod, path, dset, ainfo);
1036 if (verbose && !armod.errcode) {
1037 print_iterations(path, ainfo->prn);
1038 }
1039 if (!armod.errcode) {
1040 if (gretl_in_gui_mode()) {
1041 add_unique_output_file(&armod, path);
1042 }
1043 gretl_model_set_int(&armod, "arma_flags", (int) ainfo->flags);
1044 }
1045 } else {
1046 armod.errcode = E_UNSPEC;
1047 gretl_errmsg_set(_("Failed to execute x12arima"));
1048 }
1049
1050 if (armod.errcode && ainfo->prn != NULL) {
1051 sprintf(path, "%s%c%s.err", workdir, SLASH, yname);
1052 print_x12a_error_file(path, ainfo->prn);
1053 }
1054
1055 if (armod.errcode == 0) {
1056 transcribe_extra_info(ainfo, &armod);
1057 }
1058
1059 bailout:
1060
1061 arma_info_cleanup(ainfo);
1062
1063 return armod;
1064 }
1065