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 "usermat.h"
20
21 #define MAX_ARMA_ORDER 128
22 #define MAX_ARIMA_DIFF 2
23
24 #define SAMPLE_DEBUG 0
25
26 static void
27 real_arima_difference_series (double *dx, const double *x,
28 int t1, int t2, int *delta,
29 int k);
30
31 static void
arma_info_init(arma_info * ainfo,gretlopt opt,const int * pqspec,const DATASET * dset)32 arma_info_init (arma_info *ainfo, gretlopt opt,
33 const int *pqspec, const DATASET *dset)
34 {
35 ainfo->yno = 0;
36 ainfo->flags = 0;
37 ainfo->pflags = 0;
38 ainfo->init = 0;
39 ainfo->alist = NULL;
40
41 if (opt & OPT_X) {
42 /* we got --x-12-arima */
43 ainfo->flags |= ARMA_X12A;
44 }
45
46 if (!(opt & OPT_C)) {
47 /* we didn't get --conditional */
48 ainfo->flags |= ARMA_EXACT;
49 }
50
51 ainfo->ll = NADBL;
52
53 ainfo->pqspec = pqspec;
54 ainfo->pmask = NULL;
55 ainfo->qmask = NULL;
56
57 ainfo->p = 0;
58 ainfo->d = 0;
59 ainfo->q = 0;
60 ainfo->P = 0;
61 ainfo->D = 0;
62 ainfo->Q = 0;
63
64 ainfo->np = 0;
65 ainfo->nq = 0;
66
67 ainfo->maxlag = 0;
68 ainfo->ifc = 0;
69 ainfo->nexo = 0;
70 ainfo->nc = 0;
71
72 ainfo->t1 = dset->t1;
73 ainfo->t2 = dset->t2;
74 ainfo->pd = dset->pd;
75 ainfo->T = 0;
76 ainfo->r0 = 0;
77
78 ainfo->fncount = 0;
79 ainfo->grcount = 0;
80
81 ainfo->y = NULL;
82 ainfo->e = NULL;
83 ainfo->Z = NULL;
84 ainfo->yscale = 1.0;
85 ainfo->yshift = 0.0;
86
87 ainfo->xlist = NULL;
88 ainfo->misslist = NULL;
89 ainfo->xstats = NULL;
90 ainfo->dX = NULL;
91 ainfo->G = NULL;
92 ainfo->V = NULL;
93
94 ainfo->n_aux = 0;
95 ainfo->aux = NULL;
96
97 ainfo->prn = NULL;
98 }
99
arma_info_cleanup(arma_info * ainfo)100 static void arma_info_cleanup (arma_info *ainfo)
101 {
102 free(ainfo->alist);
103 free(ainfo->pmask);
104 free(ainfo->qmask);
105 free(ainfo->e);
106 free(ainfo->Z);
107 free(ainfo->xlist);
108 free(ainfo->misslist);
109
110 gretl_matrix_free(ainfo->xstats);
111 gretl_matrix_free(ainfo->dX);
112 gretl_matrix_free(ainfo->G);
113 gretl_matrix_free(ainfo->V);
114
115 if (arima_ydiff(ainfo)) {
116 free(ainfo->y);
117 }
118
119 doubles_array_free(ainfo->aux, ainfo->n_aux);
120 }
121
122 enum {
123 AR_MASK,
124 MA_MASK
125 };
126
127 /* Create a mask for skipping certain intermediate lags,
128 AR or MA. This function also sets ainfo->np and ainfo->nq,
129 which record the actual number of non-seasonal AR and MA
130 lags used.
131 */
132
mask_from_list(const int * list,arma_info * ainfo,int m,int * err)133 static char *mask_from_list (const int *list,
134 arma_info *ainfo,
135 int m, int *err)
136 {
137 int mlen = (m == AR_MASK)? ainfo->p : ainfo->q;
138 int nv = 0, nmax = 0;
139 char *mask;
140 int i, k;
141
142 mask = malloc(mlen + 1);
143 if (mask == NULL) {
144 *err = E_ALLOC;
145 return NULL;
146 }
147
148 for (i=0; i<mlen; i++) {
149 mask[i] = '0';
150 }
151 mask[mlen] = '\0';
152
153 for (i=1; i<=list[0]; i++) {
154 k = list[i];
155 if (k > 0) {
156 mask[k-1] = '1';
157 nv++;
158 if (k > nmax) {
159 nmax = k;
160 }
161 }
162 }
163
164 if (m == AR_MASK) {
165 ainfo->p = nmax;
166 ainfo->np = nv;
167 } else {
168 ainfo->q = nmax;
169 ainfo->nq = nv;
170 }
171
172 if (nv == 0) {
173 free(mask);
174 mask = NULL;
175 }
176
177 return mask;
178 }
179
arma_make_masks(arma_info * ainfo)180 static int arma_make_masks (arma_info *ainfo)
181 {
182 int *plist = NULL, *qlist = NULL;
183 int err = 0;
184
185 if (ainfo->pqspec != NULL) {
186 if (gretl_list_has_separator(ainfo->pqspec)) {
187 gretl_list_split_on_separator(ainfo->pqspec, &plist, &qlist);
188 } else {
189 plist = gretl_list_copy(ainfo->pqspec);
190 }
191 }
192
193 if (ainfo->p > 0) {
194 ainfo->np = ainfo->p;
195 if (plist != NULL && plist[0] > 0) {
196 ainfo->pmask = mask_from_list(plist, ainfo, AR_MASK, &err);
197 }
198 }
199
200 if (ainfo->q > 0 && !err) {
201 ainfo->nq = ainfo->q;
202 if (qlist != NULL && qlist[0] > 0) {
203 ainfo->qmask = mask_from_list(qlist, ainfo, MA_MASK, &err);
204 }
205 }
206
207 free(plist);
208 free(qlist);
209
210 return err;
211 }
212
arma_list_y_position(arma_info * ainfo)213 int arma_list_y_position (arma_info *ainfo)
214 {
215 int ypos;
216
217 if (arma_is_arima(ainfo)) {
218 ypos = arma_has_seasonal(ainfo) ? 9 : 5;
219 } else {
220 ypos = arma_has_seasonal(ainfo) ? 7 : 4;
221 }
222
223 return ypos;
224 }
225
arima_integrate(double * dx,const double * x,int t1,int t2,int d,int D,int s)226 static int arima_integrate (double *dx, const double *x,
227 int t1, int t2, int d, int D, int s)
228 {
229 double *ix;
230 int *c;
231 int k = d + s * D;
232 int i, t;
233
234 ix = malloc((t2 + 1) * sizeof *ix);
235 if (ix == NULL) {
236 return E_ALLOC;
237 }
238
239 c = arima_delta_coeffs(d, D, s);
240 if (c == NULL) {
241 free(ix);
242 return E_ALLOC;
243 }
244
245 for (t=0; t<t1; t++) {
246 ix[t] = 0.0;
247 }
248
249 for (t=t1; t<=t2; t++) {
250 ix[t] = dx[t];
251 for (i=0; i<k; i++) {
252 if (c[i] != 0) {
253 ix[t] += c[i] * x[t-i-1];
254 }
255 }
256 }
257
258 /* transcribe integrated result back into "dx" */
259 for (t=0; t<=t2; t++) {
260 if (t < t1) {
261 dx[t] = NADBL;
262 } else {
263 dx[t] = ix[t];
264 }
265 }
266
267 free(ix);
268 free(c);
269
270 return 0;
271 }
272
ainfo_data_to_model(arma_info * ainfo,MODEL * pmod)273 static void ainfo_data_to_model (arma_info *ainfo, MODEL *pmod)
274 {
275 pmod->ifc = ainfo->ifc;
276 pmod->dfn = ainfo->nc - pmod->ifc;
277 pmod->dfd = pmod->nobs - pmod->dfn;
278 pmod->ncoeff = ainfo->nc;
279
280 if (arma_has_seasonal(ainfo)) {
281 gretl_model_set_int(pmod, "arma_P", ainfo->P);
282 gretl_model_set_int(pmod, "arma_Q", ainfo->Q);
283 gretl_model_set_int(pmod, "arma_pd", ainfo->pd);
284 }
285
286 if (ainfo->d > 0 || ainfo->D > 0) {
287 gretl_model_set_int(pmod, "arima_d", ainfo->d);
288 gretl_model_set_int(pmod, "arima_D", ainfo->D);
289 }
290
291 if (ainfo->nexo > 0) {
292 gretl_model_set_int(pmod, "armax", 1);
293 }
294
295 if (ainfo->pmask != NULL) {
296 gretl_model_set_string_as_data(pmod, "pmask",
297 gretl_strdup(ainfo->pmask));
298 }
299
300 if (ainfo->qmask != NULL) {
301 gretl_model_set_string_as_data(pmod, "qmask",
302 gretl_strdup(ainfo->qmask));
303 }
304 }
305
arma_depvar_stats(MODEL * pmod,arma_info * ainfo,const DATASET * dset)306 static void arma_depvar_stats (MODEL *pmod, arma_info *ainfo,
307 const DATASET *dset)
308 {
309 if (arma_is_arima(ainfo) && !arima_ydiff(ainfo)) {
310 /* calculate differenced y for stats */
311 int d = ainfo->d, D = ainfo->D;
312 int T = pmod->t2 - pmod->t1 + 1;
313 double *dy = malloc(T * sizeof *dy);
314 int *delta = arima_delta_coeffs(d, D, ainfo->pd);
315
316 if (dy != NULL && delta != NULL) {
317 int k = d + ainfo->pd * D;
318
319 real_arima_difference_series(dy, dset->Z[ainfo->yno],
320 pmod->t1, pmod->t2, delta, k);
321 pmod->ybar = gretl_mean(0, T - 1, dy);
322 pmod->sdy = gretl_stddev(0, T - 1, dy);
323 }
324 free(dy);
325 free(delta);
326 } else {
327 pmod->ybar = gretl_mean(pmod->t1, pmod->t2, ainfo->y);
328 pmod->sdy = gretl_stddev(pmod->t1, pmod->t2, ainfo->y);
329 }
330 }
331
handle_null_model(MODEL * pmod,arma_info * ainfo)332 static void handle_null_model (MODEL *pmod, arma_info *ainfo)
333 {
334 int full_n = pmod->full_n;
335
336 pmod->ncoeff = 1;
337 pmod->full_n = 0;
338 pmod->errcode = gretl_model_allocate_storage(pmod);
339 pmod->full_n = full_n;
340
341 if (!pmod->errcode) {
342 gretl_model_set_int(pmod, "null-model", 1);
343 pmod->coeff[0] = 0.0;
344 pmod->sigma = pmod->sdy;
345 }
346 }
347
348 #define USE_ARIMA_INTEGRATE 1
349
350 /* write the various statistics from ARMA estimation into
351 a gretl MODEL struct */
352
write_arma_model_stats(MODEL * pmod,arma_info * ainfo,const DATASET * dset)353 void write_arma_model_stats (MODEL *pmod, arma_info *ainfo,
354 const DATASET *dset)
355 {
356 double mean_error;
357 int do_criteria = 1;
358 int t;
359
360 pmod->ci = ARMA;
361 ainfo_data_to_model(ainfo, pmod);
362
363 free(pmod->list);
364 pmod->list = gretl_list_copy(ainfo->alist);
365
366 if (!arma_least_squares(ainfo)) {
367 arma_depvar_stats(pmod, ainfo, dset);
368 }
369
370 mean_error = pmod->ess = 0.0;
371
372 for (t=pmod->t1; t<=pmod->t2; t++) {
373 if (!na(ainfo->y[t]) && !na(pmod->uhat[t])) {
374 #if USE_ARIMA_INTEGRATE == 0
375 if (arma_is_arima(ainfo) && arima_ydiff(ainfo)) {
376 pmod->yhat[t] = Z[ainfo->yno][t] - pmod->uhat[t];
377 }
378 #else
379 pmod->yhat[t] = ainfo->y[t] - pmod->uhat[t];
380 #endif
381 pmod->ess += pmod->uhat[t] * pmod->uhat[t];
382 mean_error += pmod->uhat[t];
383 }
384 }
385
386 #if USE_ARIMA_INTEGRATE
387 if (arma_is_arima(ainfo) && arima_ydiff(ainfo)) {
388 arima_integrate(pmod->yhat, dset->Z[ainfo->yno],
389 pmod->t1, pmod->t2,
390 ainfo->d, ainfo->D, ainfo->pd);
391 }
392 #endif
393
394 mean_error /= pmod->nobs;
395 if (arma_least_squares(ainfo) && pmod->ifc &&
396 mean_error < 1.0e-15) {
397 mean_error = 0.0;
398 }
399 gretl_model_set_double(pmod, "mean_error", mean_error);
400
401 if (na(pmod->sigma)) {
402 /* in x12a or native exact cases this is already done */
403 pmod->sigma = sqrt(pmod->ess / pmod->nobs);
404 }
405
406 /* R-squared as suggested by J\"org Breitung */
407 pmod->rsq = gretl_corr_rsq(pmod->t1, pmod->t2, dset->Z[ainfo->yno], pmod->yhat);
408 pmod->adjrsq = 1.0 - ((1.0 - pmod->rsq) * (pmod->t2 - pmod->t1) /
409 (double) pmod->dfd);
410 pmod->fstt = pmod->chisq = NADBL;
411 pmod->tss = NADBL;
412
413 if (arma_least_squares(ainfo)) {
414 /* not applicable */
415 do_criteria = 0;
416 } else if (arma_by_x12a(ainfo) && !na(pmod->criterion[C_AIC])) {
417 /* already given by x12a */
418 do_criteria = 0;
419 }
420
421 if (do_criteria) {
422 mle_criteria(pmod, 1);
423 }
424
425 if (!pmod->errcode && pmod->ncoeff == 0) {
426 handle_null_model(pmod, ainfo);
427 }
428
429 if (!pmod->errcode) {
430 gretl_model_add_arma_varnames(pmod, dset, ainfo->yno,
431 ainfo->p, ainfo->q,
432 ainfo->pmask, ainfo->qmask,
433 ainfo->P, ainfo->Q,
434 ainfo->nexo);
435 }
436 }
437
calc_max_lag(arma_info * ainfo)438 static void calc_max_lag (arma_info *ainfo)
439 {
440 if (arma_exact_ml(ainfo)) {
441 ainfo->maxlag = ainfo->d + ainfo->D * ainfo->pd;
442 } else {
443 /* conditional ML */
444 int pmax = ainfo->p + ainfo->P * ainfo->pd;
445 int dmax = ainfo->d + ainfo->D * ainfo->pd;
446
447 ainfo->maxlag = pmax + dmax;
448 }
449
450 #if SAMPLE_DEBUG
451 fprintf(stderr, "calc_max_lag: ainfo->maxlag = %d\n", ainfo->maxlag);
452 #endif
453 }
454
arma_adjust_sample(arma_info * ainfo,const DATASET * dset,int * missv,int * misst)455 static int arma_adjust_sample (arma_info *ainfo,
456 const DATASET *dset,
457 int *missv, int *misst)
458 {
459 int *list = ainfo->alist;
460 int ypos = arma_list_y_position(ainfo);
461 int t0, t1 = dset->t1, t2 = dset->t2;
462 int i, vi, vlmax, k, t;
463 int missing;
464 int err = 0;
465
466 #if SAMPLE_DEBUG
467 fprintf(stderr, "arma_adjust_sample: at start, t1=%d, t2=%d, maxlag = %d\n",
468 t1, t2, ainfo->maxlag);
469 #endif
470
471 t0 = t1 - ainfo->maxlag;
472 if (t0 < 0) {
473 t1 -= t0;
474 }
475
476 /* list position of last var to check for lags */
477 if (arma_xdiff(ainfo)) {
478 vlmax = list[0];
479 } else {
480 vlmax = ypos;
481 }
482
483 /* advance the starting point if need be */
484
485 for (t=t1; t<=t2; t++) {
486 missing = 0;
487 for (i=ypos; i<=list[0] && !missing; i++) {
488 vi = list[i];
489 if (na(dset->Z[vi][t])) {
490 /* current value missing */
491 missing = 1;
492 }
493 if (i <= vlmax) {
494 for (k=1; k<=ainfo->maxlag && !missing; k++) {
495 if (na(dset->Z[vi][t-k])) {
496 /* lagged value missing */
497 missing = 1;
498 }
499 }
500 }
501 }
502 if (missing) {
503 t1++;
504 } else {
505 break;
506 }
507 }
508
509 /* retard the ending point if need be */
510
511 for (t=t2; t>=t1; t--) {
512 missing = 0;
513 for (i=ypos; i<=list[0] && !missing; i++) {
514 vi = list[i];
515 if (na(dset->Z[vi][t])) {
516 missing = 1;
517 }
518 }
519 if (missing) {
520 t2--;
521 } else {
522 break;
523 }
524 }
525
526 if (t2 < t1) {
527 gretl_errmsg_set(_("No usable data were found"));
528 return E_MISSDATA;
529 }
530
531 missing = 0;
532
533 /* check for missing obs within the adjusted sample range */
534 for (t=t1; t<t2; t++) {
535 int tmiss = 0;
536
537 for (i=ypos; i<=list[0]; i++) {
538 vi = list[i];
539 if (na(dset->Z[vi][t])) {
540 if (missv != NULL && misst != NULL && *missv == 0) {
541 /* record info on first missing obs */
542 *missv = vi;
543 *misst = t + 1;
544 }
545 tmiss = 1;
546 }
547 }
548 if (tmiss) {
549 missing++;
550 }
551 }
552
553 if (missing > 0 && !arma_na_ok(ainfo)) {
554 err = E_MISSDATA;
555 }
556
557 if (!err) {
558 ainfo->fullT = t2 - t1 + 1;
559 ainfo->T = ainfo->fullT - missing;
560 if (ainfo->T <= ainfo->nc) {
561 /* insufficient observations */
562 err = E_DF;
563 }
564 }
565
566 if (!err) {
567 #if SAMPLE_DEBUG
568 fprintf(stderr, "arma_adjust_sample: at end, t1=%d, t2=%d\n",
569 t1, t2);
570 #endif
571 ainfo->t1 = t1;
572 ainfo->t2 = t2;
573 }
574
575 return err;
576 }
577
578 /* remove the intercept from list of regressors */
579
arma_remove_const(arma_info * ainfo,const DATASET * dset)580 static int arma_remove_const (arma_info *ainfo,
581 const DATASET *dset)
582 {
583 int *list = ainfo->alist;
584 int seasonal = arma_has_seasonal(ainfo);
585 int diffs = arma_is_arima(ainfo);
586 int xstart, ret = 0;
587 int i, j;
588
589 if (diffs) {
590 xstart = (seasonal)? 10 : 6;
591 } else {
592 xstart = (seasonal)? 8 : 5;
593 }
594
595 for (i=xstart; i<=list[0]; i++) {
596 if (list[i] == 0 || true_const(list[i], dset)) {
597 for (j=i; j<list[0]; j++) {
598 list[j] = list[j+1];
599 }
600 list[0] -= 1;
601 ret = 1;
602 break;
603 }
604 }
605
606 return ret;
607 }
608
check_arma_sep(arma_info * ainfo,int sep1)609 static int check_arma_sep (arma_info *ainfo, int sep1)
610 {
611 int *list = ainfo->alist;
612 int sep2 = (sep1 == 3)? 6 : 8;
613 int i, err = 0;
614
615 for (i=sep1+1; i<=list[0]; i++) {
616 if (list[i] == LISTSEP) {
617 if (i == sep2) {
618 /* there's a second list separator in the right place:
619 we've got a seasonal specification */
620 set_arma_has_seasonal(ainfo);
621 } else {
622 err = 1;
623 }
624 }
625 }
626
627 if (!err && sep1 == 4) {
628 /* check for apparent but not "real" arima spec */
629 if (arma_has_seasonal(ainfo)) {
630 if (list[2] == 0 && list[6] == 0) {
631 gretl_list_delete_at_pos(list, 2);
632 gretl_list_delete_at_pos(list, 5);
633 unset_arma_is_arima(ainfo);
634 }
635 } else {
636 if (list[2] == 0) {
637 gretl_list_delete_at_pos(list, 2);
638 unset_arma_is_arima(ainfo);
639 }
640 }
641 }
642
643 return err;
644 }
645
arma_add_xlist(arma_info * ainfo,int ypos)646 static int arma_add_xlist (arma_info *ainfo, int ypos)
647 {
648 int i, err = 0;
649
650 ainfo->xlist = gretl_list_new(ainfo->nexo);
651
652 if (ainfo->xlist == NULL) {
653 err = E_ALLOC;
654 } else {
655 for (i=1; i<=ainfo->nexo; i++) {
656 ainfo->xlist[i] = ainfo->alist[ypos + i];
657 }
658 }
659
660 return err;
661 }
662
663 #define count_arma_coeffs(a) (a->ifc + a->np + a->nq + a->P + a->Q + a->nexo)
664
check_arma_list(arma_info * ainfo,const DATASET * dset,gretlopt opt)665 static int check_arma_list (arma_info *ainfo,
666 const DATASET *dset,
667 gretlopt opt)
668 {
669 int *list = ainfo->alist;
670 int ypos = arma_has_seasonal(ainfo) ? 7 : 4;
671 int armax = (list[0] > ypos);
672 int hadconst = 0;
673 int err = 0;
674
675 if (list[1] < 0 || list[1] > MAX_ARMA_ORDER) {
676 /* non-seasonal AR order out of bounds */
677 err = 1;
678 } else if (list[2] < 0 || list[2] > MAX_ARMA_ORDER) {
679 /* non-seasonal MA order out of bounds */
680 err = 1;
681 }
682
683 if (!err) {
684 ainfo->p = list[1];
685 ainfo->q = list[2];
686 }
687
688 if (!err && arma_has_seasonal(ainfo)) {
689 if (list[0] < 7) {
690 err = 1;
691 } else if (list[4] < 0 || list[4] > MAX_ARMA_ORDER) {
692 /* seasonal AR order out of bounds */
693 err = 1;
694 } else if (list[5] < 0 || list[5] > MAX_ARMA_ORDER) {
695 /* seasonal MA order out of bounds */
696 err = 1;
697 }
698 }
699
700 if (!err && arma_has_seasonal(ainfo)) {
701 ainfo->P = list[4];
702 ainfo->Q = list[5];
703 }
704
705 /* now that we have p and q we can check for masked lags */
706
707 if (!err) {
708 err = arma_make_masks(ainfo);
709 }
710
711 /* If there's an explicit constant in the list here, we'll remove
712 it, since it is added implicitly later. But if we're supplied
713 with OPT_N (meaning: no intercept) we'll flag this by
714 setting ifc = 0. Also, if the user gave an armax list
715 (specifying regressors) we'll respect the absence of a constant
716 from that list by setting ifc = 0.
717 */
718
719 if (!err) {
720 if (armax) {
721 hadconst = arma_remove_const(ainfo, dset);
722 }
723 if ((opt & OPT_N) || (armax && !hadconst)) {
724 ; /* no constant present */
725 } else {
726 ainfo->ifc = 1;
727 }
728 }
729
730 if (err) {
731 gretl_errmsg_set(_("Error in arma command"));
732 } else {
733 ainfo->nexo = list[0] - ypos;
734 ainfo->nc = count_arma_coeffs(ainfo);
735 ainfo->yno = list[ypos];
736 if (ainfo->nexo > 0) {
737 err = arma_add_xlist(ainfo, ypos);
738 }
739 }
740
741 return err;
742 }
743
check_arima_list(arma_info * ainfo,const DATASET * dset,gretlopt opt)744 static int check_arima_list (arma_info *ainfo,
745 const DATASET *dset,
746 gretlopt opt)
747 {
748 int *list = ainfo->alist;
749 int ypos = arma_has_seasonal(ainfo) ? 9 : 5;
750 int armax = (list[0] > ypos);
751 int hadconst = 0;
752 int err = 0;
753
754 if (list[1] < 0 || list[1] > MAX_ARMA_ORDER) {
755 err = 1;
756 } else if (list[2] < 0 || list[2] > MAX_ARIMA_DIFF) {
757 err = 1;
758 } else if (list[3] < 0 || list[3] > MAX_ARMA_ORDER) {
759 err = 1;
760 }
761
762 if (!err) {
763 ainfo->p = list[1];
764 ainfo->d = list[2];
765 ainfo->q = list[3];
766 }
767
768 if (!err && arma_has_seasonal(ainfo)) {
769 if (list[0] < 9) {
770 err = 1;
771 } else if (list[5] < 0 || list[5] > MAX_ARMA_ORDER) {
772 err = 1;
773 } else if (list[6] < 0 || list[6] > MAX_ARIMA_DIFF) {
774 err = 1;
775 } else if (list[7] < 0 || list[7] > MAX_ARMA_ORDER) {
776 err = 1;
777 }
778 }
779
780 if (!err && arma_has_seasonal(ainfo)) {
781 ainfo->P = list[5];
782 ainfo->D = list[6];
783 ainfo->Q = list[7];
784 }
785
786 /* now that we have p and q we can check for masked lags */
787
788 if (!err) {
789 err = arma_make_masks(ainfo);
790 }
791
792 /* If there's an explicit constant in the list here, we'll remove
793 it, since it is added implicitly later. But if we're supplied
794 with OPT_N (meaning: no intercept) we'll flag this by
795 setting ifc = 0. Also, if the user gave an armax list
796 (specifying regressors) we'll respect the absence of a constant
797 from that list by setting ifc = 0.
798 */
799
800 if (!err) {
801 if (armax) {
802 hadconst = arma_remove_const(ainfo, dset);
803 }
804 if ((opt & OPT_N) || (armax && !hadconst)) {
805 ;
806 } else {
807 ainfo->ifc = 1;
808 }
809 }
810
811 if (err) {
812 gretl_errmsg_set(_("Error in arma command"));
813 } else {
814 ainfo->nexo = list[0] - ypos;
815 ainfo->nc = count_arma_coeffs(ainfo);
816 ainfo->yno = list[ypos];
817 if (ainfo->nexo > 0) {
818 err = arma_add_xlist(ainfo, ypos);
819 }
820 }
821
822 return err;
823 }
824
arma_check_list(arma_info * ainfo,const DATASET * dset,gretlopt opt)825 static int arma_check_list (arma_info *ainfo,
826 const DATASET *dset,
827 gretlopt opt)
828 {
829 int *list = ainfo->alist;
830 int sep1 = gretl_list_separator_position(list);
831 int err = 0;
832
833 if (sep1 == 3) {
834 if (list[0] < 4) {
835 err = E_PARSE;
836 }
837 } else if (sep1 == 4) {
838 if (list[0] < 5) {
839 err = E_PARSE;
840 } else {
841 set_arma_is_arima(ainfo);
842 }
843 } else {
844 err = E_PARSE;
845 }
846
847 if (!err) {
848 err = check_arma_sep(ainfo, sep1);
849 }
850
851 if (!err) {
852 if (arma_is_arima(ainfo)) {
853 /* check for arima spec */
854 err = check_arima_list(ainfo, dset, opt);
855 } else {
856 /* check for simple arma spec */
857 err = check_arma_list(ainfo, dset, opt);
858 /* catch null model */
859 if (ainfo->nc == 0) {
860 err = E_ARGS;
861 }
862 }
863 }
864
865 #if 0
866 /* catch null model */
867 if (ainfo->nc == 0) {
868 err = E_ARGS;
869 }
870 #endif
871
872 return err;
873 }
874
875 static void
real_arima_difference_series(double * dx,const double * x,int t1,int t2,int * delta,int k)876 real_arima_difference_series (double *dx, const double *x,
877 int t1, int t2, int *delta,
878 int k)
879 {
880 int i, p, t, s = 0;
881
882 for (t=t1; t<=t2; t++) {
883 dx[s] = x[t];
884 for (i=0; i<k && !na(dx[s]); i++) {
885 if (delta[i] != 0) {
886 p = t - i - 1;
887 if (p < 0 || na(x[p])) {
888 dx[s] = NADBL;
889 } else {
890 dx[s] -= delta[i] * x[p];
891 }
892 }
893 }
894 s++;
895 }
896 }
897
transcribe_extra_info(arma_info * ainfo,MODEL * armod)898 static int transcribe_extra_info (arma_info *ainfo, MODEL *armod)
899 {
900 int *ainfo_list = gretl_list_new(9);
901 int err = 0;
902
903 if (ainfo_list == NULL) {
904 armod->errcode = err = E_ALLOC;
905 } else {
906 /* wrap up a bunch of relevant integers */
907 ainfo_list[1] = ainfo->p;
908 ainfo_list[2] = ainfo->q;
909 ainfo_list[3] = ainfo->P;
910 ainfo_list[4] = ainfo->Q;
911 ainfo_list[5] = ainfo->np;
912 ainfo_list[6] = ainfo->nq;
913 ainfo_list[7] = ainfo->d;
914 ainfo_list[8] = ainfo->D;
915 ainfo_list[9] = ainfo->pd;
916 err = gretl_model_set_list_as_data(armod, "ainfo", ainfo_list);
917 }
918
919 if (!err && arma_xdiff(ainfo)) {
920 gretl_model_set_int(armod, "xdiff", 1);
921 }
922
923 return err;
924 }
925
926 #ifndef X12A_CODE
927
928 /* Add to the ainfo struct a full-length series y holding
929 the differenced version of the dependent variable.
930 If the "xdiff" flag is set on ainfo, in addition
931 create a matrix dX holding the differenced regressors;
932 in that case the time-series length of dX depends on
933 the @fullX flag -- if fullX = 0, this equals
934 ainfo->T but if fullX = 0 it equals ainfo->t2 + 1.
935 */
936
arima_difference(arma_info * ainfo,const DATASET * dset,int fullX)937 int arima_difference (arma_info *ainfo,
938 const DATASET *dset,
939 int fullX)
940 {
941 const double *y = dset->Z[ainfo->yno];
942 double *dy = NULL;
943 int *delta = NULL;
944 int s = ainfo->pd;
945 int k, t, t1 = 0;
946 int err = 0;
947
948 #if ARMA_DEBUG
949 fprintf(stderr, "doing arima_difference: d = %d, D = %d\n",
950 ainfo->d, ainfo->D);
951 fprintf(stderr, "ainfo->t1 = %d, ainfo->t2 = %d\n", ainfo->t1,
952 ainfo->t2);
953 #endif
954
955 /* note: dy is a full length series (dset->n) */
956
957 dy = malloc(dset->n * sizeof *dy);
958 if (dy == NULL) {
959 return E_ALLOC;
960 }
961
962 delta = arima_delta_coeffs(ainfo->d, ainfo->D, s);
963 if (delta == NULL) {
964 free(dy);
965 return E_ALLOC;
966 }
967
968 for (t=0; t<dset->n; t++) {
969 dy[t] = NADBL;
970 }
971
972 for (t=0; t<dset->n; t++) {
973 if (na(y[t])) {
974 t1++;
975 } else {
976 break;
977 }
978 }
979
980 t1 += ainfo->d + ainfo->D * s;
981 k = ainfo->d + s * ainfo->D;
982
983 real_arima_difference_series(dy + t1, y, t1, ainfo->t2, delta, k);
984
985 #if ARMA_DEBUG > 1
986 for (t=0; t<dset->n; t++) {
987 fprintf(stderr, "dy[%d] = % 12.7g\n", t, dy[t]);
988 }
989 #endif
990
991 ainfo->y = dy;
992 set_arima_ydiff(ainfo);
993
994 if (arma_xdiff(ainfo)) {
995 /* also difference the ARIMAX regressors */
996 int xt1 = ainfo->t1, xT = ainfo->T;
997
998 if (fullX) {
999 xt1 = 0;
1000 xT = ainfo->t2 + 1;
1001 }
1002
1003 ainfo->dX = gretl_matrix_alloc(xT, ainfo->nexo);
1004
1005 if (ainfo->dX == NULL) {
1006 err = E_ALLOC;
1007 } else {
1008 double *val = ainfo->dX->val;
1009 int i, vi;
1010
1011 for (i=0; i<ainfo->nexo; i++) {
1012 vi = ainfo->xlist[i+1];
1013 real_arima_difference_series(val, dset->Z[vi], xt1,
1014 ainfo->t2, delta, k);
1015 val += xT;
1016 }
1017 }
1018 }
1019
1020 free(delta);
1021
1022 return err;
1023 }
1024
arima_difference_undo(arma_info * ainfo,const DATASET * dset)1025 void arima_difference_undo (arma_info *ainfo, const DATASET *dset)
1026 {
1027 free(ainfo->y);
1028 ainfo->y = (double *) dset->Z[ainfo->yno];
1029
1030 if (ainfo->dX != NULL) {
1031 gretl_matrix_free(ainfo->dX);
1032 ainfo->dX = NULL;
1033 }
1034
1035 unset_arima_ydiff(ainfo);
1036 }
1037
1038 #endif /* X12A_CODE not defined */
1039