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 "libgretl.h"
20 #include "version.h"
21 #include "libset.h"
22 #include "kalman.h"
23 #include "matrix_extra.h"
24 #include "gretl_bfgs.h"
25 #include "arma_priv.h"
26
27 #include "../cephes/libprob.h"
28
29 #define ARMA_DEBUG 0
30 #define ARMA_MDEBUG 0
31 #define SHOW_INIT 0
32
33 #include "arma_common.c"
34
35 #define KALMAN_ALL 999
36
37 static const double *as197_llt_callback (const double *b,
38 int i, void *data);
39
40 static const double *as154_llt_callback (const double *b,
41 int i, void *data);
42
maybe_correct_MA(arma_info * ainfo,double * theta,double * Theta)43 int maybe_correct_MA (arma_info *ainfo,
44 double *theta,
45 double *Theta)
46 {
47 int err = 0;
48
49 if (ainfo->q > 0) {
50 err = flip_poly(theta, ainfo, 0, 0);
51 }
52 if (!err && ainfo->Q > 0) {
53 err = flip_poly(Theta, ainfo, 0, 1);
54 }
55
56 return err;
57 }
58
59 /*
60 Given an ARMA process $A(L)B(L) y_t = C(L)D(L) \epsilon_t$, finds the
61 roots of the four polynomials -- or just two polynomials if seasonal
62 AR and MA effects, B(L) and D(L) are not present -- and attaches
63 this information to the ARMA model.
64
65 pmod: MODEL pointer to which the roots info should be attached.
66
67 ainfo: gives various pieces of information on the ARMA model,
68 including seasonal and non-seasonal AR and MA orders.
69
70 coeff: ifc + p + q + P + Q vector of coefficients (if an intercept
71 is present it is element 0 and is ignored)
72
73 returns: zero on success, non-zero on failure
74 */
75
arma_model_add_roots(MODEL * pmod,arma_info * ainfo,const double * coeff)76 int arma_model_add_roots (MODEL *pmod, arma_info *ainfo,
77 const double *coeff)
78 {
79 const double *phi = coeff + ainfo->ifc;
80 const double *Phi = phi + ainfo->np;
81 const double *theta = Phi + ainfo->P;
82 const double *Theta = theta + ainfo->nq;
83 int nr = ainfo->p + ainfo->P + ainfo->q + ainfo->Q;
84 int pmax, qmax, lmax;
85 double *temp = NULL, *tmp2 = NULL;
86 cmplx *rptr, *roots = NULL;
87 int i, k, cerr = 0;
88
89 pmax = (ainfo->p > ainfo->P)? ainfo->p : ainfo->P;
90 qmax = (ainfo->q > ainfo->Q)? ainfo->q : ainfo->Q;
91 lmax = (pmax > qmax)? pmax : qmax;
92
93 if (pmax == 0 && qmax == 0) {
94 return 0;
95 }
96
97 temp = malloc((lmax + 1) * sizeof *temp);
98 tmp2 = malloc((lmax + 1) * sizeof *tmp2);
99 roots = malloc(nr * sizeof *roots);
100
101 if (temp == NULL || tmp2 == NULL || roots == NULL) {
102 free(temp);
103 free(tmp2);
104 free(roots);
105 return E_ALLOC;
106 }
107
108 temp[0] = 1.0;
109 rptr = roots;
110
111 if (ainfo->p > 0) {
112 /* A(L), non-seasonal */
113 k = 0;
114 for (i=0; i<ainfo->p; i++) {
115 if (AR_included(ainfo, i)) {
116 temp[i+1] = -phi[k++];
117 } else {
118 temp[i+1] = 0;
119 }
120 }
121 cerr = polrt(temp, tmp2, ainfo->p, rptr);
122 rptr += ainfo->p;
123 }
124
125 if (!cerr && ainfo->P > 0) {
126 /* B(L), seasonal */
127 for (i=0; i<ainfo->P; i++) {
128 temp[i+1] = -Phi[i];
129 }
130 cerr = polrt(temp, tmp2, ainfo->P, rptr);
131 rptr += ainfo->P;
132 }
133
134 if (!cerr && ainfo->q > 0) {
135 /* C(L), non-seasonal */
136 k = 0;
137 for (i=0; i<ainfo->q; i++) {
138 if (MA_included(ainfo, i)) {
139 temp[i+1] = theta[k++];
140 } else {
141 temp[i+1] = 0;
142 }
143 }
144 cerr = polrt(temp, tmp2, ainfo->q, rptr);
145 rptr += ainfo->q;
146 }
147
148 if (!cerr && ainfo->Q > 0) {
149 /* D(L), seasonal */
150 for (i=0; i<ainfo->Q; i++) {
151 temp[i+1] = Theta[i];
152 }
153 cerr = polrt(temp, tmp2, ainfo->Q, rptr);
154 }
155
156 free(temp);
157 free(tmp2);
158
159 if (cerr) {
160 free(roots);
161 } else {
162 gretl_model_set_data(pmod, "roots", roots, GRETL_TYPE_CMPLX_ARRAY,
163 nr * sizeof *roots);
164 }
165
166 return 0;
167 }
168
169 /* below: exact ML using Kalman filter apparatus */
170
171 typedef struct kalman_helper_ khelper;
172
173 struct kalman_helper_ {
174 gretl_matrix_block *B;
175 gretl_matrix *S;
176 gretl_matrix *P;
177 gretl_matrix *F;
178 gretl_matrix *A;
179 gretl_matrix *H;
180 gretl_matrix *Q;
181 gretl_matrix *E;
182 gretl_matrix *Svar;
183
184 gretl_matrix *Svar2;
185 gretl_matrix *vQ;
186
187 gretl_matrix *F_; /* used only for ARIMA via levels */
188 gretl_matrix *Q_; /* ditto */
189 gretl_matrix *P_; /* ditto */
190
191 arma_info *kainfo;
192 };
193
kalman_helper_free(khelper * kh)194 static void kalman_helper_free (khelper *kh)
195 {
196 if (kh != NULL) {
197 gretl_matrix_block_destroy(kh->B);
198 gretl_matrix_free(kh->Svar2);
199 gretl_matrix_free(kh->vQ);
200 gretl_matrix_free(kh->F_);
201 gretl_matrix_free(kh->Q_);
202 gretl_matrix_free(kh->P_);
203 free(kh);
204 }
205 }
206
kalman_helper_new(arma_info * ainfo,int r,int k)207 static khelper *kalman_helper_new (arma_info *ainfo,
208 int r, int k)
209 {
210 khelper *kh;
211 int r0, r2;
212 int err = 0;
213
214 kh = malloc(sizeof *kh);
215 if (kh == NULL) {
216 return NULL;
217 }
218
219 r0 = ainfo->r0;
220 r2 = r0 * r0;
221
222 kh->Svar2 = kh->vQ = NULL;
223 kh->F_ = kh->Q_ = kh->P_ = NULL;
224
225 kh->B = gretl_matrix_block_new(&kh->S, r, 1,
226 &kh->P, r, r,
227 &kh->F, r, r,
228 &kh->A, k, 1,
229 &kh->H, r, 1,
230 &kh->Q, r, r,
231 &kh->E, ainfo->fullT, 1,
232 &kh->Svar, r2, r2,
233 NULL);
234
235 if (kh->B == NULL) {
236 err = E_ALLOC;
237 } else if (arma_using_vech(ainfo)) {
238 int m = r0 * (r0 + 1) / 2;
239
240 kh->Svar2 = gretl_matrix_alloc(m, m);
241 kh->vQ = gretl_column_vector_alloc(m);
242 if (kh->Svar2 == NULL || kh->vQ == NULL) {
243 err = E_ALLOC;
244 }
245 } else {
246 kh->vQ = gretl_column_vector_alloc(r2);
247 if (kh->vQ == NULL) {
248 err = E_ALLOC;
249 }
250 }
251
252 if (!err && arima_levels(ainfo)) {
253 kh->F_ = gretl_matrix_alloc(r0, r0);
254 kh->Q_ = gretl_matrix_alloc(r0, r0);
255 kh->P_ = gretl_matrix_alloc(r0, r0);
256 if (kh->F_ == NULL || kh->Q_ == NULL || kh->P_ == NULL) {
257 err = E_ALLOC;
258 }
259 }
260
261 if (err) {
262 kalman_helper_free(kh);
263 kh = NULL;
264 } else {
265 kh->kainfo = ainfo;
266 }
267
268 return kh;
269 }
270
271 /* Get the dimension of the state-space representation: note
272 that this is augmented if we're estimating an ARIMA
273 model using the levels formulation in order to handle
274 missing values -- see Harvey and Pierse, "Estimating
275 Missing Observations in Economic Time Series", JASA 1984.
276 */
277
ainfo_get_state_size(arma_info * ainfo)278 static int ainfo_get_state_size (arma_info *ainfo)
279 {
280 int plen = ainfo->p + ainfo->pd * ainfo->P;
281 int qlen = ainfo->q + ainfo->pd * ainfo->Q;
282 int r = (plen > qlen + 1)? plen : qlen + 1;
283
284 ainfo->r0 = r;
285
286 if (arima_levels(ainfo)) {
287 r += ainfo->d + ainfo->pd * ainfo->D;
288 }
289
290 return r;
291 }
292
allocate_ac_mc(arma_info * ainfo)293 static int allocate_ac_mc (arma_info *ainfo)
294 {
295 int m = (ainfo->P > 0) + (ainfo->Q > 0);
296 int err = 0;
297
298 if (m > 0) {
299 double *ac = NULL, *mc = NULL;
300 int n, i = 0;
301
302 ainfo->aux = doubles_array_new(m, 0);
303 if (ainfo->aux == NULL) {
304 return E_ALLOC;
305 }
306
307 if (ainfo->P > 0) {
308 n = 1 + ainfo->p + ainfo->pd * ainfo->P;
309 ac = malloc(n * sizeof *ac);
310 if (ac == NULL) {
311 err = E_ALLOC;
312 } else {
313 ainfo->aux[i++] = ac;
314 }
315 }
316
317 if (!err && ainfo->Q > 0) {
318 n = 1 + ainfo->q + ainfo->pd * ainfo->Q;
319 mc = malloc(n * sizeof *mc);
320 if (mc == NULL) {
321 err = E_ALLOC;
322 } else {
323 ainfo->aux[i++] = mc;
324 }
325 }
326
327 if (err) {
328 doubles_array_free(ainfo->aux, m);
329 } else {
330 ainfo->n_aux = m;
331 }
332 }
333
334 return err;
335 }
336
write_big_phi(const double * phi,const double * Phi,arma_info * ainfo,gretl_matrix * F)337 static void write_big_phi (const double *phi,
338 const double *Phi,
339 arma_info *ainfo,
340 gretl_matrix *F)
341 {
342 int pmax = ainfo->p + ainfo->pd * ainfo->P;
343 double *ac = ainfo->aux[0];
344 double x, y;
345 int i, j, k, ii;
346
347 for (i=0; i<=pmax; i++) {
348 ac[i] = 0.0;
349 }
350
351 for (j=-1; j<ainfo->P; j++) {
352 x = (j < 0)? -1 : Phi[j];
353 k = 0.0;
354 for (i=-1; i<ainfo->p; i++) {
355 if (i < 0) {
356 y = -1;
357 } else if (AR_included(ainfo, i)) {
358 y = phi[k++];
359 } else {
360 y = 0.0;
361 }
362 ii = (j+1) * ainfo->pd + (i+1);
363 ac[ii] -= x * y;
364 }
365 }
366
367 for (i=0; i<pmax; i++) {
368 gretl_matrix_set(F, 0, i, ac[i+1]);
369 }
370 }
371
write_big_theta(const double * theta,const double * Theta,arma_info * ainfo,gretl_matrix * H,gretl_matrix * F)372 static void write_big_theta (const double *theta,
373 const double *Theta,
374 arma_info *ainfo,
375 gretl_matrix *H,
376 gretl_matrix *F)
377 {
378 int qmax = ainfo->q + ainfo->pd * ainfo->Q;
379 int i = (ainfo->P > 0)? 1 : 0;
380 double *mc = ainfo->aux[i];
381 double x, y;
382 int j, k, ii;
383
384 for (i=0; i<=qmax; i++) {
385 mc[i] = 0.0;
386 }
387
388 for (j=-1; j<ainfo->Q; j++) {
389 x = (j < 0)? 1 : Theta[j];
390 k = 0;
391 for (i=-1; i<ainfo->q; i++) {
392 if (i < 0) {
393 y = 1;
394 } else if (MA_included(ainfo, i)) {
395 y = theta[k++];
396 } else {
397 y = 0.0;
398 }
399 ii = (j+1) * ainfo->pd + (i+1);
400 mc[ii] += x * y;
401 }
402 }
403
404 for (i=1; i<=qmax; i++) {
405 if (H != NULL) {
406 H->val[i] = mc[i];
407 } else {
408 gretl_matrix_set(F, ainfo->r0, i, mc[i]);
409 }
410 }
411 }
412
condense_row(gretl_matrix * targ,const gretl_matrix * src,int targrow,int srcrow,int n)413 static void condense_row (gretl_matrix *targ,
414 const gretl_matrix *src,
415 int targrow, int srcrow,
416 int n)
417 {
418 double x;
419 int i, j, k, g;
420 int targcol = 0;
421
422 for (j=0; j<n; j++) {
423 for (i=j; i<n; i++) {
424 k = j * n + i;
425 g = (k % n) * n + k / n;
426 x = gretl_matrix_get(src, srcrow, k);
427 if (g != k) {
428 x += gretl_matrix_get(src, srcrow, g);
429 }
430 gretl_matrix_set(targ, targrow, targcol++, x);
431 }
432 }
433 }
434
condense_state_vcv(gretl_matrix * targ,const gretl_matrix * src,int n)435 static void condense_state_vcv (gretl_matrix *targ,
436 const gretl_matrix *src,
437 int n)
438 {
439 int posr = 0, posc = 0;
440 int i, j;
441
442 for (i=0; i<n; i++) {
443 for (j=0; j<n; j++) {
444 if (j >= i) {
445 condense_row(targ, src, posr++, posc, n);
446 }
447 posc++;
448 }
449 }
450 }
451
kalman_matrices_init(arma_info * ainfo,khelper * kh,const double * y)452 static int kalman_matrices_init (arma_info *ainfo,
453 khelper *kh,
454 const double *y)
455 {
456 int r0 = ainfo->r0;
457 int r = kh->F->rows;
458
459 gretl_matrix_zero(kh->A);
460 gretl_matrix_zero(kh->S);
461 gretl_matrix_zero(kh->P);
462 gretl_matrix_zero(kh->F);
463 gretl_matrix_inscribe_I(kh->F, 1, 0, r0 - 1);
464
465 gretl_matrix_zero(kh->Q);
466 gretl_matrix_set(kh->Q, 0, 0, 1.0);
467
468 gretl_matrix_zero(kh->H);
469 gretl_vector_set(kh->H, 0, 1.0);
470
471 if (arima_levels(ainfo)) {
472 /* write additional constant elements of F, H and S */
473 int d = ainfo->d, D = ainfo->D;
474 int s = ainfo->pd;
475 int i, k = d + s * D;
476 int *c = arima_delta_coeffs(d, D, s);
477 double y0;
478
479 if (c == NULL) {
480 return E_ALLOC;
481 }
482 for (i=0; i<k; i++) {
483 gretl_matrix_set(kh->F, r0, r0 + i, c[i]);
484 }
485 gretl_matrix_set(kh->F, r0, 0, 1.0);
486 if (r - r0 > 1) {
487 gretl_matrix_inscribe_I(kh->F, r0 + 1, r0, k - 1);
488 }
489 for (i=0; i<k; i++) {
490 gretl_vector_set(kh->H, r0 + i, c[i]);
491 /* lagged data */
492 y0 = y[ainfo->t1 - 1 - i];
493 if (ainfo->yscale != 1.0 && !na(y0)) {
494 y0 -= ainfo->yshift;
495 y0 *= ainfo->yscale;
496 }
497 gretl_vector_set(kh->S, r0 + i, y0);
498 }
499 free(c);
500
501 #if ARMA_DEBUG
502 gretl_matrix_print(kh->S, "S0 (arima via levels)");
503 #endif
504 /* initialize the plain-arma "shadow" matrices */
505 gretl_matrix_zero(kh->F_);
506 gretl_matrix_inscribe_I(kh->F_, 1, 0, r0 - 1);
507 gretl_matrix_zero(kh->Q_);
508 gretl_matrix_set(kh->Q_, 0, 0, 1.0);
509 gretl_matrix_zero(kh->P_);
510 } else if (ainfo->np == 0 && ainfo->P == 0) {
511 /* initialize P to identity matrix */
512 gretl_matrix_inscribe_I(kh->P, 0, 0, kh->P->rows);
513 }
514
515 return 0;
516 }
517
518 #define PRINT_P_INFO 0
519
520 #if PRINT_P_INFO
521
print_P_info(const gretl_matrix * m,arma_info * ainfo)522 static void print_P_info (const gretl_matrix *m,
523 arma_info *ainfo)
524 {
525 double x, x0 = m->val[0];
526 int i, j, id = (x0 == 1.0)? 2 : 1;
527
528 for (j=0; j<m->cols && id > 0; j++) {
529 for (i=0; i<m->rows; i++) {
530 x = gretl_matrix_get(m, i, j);
531 if (i == j && x != x0) {
532 id = 0;
533 } else if (i != j && x != 0.0) {
534 id = 0;
535 }
536 }
537 }
538
539 fprintf(stderr, "%d,%d,%d,%d,%d,%d,%d,%d,%d\n",
540 id, ainfo->np, ainfo->d, ainfo->nq,
541 ainfo->P, ainfo->D, ainfo->Q,
542 ainfo->ifc, ainfo->nexo);
543 }
544
545 #endif
546
write_kalman_matrices(khelper * kh,const double * b,int idx)547 static int write_kalman_matrices (khelper *kh,
548 const double *b,
549 int idx)
550 {
551 arma_info *ainfo = kh->kainfo;
552 const double *phi = b + ainfo->ifc;
553 const double *Phi = phi + ainfo->np;
554 const double *theta = Phi + ainfo->P;
555 const double *Theta = theta + ainfo->nq;
556 const double *beta = Theta + ainfo->Q;
557 double mu = (ainfo->ifc)? b[0] : 0.0;
558 int rewrite_A = 0;
559 int rewrite_F = 0;
560 int rewrite_H = 0;
561 int i, k, err = 0;
562
563 if (idx == KALMAN_ALL) {
564 rewrite_A = rewrite_F = rewrite_H = 1;
565 } else {
566 /* called in context of calculating score, for OPG matrix */
567 int pmax = ainfo->ifc + ainfo->np + ainfo->P;
568 int tmax = pmax + ainfo->nq + ainfo->Q;
569
570 if (ainfo->ifc && idx == 0) {
571 rewrite_A = 1;
572 } else if (idx >= ainfo->ifc && idx < pmax) {
573 rewrite_F = 1;
574 } else if (idx >= ainfo->ifc && idx < tmax) {
575 rewrite_H = 1;
576 } else {
577 rewrite_A = 1;
578 }
579 }
580
581 /* revise for pure MA model */
582 if (ainfo->np == 0 && ainfo->P == 0 && !arima_levels(ainfo)) {
583 rewrite_F = 0;
584 }
585 /* and for case of no constant or other regressors */
586 if (ainfo->ifc == 0 && ainfo->nexo == 0) {
587 rewrite_A = 0;
588 }
589
590 #if ARMA_MDEBUG
591 fprintf(stderr, "write_kalman_matrices: rewrites: A=%d, F=%d, H=%d\n",
592 rewrite_A, rewrite_F, rewrite_H);
593 # if ARMA_MDEBUG > 1
594 fprintf(stderr, "\n*** write_kalman_matrices: before\n");
595 gretl_matrix_print(kh->A, "A");
596 gretl_matrix_print(kh->F, "F");
597 gretl_matrix_print(kh->H, "H");
598 gretl_matrix_print(kh->P, "P");
599 # endif
600 #endif
601
602 /* See Hamilton, Time Series Analysis, ch 13, p. 375 */
603
604 if (rewrite_A) {
605 /* const and coeffs on exogenous vars */
606 gretl_vector_set(kh->A, 0, mu);
607 for (i=0; i<ainfo->nexo; i++) {
608 gretl_vector_set(kh->A, i + 1, beta[i]);
609 }
610 }
611
612 if (rewrite_H) {
613 /* form the H vector using theta and/or Theta */
614 if (ainfo->Q > 0) {
615 write_big_theta(theta, Theta, ainfo, kh->H, NULL);
616 } else {
617 k = 0;
618 for (i=0; i<ainfo->q; i++) {
619 if (MA_included(ainfo, i)) {
620 gretl_vector_set(kh->H, i+1, theta[k++]);
621 } else {
622 gretl_vector_set(kh->H, i+1, 0.0);
623 }
624 }
625 }
626 }
627
628 if (rewrite_F) {
629 /* form the F matrix using phi and/or Phi */
630 gretl_matrix *F = (kh->F_ != NULL)? kh->F_ : kh->F;
631 gretl_matrix *Q = (kh->Q_ != NULL)? kh->Q_ : kh->Q;
632 gretl_matrix *P = (kh->P_ != NULL)? kh->P_ : kh->P;
633
634 if (ainfo->P > 0) {
635 write_big_phi(phi, Phi, ainfo, F);
636 } else {
637 k = 0;
638 for (i=0; i<ainfo->p; i++) {
639 if (AR_included(ainfo, i)) {
640 gretl_matrix_set(F, 0, i, phi[k++]);
641 } else {
642 gretl_matrix_set(F, 0, i, 0.0);
643 }
644 }
645 }
646
647 if (arima_levels(ainfo)) {
648 /* the full F matrix incorporates \theta */
649 if (ainfo->Q > 0) {
650 write_big_theta(theta, Theta, ainfo, NULL, kh->F);
651 } else {
652 k = 0;
653 for (i=0; i<ainfo->q; i++) {
654 if (MA_included(ainfo, i)) {
655 gretl_matrix_set(kh->F, ainfo->r0, i+1, theta[k++]);
656 } else {
657 gretl_matrix_set(kh->F, ainfo->r0, i+1, 0.0);
658 }
659 }
660 }
661 }
662
663 /* form $P_{1|0}$ (MSE) matrix, as per Hamilton, ch 13, p. 378. */
664
665 gretl_matrix_kronecker_product(F, F, kh->Svar);
666 gretl_matrix_I_minus(kh->Svar);
667 if (arma_using_vech(ainfo)) {
668 condense_state_vcv(kh->Svar2, kh->Svar, gretl_matrix_rows(F));
669 gretl_matrix_vectorize_h(kh->vQ, Q);
670 err = gretl_LU_solve(kh->Svar2, kh->vQ);
671 if (!err) {
672 gretl_matrix_unvectorize_h(P, kh->vQ);
673 }
674 } else {
675 gretl_matrix_vectorize(kh->vQ, Q);
676 err = gretl_LU_solve(kh->Svar, kh->vQ);
677 if (!err) {
678 gretl_matrix_unvectorize(P, kh->vQ);
679 }
680 }
681 # if PRINT_P_INFO
682 print_P_info(P, ainfo);
683 # endif
684 }
685
686 if (arima_levels(ainfo)) {
687 /* complete the job on F, Q, P */
688 gretl_matrix_inscribe_matrix(kh->F, kh->F_, 0, 0, GRETL_MOD_NONE);
689 gretl_matrix_inscribe_matrix(kh->Q, kh->Q_, 0, 0, GRETL_MOD_NONE);
690 gretl_matrix_inscribe_matrix(kh->P, kh->P_, 0, 0, GRETL_MOD_NONE);
691 }
692
693 #if ARMA_MDEBUG
694 fprintf(stderr, "\n*** after\n");
695 gretl_matrix_print(kh->A, "A");
696 gretl_matrix_print(kh->F, "F");
697 gretl_matrix_print(kh->H, "H");
698 gretl_matrix_print(kh->P, "P");
699 #endif
700
701 return err;
702 }
703
rewrite_kalman_matrices(kalman * K,const double * b,int i)704 static int rewrite_kalman_matrices (kalman *K, const double *b, int i)
705 {
706 khelper *kh = (khelper *) kalman_get_data(K);
707 int err = write_kalman_matrices(kh, b, i);
708
709 if (!err) {
710 kalman_set_initial_state_vector(K, kh->S);
711 kalman_set_initial_MSE_matrix(K, kh->P);
712 }
713
714 return err;
715 }
716
kalman_arma_llt_callback(const double * b,int i,void * data)717 static const double *kalman_arma_llt_callback (const double *b, int i,
718 void *data)
719 {
720 kalman *K = (kalman *) data;
721 khelper *kh = kalman_get_data(K);
722 int err;
723
724 rewrite_kalman_matrices(K, b, i);
725 err = kalman_forecast(K, NULL);
726
727 #if ARMA_DEBUG
728 fprintf(stderr, "kalman_arma_llt: kalman f'cast gave "
729 "err = %d, ll = %#.12g\n", err, kalman_get_loglik(K));
730 #endif
731
732 return (err)? NULL : kh->E->val;
733 }
734
735 /* add covariance matrix and standard errors based on Outer Product of
736 Gradient
737 */
738
arma_OPG_vcv(MODEL * pmod,void * data,int algo,double * b,double s2,int k,int T,PRN * prn)739 static int arma_OPG_vcv (MODEL *pmod, void *data, int algo,
740 double *b, double s2,
741 int k, int T,
742 PRN *prn)
743 {
744 gretl_matrix *G = NULL;
745 gretl_matrix *V = NULL;
746 int err = 0;
747
748 if (algo == 154) {
749 G = numerical_score_matrix(b, T, k, as154_llt_callback,
750 data, &err);
751 } else if (algo == 197) {
752 G = numerical_score_matrix(b, T, k, as197_llt_callback,
753 data, &err);
754 } else {
755 G = numerical_score_matrix(b, T, k, kalman_arma_llt_callback,
756 data, &err);
757 }
758
759 if (!err) {
760 V = gretl_matrix_XTX_new(G);
761 if (V == NULL) {
762 err = E_ALLOC;
763 }
764 }
765
766 if (!err) {
767 double rcond = gretl_symmetric_matrix_rcond(V, &err);
768
769 if (!err && rcond < 1.0E-10) {
770 pprintf(prn, "OPG: rcond = %g; will try Hessian\n", rcond);
771 err = 1;
772 }
773 }
774
775 if (!err) {
776 err = gretl_invert_symmetric_matrix(V);
777 }
778
779 if (!err) {
780 gretl_matrix_multiply_by_scalar(V, s2);
781 err = gretl_model_write_vcv(pmod, V);
782 }
783
784 gretl_matrix_free(G);
785 gretl_matrix_free(V);
786
787 return err;
788 }
789
arma_QML_vcv(MODEL * pmod,gretl_matrix * H,void * data,int algo,double * b,double s2,int k,int T,PRN * prn)790 static int arma_QML_vcv (MODEL *pmod, gretl_matrix *H,
791 void *data, int algo,
792 double *b, double s2, int k, int T,
793 PRN *prn)
794 {
795 gretl_matrix *G;
796 int err = 0;
797
798 if (algo == 154) {
799 G = numerical_score_matrix(b, T, k, as154_llt_callback,
800 data, &err);
801 } else if (algo == 197) {
802 G = numerical_score_matrix(b, T, k, as197_llt_callback,
803 data, &err);
804 } else {
805 G = numerical_score_matrix(b, T, k, kalman_arma_llt_callback,
806 data, &err);
807 }
808
809 if (!err) {
810 gretl_matrix_divide_by_scalar(G, sqrt(s2));
811 err = gretl_model_add_QML_vcv(pmod, ARMA, H, G,
812 NULL, OPT_NONE, NULL);
813 }
814
815 gretl_matrix_free(G);
816
817 return err;
818 }
819
820 #if ARMA_DEBUG
821
debug_print_theta(const double * theta,const double * Theta,arma_info * ainfo)822 static void debug_print_theta (const double *theta,
823 const double *Theta,
824 arma_info *ainfo)
825 {
826 int i, k = 0;
827
828 fprintf(stderr, "kalman_arma_ll():\n");
829
830 for (i=0; i<ainfo->q; i++) {
831 if (MA_included(ainfo, i)) {
832 fprintf(stderr, "theta[%d] = %#.12g\n", i+1, theta[k++]);
833 }
834 }
835
836 for (i=0; i<ainfo->Q; i++) {
837 fprintf(stderr, "Theta[%d] = %#.12g\n", i, Theta[i]);
838 }
839 }
840
841 #endif
842
843 static int kalman_do_ma_check = 1;
844
kalman_arma_ll(const double * b,void * data)845 static double kalman_arma_ll (const double *b, void *data)
846 {
847 kalman *K = (kalman *) data;
848 khelper *kh = kalman_get_data(K);
849 arma_info *ainfo = kh->kainfo;
850 int offset = ainfo->ifc + ainfo->np + ainfo->P;
851 double *theta = (double *) b + offset;
852 double *Theta = theta + ainfo->nq;
853 double ll = NADBL;
854 int err = 0;
855
856 #if ARMA_DEBUG
857 if (ainfo->q > 0 || ainfo->Q > 0) {
858 debug_print_theta(theta, Theta, ainfo);
859 }
860 #endif
861
862 if (kalman_do_ma_check && maybe_correct_MA(ainfo, theta, Theta)) {
863 pputs(kalman_get_printer(K), _("MA estimate(s) out of bounds\n"));
864 return NADBL;
865 }
866
867 err = rewrite_kalman_matrices(K, b, KALMAN_ALL);
868
869 if (!err) {
870 err = kalman_forecast(K, NULL);
871 ll = kalman_get_loglik(K);
872 #if ARMA_DEBUG > 1
873 fprintf(stderr, "kalman_arma_ll: kalman_forecast gave %d, "
874 "loglik = %#.12g\n", err, ll);
875 #endif
876 }
877
878 return ll;
879 }
880
arima_ydiff_only(arma_info * ainfo)881 static int arima_ydiff_only (arma_info *ainfo)
882 {
883 if ((ainfo->d > 0 || ainfo->D > 0) &&
884 ainfo->nexo > 0 && !arma_xdiff(ainfo)) {
885 return 1;
886 } else {
887 return 0;
888 }
889 }
890
arma_use_opg(gretlopt opt)891 static int arma_use_opg (gretlopt opt)
892 {
893 int ret = 0; /* use of the Hessian is the default */
894
895 if (opt & OPT_G) {
896 ret = 1;
897 } else if (libset_get_int(ARMA_VCV) == ML_OP) {
898 ret = 1;
899 }
900
901 return ret;
902 }
903
904 /* The following is now basically functionless, other
905 than for backward compatibility: it duplicates the
906 model's $uhat array as the $ehat vector, just in
907 case any scripts want it under that name.
908 */
909
arma_add_ehat(MODEL * pmod,arma_info * ainfo,kalman * K,double * b)910 static void arma_add_ehat (MODEL *pmod, arma_info *ainfo,
911 kalman *K, double *b)
912 {
913 khelper *kh = kalman_get_data(K);
914 gretl_matrix *ehat = gretl_matrix_copy(kh->E);
915
916 if (ehat != NULL) {
917 gretl_matrix_set_t1(ehat, pmod->t1);
918 gretl_matrix_set_t2(ehat, pmod->t2);
919 gretl_model_set_matrix_as_data(pmod, "ehat", ehat);
920 }
921 }
922
kalman_arma_finish(MODEL * pmod,arma_info * ainfo,const DATASET * dset,kalman * K,double * b,gretlopt opt,PRN * prn)923 static int kalman_arma_finish (MODEL *pmod,
924 arma_info *ainfo,
925 const DATASET *dset,
926 kalman *K, double *b,
927 gretlopt opt, PRN *prn)
928 {
929 khelper *kh = kalman_get_data(K);
930 int do_opg = arma_use_opg(opt);
931 int kopt, i, t, k = ainfo->nc;
932 int QML = (opt & OPT_R);
933 double s2;
934 int err;
935
936 pmod->t1 = ainfo->t1;
937 pmod->t2 = ainfo->t2;
938 pmod->nobs = ainfo->T;
939 pmod->ncoeff = ainfo->nc;
940 pmod->full_n = dset->n;
941
942 /* in the Kalman case the basic model struct is empty, so we
943 have to allocate for coefficients, residuals and so on
944 */
945 err = gretl_model_allocate_storage(pmod);
946 if (err) {
947 return err;
948 }
949
950 for (i=0; i<k; i++) {
951 pmod->coeff[i] = b[i];
952 }
953
954 i = 0;
955 for (t=pmod->t1; t<=pmod->t2; t++) {
956 pmod->uhat[t] = gretl_vector_get(kh->E, i++);
957 }
958
959 s2 = kalman_get_arma_variance(K);
960 pmod->sigma = sqrt(s2);
961
962 pmod->lnL = kalman_get_loglik(K);
963 kopt = kalman_get_options(K);
964
965 /* rescale if we're using average loglikelihood */
966 if (kopt & KALMAN_AVG_LL) {
967 pmod->lnL *= ainfo->T;
968 }
969
970 if (opt & OPT_E) {
971 arma_add_ehat(pmod, ainfo, K, b);
972 }
973
974 #if ARMA_DEBUG
975 fprintf(stderr, "kalman_arma_finish: doing VCV, method %s\n",
976 (opt & OPT_R)? "QML" : (do_opg)? "OPG" : "Hessian");
977 #endif
978
979 if (!do_opg) {
980 /* base covariance matrix on Hessian (perhaps QML) */
981 gretl_matrix *Hinv;
982 double d = 0.0; /* adjust? */
983
984 kalman_do_ma_check = 0;
985 Hinv = numerical_hessian_inverse(b, ainfo->nc, kalman_arma_ll,
986 K, d, &err);
987 kalman_do_ma_check = 1;
988 if (!err) {
989 if (kopt & KALMAN_AVG_LL) {
990 gretl_matrix_divide_by_scalar(Hinv, ainfo->T);
991 }
992 if (QML) {
993 err = arma_QML_vcv(pmod, Hinv, K, 0, b, s2, k, ainfo->T, prn);
994 } else {
995 err = gretl_model_write_vcv(pmod, Hinv);
996 if (!err) {
997 gretl_model_set_vcv_info(pmod, VCV_ML, ML_HESSIAN);
998 }
999 }
1000 } else if (!(opt & OPT_H)) {
1001 /* fallback when Hessian not explicitly requested */
1002 err = 0;
1003 do_opg = 1;
1004 gretl_model_set_int(pmod, "hess-error", 1);
1005 }
1006 gretl_matrix_free(Hinv);
1007 }
1008
1009 if (do_opg) {
1010 err = arma_OPG_vcv(pmod, K, 0, b, s2, k, ainfo->T, prn);
1011 if (!err) {
1012 gretl_model_set_vcv_info(pmod, VCV_ML, ML_OP);
1013 pmod->opt |= OPT_G;
1014 }
1015 }
1016
1017 if (!err) {
1018 write_arma_model_stats(pmod, ainfo, dset);
1019 arma_model_add_roots(pmod, ainfo, b);
1020 gretl_model_set_int(pmod, "arma_flags", ARMA_EXACT);
1021 if (arma_lbfgs(ainfo)) {
1022 pmod->opt |= OPT_L;
1023 }
1024 if (arima_ydiff_only(ainfo)) {
1025 pmod->opt |= OPT_Y;
1026 }
1027 }
1028
1029 #if 0 /* not ready yet: after estimation of ARMA with missing
1030 values, use the Kalman smoother to estimate the
1031 missing data
1032 */
1033 if (!err && arma_missvals(ainfo)) {
1034 gretl_matrix *ys = kalman_arma_smooth(K, &err);
1035 int t;
1036
1037 i = 0;
1038 for (t=ainfo->t1; t<=ainfo->t2; t++) {
1039 if (na(pmod->yhat[t])) {
1040 pmod->yhat[t] = gretl_vector_get(ys, i);
1041 } else {
1042 fprintf(stderr, "%g vs %g\n", pmod->yhat[t], ys->val[i]);
1043 }
1044 i++;
1045 }
1046 gretl_matrix_free(ys);
1047 }
1048 #endif
1049
1050 return err;
1051 }
1052
kalman_rescale_y(gretl_vector * y,arma_info * ainfo)1053 static void kalman_rescale_y (gretl_vector *y, arma_info *ainfo)
1054 {
1055 int i;
1056
1057 #if ARMA_DEBUG
1058 fprintf(stderr, "kalman_rescale_y: multiplying by %g\n",
1059 ainfo->yscale);
1060 #endif
1061
1062 for (i=0; i<y->rows; i++) {
1063 if (!isnan(y->val[i])) {
1064 y->val[i] -= ainfo->yshift;
1065 y->val[i] *= ainfo->yscale;
1066 }
1067 }
1068 }
1069
1070 /* for Kalman: convert from full-length y series to
1071 y vector of length ainfo->T */
1072
form_arma_y_vector(arma_info * ainfo,int * err)1073 static gretl_matrix *form_arma_y_vector (arma_info *ainfo,
1074 int *err)
1075 {
1076 gretl_matrix *yvec;
1077
1078 yvec = gretl_vector_from_series(ainfo->y, ainfo->t1, ainfo->t2);
1079
1080 if (yvec == NULL) {
1081 *err = E_ALLOC;
1082 } else {
1083 if (ainfo->yscale != 1.0) {
1084 kalman_rescale_y(yvec, ainfo);
1085 }
1086 #if ARMA_DEBUG
1087 gretl_matrix_print(yvec, "arma y vector");
1088 #endif
1089 }
1090
1091 return yvec;
1092 }
1093
form_arma_X_matrix(arma_info * ainfo,const DATASET * dset,int * err)1094 static gretl_matrix *form_arma_X_matrix (arma_info *ainfo,
1095 const DATASET *dset,
1096 int *err)
1097 {
1098 gretl_matrix *X;
1099 int missop;
1100
1101 #if ARMA_DEBUG
1102 printlist(ainfo->xlist, "ainfo->xlist (exog vars)");
1103 #endif
1104
1105 if (arma_na_ok(ainfo)) {
1106 missop = M_MISSING_OK;
1107 } else {
1108 missop = M_MISSING_ERROR;
1109 }
1110
1111 X = gretl_matrix_data_subset(ainfo->xlist, dset,
1112 ainfo->t1, ainfo->t2,
1113 missop, err);
1114
1115 #if ARMA_DEBUG
1116 gretl_matrix_print(X, "X");
1117 #endif
1118
1119 return X;
1120 }
1121
kalman_undo_y_scaling(arma_info * ainfo,gretl_matrix * y,double * b,kalman * K)1122 static int kalman_undo_y_scaling (arma_info *ainfo,
1123 gretl_matrix *y, double *b,
1124 kalman *K)
1125 {
1126 double *beta = b + ainfo->ifc + ainfo->np + ainfo->P +
1127 ainfo->nq + ainfo->Q;
1128 int i, t, T = ainfo->t2 - ainfo->t1 + 1;
1129 int err = 0;
1130
1131 if (ainfo->ifc) {
1132 b[0] /= ainfo->yscale;
1133 b[0] += ainfo->yshift;
1134 }
1135
1136 for (i=0; i<ainfo->nexo; i++) {
1137 beta[i] /= ainfo->yscale;
1138 }
1139
1140 i = ainfo->t1;
1141 for (t=0; t<T; t++) {
1142 y->val[t] /= ainfo->yscale;
1143 y->val[t] += ainfo->yshift;
1144 }
1145
1146 if (na(kalman_arma_ll(b, K))) {
1147 err = 1;
1148 }
1149
1150 return err;
1151 }
1152
free_arma_X_matrix(arma_info * ainfo,gretl_matrix * X)1153 static void free_arma_X_matrix (arma_info *ainfo, gretl_matrix *X)
1154 {
1155 if (X == ainfo->dX) {
1156 gretl_matrix_free(ainfo->dX);
1157 ainfo->dX = NULL;
1158 } else {
1159 gretl_matrix_free(X);
1160 }
1161 }
1162
kalman_arma(const double * coeff,const DATASET * dset,arma_info * ainfo,MODEL * pmod,gretlopt opt)1163 static int kalman_arma (const double *coeff,
1164 const DATASET *dset,
1165 arma_info *ainfo,
1166 MODEL *pmod,
1167 gretlopt opt)
1168 {
1169 kalman *K = NULL;
1170 khelper *kh = NULL;
1171 gretl_matrix *y = NULL;
1172 gretl_matrix *X = NULL;
1173 int r, k = 1 + ainfo->nexo; /* number of exog vars plus space for const */
1174 int use_newton = 0;
1175 double *b;
1176 int err = 0;
1177
1178 b = copyvec(coeff, ainfo->nc);
1179 if (b == NULL) {
1180 return E_ALLOC;
1181 }
1182
1183 #if ARMA_DEBUG
1184 fputs("# kalman_arma: initial coefficients:\n", stderr);
1185 fprintf(stderr, "%d 1\n", ainfo->nc);
1186 int i;
1187 for (i=0; i<ainfo->nc; i++) {
1188 fprintf(stderr, "%.15g\n", b[i]);
1189 }
1190 #endif
1191
1192 y = form_arma_y_vector(ainfo, &err);
1193
1194 if (!err && ainfo->nexo > 0) {
1195 if (ainfo->dX != NULL) {
1196 X = ainfo->dX;
1197 } else {
1198 X = form_arma_X_matrix(ainfo, dset, &err);
1199 }
1200 }
1201
1202 if (!err) {
1203 err = allocate_ac_mc(ainfo);
1204 }
1205
1206 if (err) {
1207 goto bailout;
1208 }
1209
1210 r = ainfo_get_state_size(ainfo);
1211
1212 /* when should we use vech apparatus? */
1213 if (r > 4) {
1214 set_arma_use_vech(ainfo);
1215 }
1216
1217 kh = kalman_helper_new(ainfo, r, k);
1218 if (kh == NULL) {
1219 err = E_ALLOC;
1220 goto bailout;
1221 }
1222
1223 kalman_matrices_init(ainfo, kh, dset->Z[ainfo->yno]);
1224
1225 #if ARMA_DEBUG
1226 fprintf(stderr, "ready to estimate: ainfo specs:\n"
1227 "p=%d, P=%d, q=%d, Q=%d, ifc=%d, nexo=%d, t1=%d, t2=%d\n",
1228 ainfo->p, ainfo->P, ainfo->q, ainfo->Q, ainfo->ifc,
1229 ainfo->nexo, ainfo->t1, ainfo->t2);
1230 fprintf(stderr, "Kalman dims: r = %d, k = %d, T = %d, ncoeff=%d\n",
1231 r, k, ainfo->T, ainfo->nc);
1232 #endif
1233
1234 K = kalman_new(kh->S, kh->P, kh->F, kh->A, kh->H, kh->Q,
1235 NULL, y, X, NULL, kh->E, &err);
1236
1237 if (err) {
1238 fprintf(stderr, "kalman_new(): err = %d\n", err);
1239 } else {
1240 double toler;
1241 int maxit;
1242 int avg_ll;
1243
1244 kalman_attach_printer(K, ainfo->prn);
1245 kalman_attach_data(K, kh);
1246
1247 if (r > 3 && !arima_levels(ainfo)) {
1248 kalman_set_nonshift(K, 1);
1249 } else {
1250 kalman_set_nonshift(K, r);
1251 }
1252
1253 avg_ll = arma_avg_ll(ainfo);
1254 use_newton = libset_get_int(GRETL_OPTIM) == OPTIM_NEWTON;
1255
1256 if (avg_ll) {
1257 kalman_set_options(K, KALMAN_ARMA_LL | KALMAN_AVG_LL);
1258 } else {
1259 kalman_set_options(K, KALMAN_ARMA_LL);
1260 }
1261
1262 BFGS_defaults(&maxit, &toler, ARMA);
1263
1264 if (use_newton) {
1265 double crittol = 1.0e-7;
1266 double gradtol = 1.0e-7;
1267
1268 err = newton_raphson_max(b, ainfo->nc, maxit,
1269 crittol, gradtol, &ainfo->fncount,
1270 C_LOGLIK, kalman_arma_ll,
1271 NULL, NULL, K, opt,
1272 ainfo->prn);
1273 } else {
1274 int save_lbfgs = libset_get_bool(USE_LBFGS);
1275
1276 if (save_lbfgs) {
1277 ainfo->pflags |= ARMA_LBFGS;
1278 } else if (opt & OPT_L) {
1279 libset_set_bool(USE_LBFGS, 1);
1280 ainfo->pflags |= ARMA_LBFGS;
1281 }
1282
1283 err = BFGS_max(b, ainfo->nc, maxit, toler,
1284 &ainfo->fncount, &ainfo->grcount,
1285 kalman_arma_ll, C_LOGLIK,
1286 NULL, K, NULL, opt | OPT_A,
1287 ainfo->prn);
1288
1289 if (save_lbfgs == 0 && (opt & OPT_L)) {
1290 libset_set_bool(USE_LBFGS, 0);
1291 }
1292 }
1293
1294 if (err) {
1295 fprintf(stderr, "kalman_arma: optimizer returned %d\n", err);
1296 }
1297 }
1298
1299 #if ARMA_DEBUG
1300 fprintf(stderr, "undo_scaling? yscale = %g\n", ainfo->yscale);
1301 #endif
1302
1303 if (!err && ainfo->yscale != 1.0) {
1304 kalman_undo_y_scaling(ainfo, y, b, K);
1305 }
1306
1307 if (!err) {
1308 if (use_newton) {
1309 gretl_model_set_int(pmod, "iters", ainfo->fncount);
1310 } else {
1311 gretl_model_set_int(pmod, "fncount", ainfo->fncount);
1312 gretl_model_set_int(pmod, "grcount", ainfo->grcount);
1313 }
1314 err = kalman_arma_finish(pmod, ainfo, dset, K, b,
1315 opt, ainfo->prn);
1316 }
1317
1318 bailout:
1319
1320 if (err) {
1321 pmod->errcode = err;
1322 }
1323
1324 kalman_free(K);
1325 kalman_helper_free(kh);
1326
1327 gretl_matrix_free(y);
1328 free_arma_X_matrix(ainfo, X);
1329 free(b);
1330
1331 return err;
1332 }
1333
1334 /* end of Kalman-specific material */
1335
1336 #define y_missing(y) (na(y) || isnan(y))
1337
1338 /* support for AS 154 and AS 197 */
1339
1340 # include "as197.c"
1341 # include "as154.c"
1342 # include "as_driver.c"
1343
arma_init_message(arma_info * ainfo)1344 static void arma_init_message (arma_info *ainfo)
1345 {
1346 pprintf(ainfo->prn, "\n%s: ", _("ARMA initialization"));
1347
1348 if (ainfo->init == INI_USER) {
1349 pprintf(ainfo->prn, "%s\n\n", _("user-specified values"));
1350 } else if (ainfo->init == INI_HR) {
1351 pprintf(ainfo->prn, "%s\n\n", _("Hannan-Rissanen method"));
1352 } else if (ainfo->init == INI_SMALL) {
1353 pprintf(ainfo->prn, "%s\n\n", _("small MA values"));
1354 } else if (ainfo->init == INI_NLS) {
1355 pprintf(ainfo->prn, "%s\n\n", _("using nonlinear AR model"));
1356 } else if (ainfo->init == INI_OLS) {
1357 pprintf(ainfo->prn, "%s\n\n", _("using linear AR model"));
1358 }
1359 }
1360
user_arma_init(double * coeff,arma_info * ainfo)1361 static int user_arma_init (double *coeff, arma_info *ainfo)
1362 {
1363 int i, nc = n_initvals();
1364
1365 if (nc == 0) {
1366 return 0;
1367 } else if (nc < ainfo->nc) {
1368 pprintf(ainfo->prn, "ARMA initialization: need %d coeffs but got %d\n",
1369 ainfo->nc, nc);
1370 return E_DATA;
1371 }
1372
1373 if (arma_exact_ml(ainfo)) {
1374 /* user-specified initializer is handled within BFGSmax */
1375 for (i=0; i<ainfo->nc; i++) {
1376 coeff[i] = 0.0;
1377 }
1378 } else {
1379 gretl_matrix *m = get_initvals();
1380
1381 for (i=0; i<ainfo->nc; i++) {
1382 coeff[i] = gretl_vector_get(m, i);
1383 }
1384 gretl_matrix_free(m);
1385 }
1386
1387 ainfo->init = INI_USER;
1388
1389 return 0;
1390 }
1391
1392 /* Should we try Hannan-Rissanen initialization of ARMA
1393 coefficients? */
1394
prefer_hr_init(arma_info * ainfo)1395 static int prefer_hr_init (arma_info *ainfo)
1396 {
1397 int ret = 0;
1398
1399 if (ainfo->q > 1 || ainfo->Q > 0) {
1400 ret = 1;
1401 if (arma_xdiff(ainfo)) {
1402 /* don't use for ARIMAX (yet?) */
1403 ret = 0;
1404 } else if (ainfo->T < 100) {
1405 /* unlikely to work well with small sample */
1406 ret = 0;
1407 } else if (ainfo->p > 0 && ainfo->P > 0) {
1408 /* not sure about this: HR catches the MA terms, but NLS
1409 handles the seasonal/non-seasonal AR interactions
1410 better?
1411 */
1412 ret = 0;
1413 } else if ((ainfo->P > 0 && ainfo->p >= ainfo->pd) ||
1414 (ainfo->Q > 0 && ainfo->q >= ainfo->pd)) {
1415 /* overlapping seasonal/non-seasonal orders screw things up */
1416 ret = 0;
1417 } else if (ret && arma_exact_ml(ainfo)) {
1418 /* screen for cases where we'll use NLS */
1419 if (ainfo->P > 0) {
1420 ret = 0;
1421 } else if (ainfo->p + ainfo->P > 0 && ainfo->nexo > 0) {
1422 ret = 0;
1423 } else if (ainfo->Q > 0 && arma_missvals(ainfo)) {
1424 ret = 0;
1425 }
1426 }
1427 }
1428
1429 #if ARMA_DEBUG
1430 fprintf(stderr, "prefer_hr_init? %s\n", ret? "yes" : "no");
1431 #endif
1432
1433 return ret;
1434 }
1435
1436 /* estimate an ARIMA (0,d,0) x (0,D,0) model via OLS */
1437
arima_by_ls(const DATASET * dset,arma_info * ainfo,MODEL * pmod)1438 static int arima_by_ls (const DATASET *dset, arma_info *ainfo,
1439 MODEL *pmod)
1440 {
1441 gretl_matrix *X;
1442 gretl_matrix *b, *u, *V;
1443 double x, s2;
1444 int i, t, k = ainfo->dX->cols;
1445 int err = 0;
1446
1447 if (ainfo->ifc) {
1448 /* the constant will not have been included in ainfo->dX */
1449 X = gretl_matrix_alloc(ainfo->T, k + 1);
1450 if (X == NULL) {
1451 return E_ALLOC;
1452 }
1453 for (i=0; i<=k; i++) {
1454 for (t=0; t<ainfo->T; t++) {
1455 if (i == 0) {
1456 gretl_matrix_set(X, t, i, 1.0);
1457 } else {
1458 x = gretl_matrix_get(ainfo->dX, t, i-1);
1459 gretl_matrix_set(X, t, i, x);
1460 }
1461 }
1462 }
1463 k++;
1464 } else {
1465 X = ainfo->dX;
1466 }
1467
1468 b = gretl_column_vector_alloc(k);
1469 u = gretl_column_vector_alloc(ainfo->T);
1470 V = gretl_matrix_alloc(k, k);
1471
1472 if (b == NULL || u == NULL || V == NULL) {
1473 err = E_ALLOC;
1474 } else {
1475 gretl_vector y;
1476
1477 gretl_matrix_init(&y);
1478 y.rows = ainfo->T;
1479 y.cols = 1;
1480 y.val = ainfo->y + ainfo->t1;
1481 gretl_matrix_set_t1(&y, ainfo->t1);
1482 gretl_matrix_set_t2(&y, ainfo->t2);
1483
1484 err = gretl_matrix_ols(&y, X, b, V, u, &s2);
1485 }
1486
1487 if (!err) {
1488 pmod->ncoeff = k;
1489 pmod->full_n = dset->n;
1490 err = gretl_model_allocate_storage(pmod);
1491 }
1492
1493 if (!err) {
1494 for (i=0; i<k; i++) {
1495 pmod->coeff[i] = b->val[i];
1496 }
1497 for (t=0; t<ainfo->T; t++) {
1498 pmod->uhat[t + ainfo->t1] = u->val[t];
1499 }
1500 err = gretl_model_write_vcv(pmod, V);
1501 }
1502
1503 if (!err) {
1504 pmod->ybar = gretl_mean(ainfo->t1, ainfo->t2, ainfo->y);
1505 pmod->sdy = gretl_stddev(ainfo->t1, ainfo->t2, ainfo->y);
1506 pmod->nobs = ainfo->T;
1507 }
1508
1509 gretl_matrix_free(b);
1510 gretl_matrix_free(u);
1511 gretl_matrix_free(V);
1512
1513 if (X != ainfo->dX) {
1514 gretl_matrix_free(X);
1515 }
1516
1517 return err;
1518 }
1519
1520 /* calculate info criteria for compatibility with ML? */
1521 #define ML_COMPAT 1 /* 2017-03-23 */
1522
arma_via_OLS(arma_info * ainfo,const double * coeff,const DATASET * dset,MODEL * pmod)1523 static int arma_via_OLS (arma_info *ainfo, const double *coeff,
1524 const DATASET *dset, MODEL *pmod)
1525 {
1526 int err = 0;
1527
1528 ainfo->flags |= ARMA_LS;
1529
1530 if (arma_xdiff(ainfo)) {
1531 err = arima_by_ls(dset, ainfo, pmod);
1532 } else {
1533 err = arma_by_ls(coeff, dset, ainfo, pmod);
1534 }
1535
1536 if (!err) {
1537 ArmaFlags f = arma_exact_ml(ainfo) ? ARMA_OLS : ARMA_LS;
1538
1539 pmod->t1 = ainfo->t1;
1540 pmod->t2 = ainfo->t2;
1541 pmod->full_n = dset->n;
1542 write_arma_model_stats(pmod, ainfo, dset);
1543 if (arma_exact_ml(ainfo)) {
1544 #if ML_COMPAT
1545 /* In the case of ainfo->nc == 0 (no coefficients
1546 actually estimated), pmod->ncoeff will be 1, since
1547 we add a dummy constant with value 0. That "1"
1548 will account for the variance estimate so that
1549 addk should be zero. Otherwise we add 1 for the
1550 variance estimate.
1551 */
1552 int addk = ainfo->nc == 0 ? 0 : 1;
1553
1554 mle_criteria(pmod, addk);
1555 #else
1556 ls_criteria(pmod);
1557 #endif
1558 } else {
1559 arma_model_add_roots(pmod, ainfo, pmod->coeff);
1560 }
1561 gretl_model_set_int(pmod, "arma_flags", f);
1562 }
1563
1564 if (!err && pmod->errcode) {
1565 err = pmod->errcode;
1566 }
1567
1568 return err;
1569 }
1570
1571 /* Set flag to indicate differencing of exogenous regressors, in the
1572 case of an ARIMAX model using native exact ML -- unless this is
1573 forbidden by OPT_Y (--y-diff-only). Note that we don't do this
1574 when we're using conditional ML (BHHH).
1575 */
1576
maybe_set_xdiff_flag(arma_info * ainfo,gretlopt opt)1577 static void maybe_set_xdiff_flag (arma_info *ainfo, gretlopt opt)
1578 {
1579 if (arma_exact_ml(ainfo) &&
1580 (ainfo->d > 0 || ainfo->D > 0) &&
1581 ainfo->nexo > 0 && !(opt & OPT_Y)) {
1582 ainfo->pflags |= ARMA_XDIFF;
1583 }
1584 }
1585
1586 /* Respond to OPT_S (--stdx): standardize exogenous regressors */
1587
arma_standardize_x(arma_info * ainfo,DATASET * dset)1588 static int arma_standardize_x (arma_info *ainfo,
1589 DATASET *dset)
1590 {
1591 int orig_v = dset->v;
1592 int err = 0;
1593
1594 ainfo->xstats = gretl_matrix_alloc(ainfo->nexo, 2);
1595 if (ainfo->xstats == NULL) {
1596 return E_ALLOC;
1597 }
1598
1599 err = dataset_add_series(dset, ainfo->nexo);
1600
1601 if (!err) {
1602 double xbar, sdx;
1603 int i, vi, vj, t;
1604
1605 for (i=0; i<ainfo->nexo && !err; i++) {
1606 vi = ainfo->xlist[i+1];
1607 err = gretl_moments(ainfo->t1, ainfo->t2, dset->Z[vi],
1608 NULL, &xbar, &sdx, NULL, NULL, 1);
1609 if (!err) {
1610 vj = orig_v + i;
1611 for (t=0; t<dset->n; t++) {
1612 dset->Z[vj][t] = (dset->Z[vi][t] - xbar) / sdx;
1613 }
1614 /* replace x-ref with standardized version */
1615 ainfo->xlist[i+1] = vj;
1616 /* and record the stats used */
1617 gretl_matrix_set(ainfo->xstats, i, 0, xbar);
1618 gretl_matrix_set(ainfo->xstats, i, 1, sdx);
1619 }
1620 }
1621 }
1622
1623 if (!err) {
1624 set_arma_stdx(ainfo);
1625 }
1626
1627 return err;
1628 }
1629
1630 /* Set flag to allow NAs within the sample range for an
1631 ARMA model using native exact ML.
1632 */
1633
maybe_allow_missvals(arma_info * ainfo)1634 static void maybe_allow_missvals (arma_info *ainfo)
1635 {
1636 if (arma_exact_ml(ainfo)) {
1637 ainfo->pflags |= ARMA_NAOK;
1638 }
1639 }
1640
check_arma_options(gretlopt opt)1641 static int check_arma_options (gretlopt opt)
1642 {
1643 int err;
1644
1645 /* can't specify LBFGS or --robust with conditional ML */
1646 err = options_incompatible_with(opt, OPT_C, OPT_L | OPT_R);
1647
1648 if (!err) {
1649 /* nor more than one of AS, CML or Kalman */
1650 err = incompatible_options(opt, OPT_A | OPT_C | OPT_K);
1651 }
1652
1653 if (!err) {
1654 /* nor --stdx with --kalman */
1655 err = incompatible_options(opt, OPT_S | OPT_K);
1656 }
1657
1658 return err;
1659 }
1660
arma_model(const int * list,const int * pqspec,DATASET * dset,gretlopt opt,PRN * prn)1661 MODEL arma_model (const int *list, const int *pqspec,
1662 DATASET *dset, gretlopt opt, PRN *prn)
1663 {
1664 double *coeff = NULL;
1665 MODEL armod;
1666 arma_info ainfo_s, *ainfo;
1667 int missv = 0, misst = 0;
1668 int orig_v = dset->v;
1669 int err = 0;
1670
1671 ainfo = &ainfo_s;
1672 arma_info_init(ainfo, opt, pqspec, dset);
1673
1674 if (opt & OPT_V) {
1675 ainfo->prn = prn;
1676 }
1677
1678 gretl_model_init(&armod, dset);
1679
1680 err = check_arma_options(opt);
1681
1682 if (!err) {
1683 ainfo->alist = gretl_list_copy(list);
1684 if (ainfo->alist == NULL) {
1685 err = E_ALLOC;
1686 }
1687 }
1688
1689 if (!err) {
1690 err = arma_check_list(ainfo, dset, opt);
1691 }
1692
1693 if (!err) {
1694 /* calculate maximum lag */
1695 maybe_set_xdiff_flag(ainfo, opt);
1696 calc_max_lag(ainfo);
1697 }
1698
1699 if (!err) {
1700 /* adjust sample range if need be */
1701 maybe_allow_missvals(ainfo);
1702 err = arma_adjust_sample(ainfo, dset, &missv, &misst);
1703 if (err) {
1704 if (missv > 0 && misst > 0) {
1705 gretl_errmsg_sprintf(_("Missing value encountered for "
1706 "variable %d, obs %d"), missv, misst);
1707 }
1708 } else if (missv > 0) {
1709 set_arma_missvals(ainfo);
1710 }
1711 }
1712
1713 if (!err && ainfo->nexo > 0 && arma_exact_ml(ainfo)) {
1714 /* FIXME check conditionality more rigorously */
1715 if ((opt & OPT_S) && !arma_xdiff(ainfo)) {
1716 /* --stdx */
1717 err = arma_standardize_x(ainfo, dset);
1718 }
1719 }
1720
1721 if (!err) {
1722 /* allocate initial coefficient vector */
1723 coeff = malloc(ainfo->nc * sizeof *coeff);
1724 if (coeff == NULL) {
1725 err = E_ALLOC;
1726 }
1727 }
1728
1729 if (!err) {
1730 /* organize the dependent variable */
1731 ainfo->y = (double *) dset->Z[ainfo->yno];
1732 if (ainfo->d > 0 || ainfo->D > 0) {
1733 if (arma_missvals(ainfo)) {
1734 /* for now: insist on native Kalman, since only it
1735 handles the levels formulation of ARIMA */
1736 opt &= ~OPT_A;
1737 opt |= OPT_K;
1738 set_arima_levels(ainfo);
1739 } else {
1740 /* note: this replaces ainfo->y */
1741 err = arima_difference(ainfo, dset, 0);
1742 }
1743 }
1744 }
1745
1746 if (err) {
1747 goto bailout;
1748 }
1749
1750 if (ainfo->p == 0 && ainfo->P == 0 &&
1751 ainfo->q == 0 && ainfo->Q == 0 &&
1752 arma_exact_ml(ainfo) && !arma_missvals(ainfo)) {
1753 /* pure "I" model, no NAs: OLS provides the MLE */
1754 err = arma_via_OLS(ainfo, NULL, dset, &armod);
1755 goto bailout; /* estimation handled */
1756 }
1757
1758 /* start initialization of the coefficients */
1759
1760 /* first see if the user specified some values */
1761 err = user_arma_init(coeff, ainfo);
1762 if (err) {
1763 goto bailout;
1764 }
1765
1766 if (!arma_exact_ml(ainfo) && ainfo->q == 0 && ainfo->Q == 0) {
1767 /* for a pure AR model, the conditional MLE is least
1768 squares (OLS or NLS); in the NLS case a user-specified
1769 initializer may be useful, if present
1770 */
1771 const double *b = ainfo->init ? coeff : NULL;
1772
1773 err = arma_via_OLS(ainfo, b, dset, &armod);
1774 goto bailout; /* estimation handled */
1775 }
1776
1777 /* see if it may be helpful to scale the dependent
1778 variable, if we're doing exact ML */
1779 if (arma_exact_ml(ainfo) && ainfo->ifc && ainfo->init != INI_USER) {
1780 maybe_set_yscale(ainfo);
1781 #if SHOW_INIT
1782 fprintf(stderr, "yscale = %g\n", ainfo->yscale);
1783 #endif
1784 }
1785
1786 /* try Hannan-Rissanen init, if suitable */
1787 if (!ainfo->init && prefer_hr_init(ainfo)) {
1788 hr_arma_init(coeff, dset, ainfo);
1789 #if SHOW_INIT
1790 fprintf(stderr, "HR init (%d %d): %s\n", ainfo->p, ainfo->q,
1791 ainfo->init ? "success" : "fail");
1792 #endif
1793 }
1794
1795 /* initialize via AR model by OLS or NLS, adding minimal
1796 MA coefficients if needed: this is the fallback if
1797 Hannan-Rissanen fails, but also the default if the
1798 conditions of applicability of H-R are not met
1799 */
1800 if (!err && !ainfo->init) {
1801 err = ar_arma_init(coeff, dset, ainfo, &armod, opt);
1802 #if SHOW_INIT
1803 fprintf(stderr, "AR init: err = %d\n", err);
1804 #endif
1805 }
1806
1807 if (ainfo->prn != NULL && ainfo->init) {
1808 arma_init_message(ainfo);
1809 }
1810
1811 if (!err) {
1812 clear_model_xpx(&armod);
1813 if (arma_exact_ml(ainfo)) {
1814 if (opt & OPT_K) {
1815 err = kalman_arma(coeff, dset, ainfo, &armod, opt);
1816 } else {
1817 err = as_arma(coeff, dset, ainfo, &armod, opt);
1818 }
1819 } else {
1820 err = bhhh_arma(coeff, dset, ainfo, &armod, opt);
1821 }
1822 }
1823
1824 bailout:
1825
1826 if (err && !armod.errcode) {
1827 armod.errcode = err;
1828 }
1829
1830 if (!err) {
1831 transcribe_extra_info(ainfo, &armod);
1832 }
1833
1834 if (!armod.errcode) {
1835 gretl_model_smpl_init(&armod, dset);
1836 }
1837
1838 free(coeff);
1839 arma_info_cleanup(ainfo);
1840
1841 if (dset->v > orig_v) {
1842 dataset_drop_last_variables(dset, dset->v - orig_v);
1843 }
1844
1845 return armod;
1846 }
1847