1 /*
2 * gretl -- Gnu Regression, Econometrics and Time-series Library
3 * Copyright (C) 2001 Allin Cottrell and Riccardo "Jack" Lucchetti
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, either version 3 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program. If not, see <http://www.gnu.org/licenses/>.
17 *
18 */
19
20 #define FULL_XML_HEADERS 1
21
22 #include "libgretl.h"
23 #include "uservar.h"
24 #include "gretl_func.h"
25 #include "gretl_xml.h"
26 #include "matrix_extra.h"
27 #include "libset.h"
28 #include "kalman.h"
29
30 /**
31 * SECTION:kalman
32 * @short_description: The Kalman filter
33 * @title: Kalman
34 * @include: gretl/libgretl.h, gretl/kalman.h
35 *
36 */
37
38 #define KDEBUG 0
39
40 /* notation (quasi-Hamilton):
41
42 State: S_{t+1} = F*S_t + v_t, E(v_t*v_t') = Q
43
44 Obs: y_t = A'*x_t + H'*S_t + w_t, E(w_t*w_t') = R
45
46 */
47
48 typedef struct crossinfo_ crossinfo;
49
50 struct crossinfo_ {
51 gretl_matrix *BB; /* BB': r x r */
52 gretl_matrix *CC; /* CC': n x n */
53 gretl_matrix *BC; /* BC': r x n */
54 };
55
56 typedef struct stepinfo_ stepinfo;
57
58 struct stepinfo_ {
59 gretl_matrix *F; /* T x (r * r) */
60 gretl_matrix *H; /* T x (r * n) */
61 };
62
63 struct kalman_ {
64 int flags; /* for recording any options */
65 int fnlevel; /* level of function execution */
66
67 int r; /* rows of S = number of elements in state */
68 int n; /* columns of y = number of observables */
69 int k; /* columns of A = number of exogenous vars in obs eqn */
70 int p; /* length of combined disturbance vector */
71 int T; /* rows of y = number of observations */
72 int okT; /* T - number of missing observations */
73 int t; /* current time step, when filtering */
74
75 int ifc; /* boolean: obs equation includes an implicit constant? */
76
77 double SSRw; /* \sum_{t=1}^T e_t^{\prime} V_t^{-1} e_t */
78 double sumldet; /* \sum_{t=1}^T ln |V_t| */
79 double loglik; /* log-likelihood */
80 double s2; /* = SSRw / k */
81
82 int nonshift; /* When F is a companion matrix (e.g. in arma), the
83 number of rows of F that do something other than
84 shifting down the elements of the state vector
85 */
86
87 /* continuously updated matrices */
88 gretl_matrix *S0; /* r x 1: state vector, before updating */
89 gretl_matrix *S1; /* r x 1: state vector, after updating */
90 gretl_matrix *P0; /* r x r: MSE matrix, before updating */
91 gretl_matrix *P1; /* r x r: MSE matrix, after updating */
92 gretl_matrix *e; /* n x 1: one-step forecast error(s), time t */
93
94 /* input data matrices: note that the order matters for various
95 functions, including matrix_is_varying()
96 */
97 gretl_matrix *F; /* r x r: state transition matrix */
98 gretl_matrix *A; /* k x n: coeffs on exogenous vars, obs eqn */
99 gretl_matrix *H; /* r x n: coeffs on state variables, obs eqn */
100 gretl_matrix *Q; /* r x r: contemp covariance matrix, state eqn */
101 gretl_matrix *R; /* n x n: contemp covariance matrix, obs eqn */
102 gretl_matrix *mu; /* r x 1: constant term in state transition */
103 gretl_matrix *y; /* T x n: dependent variable vector (or matrix) */
104 gretl_matrix *x; /* T x k: independent variables matrix */
105 gretl_matrix *Sini; /* r x 1: S_{1|0} */
106 gretl_matrix *Pini; /* r x r: P_{1|0} */
107
108 /* user inputs for cross-correlated disturbances */
109 gretl_matrix *B; /* r x p: BB' = Q */
110 gretl_matrix *C; /* n x p: CC' = R */
111
112 /* apparatus for registering time-variation of matrices */
113 char *matcall;
114 char *varying;
115
116 /* optional matrices for recording extra info */
117 gretl_matrix *LL; /* T x 1: loglikelihood, all time-steps */
118
119 /* optional run-time export matrices */
120 gretl_matrix *E; /* T x n: forecast errors, all time-steps */
121 gretl_matrix *V; /* T x nn: MSE for observables, all time-steps */
122 gretl_matrix *S; /* T x r: state vector, all time-steps */
123 gretl_matrix *P; /* T x nr: MSE for state, all time-steps */
124 gretl_matrix *K; /* T x rn: gain matrix, all time-steps */
125 gretl_matrix *U; /* T x ??: smoothed disturbances */
126 gretl_matrix *Vsd; /* Variance of smoothed disturbance */
127
128 /* structure needed only for cross-correlated case */
129 crossinfo *cross;
130
131 /* structure needed only when smoothing in the time-varying case */
132 stepinfo *step;
133
134 /* workspace matrices */
135 gretl_matrix_block *Blk; /* holder for the following */
136 gretl_matrix *PH;
137 gretl_matrix *HPH;
138 gretl_matrix *FPH;
139 gretl_matrix *Vt;
140 gretl_matrix *Ve;
141 gretl_matrix *PHV;
142 gretl_matrix *Ax;
143 gretl_matrix *Kt;
144 gretl_matrix *Tmpnn;
145 gretl_matrix *Tmprr;
146 gretl_matrix *Tmprr_2a;
147 gretl_matrix *Tmprr_2b;
148 gretl_matrix *Tmpr1;
149
150 gretl_bundle *b; /* the bundle of which this struct is a member */
151 void *data; /* handle for attaching additional info */
152 PRN *prn; /* verbose printer */
153 };
154
155 /* max number of time-varying matrices: F, A, H, Q, R, mu */
156 #define K_N_MATCALLS 6
157
158 #define arma_ll(K) (K->flags & KALMAN_ARMA_LL)
159
160 #define set_kalman_running(K) (K->flags |= KALMAN_FORWARD)
161 #define set_kalman_stopped(K) (K->flags &= ~KALMAN_FORWARD)
162 #define kalman_is_running(K) (K->flags & KALMAN_FORWARD)
163 #define kalman_simulating(K) (K->flags & KALMAN_SIM)
164 #define kalman_checking(K) (K->flags & KALMAN_CHECK)
165 #define kalman_xcorr(K) (K->flags & KALMAN_CROSS)
166 #define kalman_ssfsim(K) (K->flags & KALMAN_SSFSIM)
167
168 #define filter_is_varying(K) (K->matcall != NULL)
169
170 static const char *kalman_matrix_name (int sym);
171 static int kalman_revise_variance (kalman *K);
172 static int check_for_matrix_updates (kalman *K, ufunc *uf);
173
174 /* symbolic identifiers for input matrices: note that potentially
175 time-varying matrices must appear first in the enumeration, and
176 the order must match the order of the "input data matrices" in
177 the Kalman struct (above).
178 */
179
180 enum {
181 K_F = 0,
182 K_A,
183 K_H,
184 K_Q,
185 K_R,
186 K_m,
187 K_y,
188 K_x,
189 K_S,
190 K_P,
191 K_MMAX /* sentinel */
192 };
193
free_stepinfo(kalman * K)194 void free_stepinfo (kalman *K)
195 {
196 if (K->step != NULL) {
197 gretl_matrix_free(K->step->F);
198 gretl_matrix_free(K->step->H);
199 free(K->step);
200 K->step = NULL;
201 }
202 }
203
free_crossinfo(crossinfo * c)204 void free_crossinfo (crossinfo *c)
205 {
206 gretl_matrix_free(c->BB);
207 gretl_matrix_free(c->CC);
208 gretl_matrix_free(c->BC);
209
210 free(c);
211 }
212
kalman_free(kalman * K)213 void kalman_free (kalman *K)
214 {
215 if (K == NULL) {
216 return;
217 }
218
219 gretl_matrix_free(K->S0);
220 gretl_matrix_free(K->P0);
221 gretl_matrix_free(K->S1);
222 gretl_matrix_free(K->P1);
223 gretl_matrix_free(K->e);
224 gretl_matrix_free(K->LL);
225
226 gretl_matrix_block_destroy(K->Blk);
227
228 if (K->flags & KALMAN_BUNDLE) {
229 gretl_matrix **mptr[] = {
230 &K->F, &K->A, &K->H, &K->Q, &K->R,
231 &K->mu, &K->y, &K->x, &K->Sini, &K->Pini
232 };
233 int i;
234
235 if (kalman_xcorr(K)) {
236 mptr[3] = &K->B;
237 mptr[4] = &K->C;
238 }
239
240 for (i=0; i<K_MMAX; i++) {
241 gretl_matrix_free(*mptr[i]);
242 }
243
244 /* we also own these "export" matrices */
245 gretl_matrix_free(K->E);
246 gretl_matrix_free(K->V);
247 gretl_matrix_free(K->S);
248 gretl_matrix_free(K->P);
249 gretl_matrix_free(K->K);
250 gretl_matrix_free(K->U);
251 gretl_matrix_free(K->Vsd);
252 }
253
254 /* time-variation info */
255 free(K->matcall);
256 free(K->varying);
257
258 if (K->cross != NULL) {
259 free_crossinfo(K->cross);
260 }
261
262 if (K->step != NULL) {
263 free_stepinfo(K);
264 }
265
266 free(K);
267 }
268
kalman_new_empty(int flags)269 static kalman *kalman_new_empty (int flags)
270 {
271 kalman *K = malloc(sizeof *K);
272
273 if (K != NULL) {
274 K->Sini = K->Pini = NULL;
275 K->S0 = K->S1 = NULL;
276 K->P0 = K->P1 = NULL;
277 K->LL = NULL;
278 K->e = NULL;
279 K->Blk = NULL;
280 K->F = K->A = K->H = NULL;
281 K->Q = K->R = NULL;
282 K->B = K->C = NULL;
283 K->E = K->V = K->S = K->P = K->K = NULL;
284 K->y = K->x = NULL;
285 K->mu = NULL;
286 K->U = NULL;
287 K->Vsd = NULL;
288 K->matcall = NULL;
289 K->varying = NULL;
290 K->cross = NULL;
291 K->step = NULL;
292 K->flags = flags;
293 K->fnlevel = 0;
294 K->t = 0;
295 K->prn = NULL;
296 K->data = NULL;
297 K->b = NULL;
298 }
299
300 return K;
301 }
302
303 #define kappa 1.0e7
304
diffuse_Pini(kalman * K)305 static void diffuse_Pini (kalman *K)
306 {
307 int i;
308
309 gretl_matrix_zero(K->P0);
310
311 for (i=0; i<K->r; i++) {
312 gretl_matrix_set(K->P0, i, i, kappa);
313 }
314 }
315
F_out_of_bounds(kalman * K)316 static int F_out_of_bounds (kalman *K)
317 {
318 gretl_matrix *evals;
319 double r, c, x;
320 int i, err = 0;
321
322 evals = gretl_general_matrix_eigenvals(K->F, &err);
323
324 for (i=0; i<evals->rows && !err; i++) {
325 r = gretl_matrix_get(evals, i, 0);
326 c = gretl_matrix_get(evals, i, 1);
327 x = sqrt(r*r + c*c);
328 if (x >= 1.0) {
329 fprintf(stderr, "F: modulus of eigenvalue %d = %g\n", i+1, x);
330 err = E_SINGULAR;
331 }
332 }
333
334 gretl_matrix_free(evals);
335
336 return err;
337 }
338
339 /* If the user has not given an initial value for P_{1|0}, compute
340 this automatically as per Hamilton, ch 13, p. 378. This works only
341 if the eigenvalues of K->F lie inside the unit circle. Failing
342 that, or if the --diffuse option is given for the user Kalman
343 filter, we apply a diffuse initialization.
344 */
345
construct_Pini(kalman * K)346 static int construct_Pini (kalman *K)
347 {
348 gretl_matrix *Svar;
349 gretl_matrix *vQ;
350 int r2, err = 0;
351
352 if (K->flags & KALMAN_DIFFUSE) {
353 diffuse_Pini(K);
354 return 0;
355 }
356
357 r2 = K->r * K->r;
358
359 Svar = gretl_matrix_alloc(r2, r2);
360 vQ = gretl_column_vector_alloc(r2);
361
362 if (Svar == NULL || vQ == NULL) {
363 gretl_matrix_free(Svar);
364 gretl_matrix_free(vQ);
365 return E_ALLOC;
366 }
367
368 gretl_matrix_kronecker_product(K->F, K->F, Svar);
369 gretl_matrix_I_minus(Svar);
370 gretl_matrix_vectorize(vQ, K->Q);
371
372 err = gretl_LU_solve(Svar, vQ);
373 if (err) {
374 /* failed: are some of the eigenvalues out of bounds? */
375 err = F_out_of_bounds(K);
376 if (err == E_SINGULAR) {
377 err = 0;
378 diffuse_Pini(K);
379 K->flags |= KALMAN_DIFFUSE;
380 }
381 } else {
382 gretl_matrix_unvectorize(K->P0, vQ);
383 }
384
385 gretl_matrix_free(Svar);
386 gretl_matrix_free(vQ);
387
388 return err;
389 }
390
check_matrix_dims(kalman * K,const gretl_matrix * m,int i)391 static int check_matrix_dims (kalman *K, const gretl_matrix *m, int i)
392 {
393 int r = 0, c = 0, symm = (i == K_Q || i == K_R);
394 int err = 0;
395
396 if (i == K_F || i == K_Q || i == K_P) {
397 r = c = K->r;
398 } else if (i == K_A) {
399 r = K->k;
400 c = K->n;
401 } else if (i == K_H) {
402 r = K->r;
403 c = K->n;
404 } else if (i == K_R) {
405 r = c = K->n;
406 } else if (i == K_S || i == K_m) {
407 r = K->r;
408 c = 1;
409 }
410
411 if (m->rows != r || m->cols != c) {
412 gretl_errmsg_sprintf("kalman: %s is %d x %d, should be %d x %d\n",
413 kalman_matrix_name(i), m->rows, m->cols, r, c);
414 err = E_NONCONF;
415 } else if (symm && !gretl_matrix_is_symmetric(m)) {
416 gretl_errmsg_sprintf("kalman: %s is not symmetric\n",
417 kalman_matrix_name(i));
418 err = E_NONCONF;
419 }
420
421 return err;
422 }
423
424 enum {
425 K_E,
426 K_V,
427 K_BIG_S,
428 K_BIG_P,
429 K_K,
430 K_LL,
431 };
432
maybe_resize_export_matrix(kalman * K,gretl_matrix * m,int i)433 static int maybe_resize_export_matrix (kalman *K, gretl_matrix *m, int i)
434 {
435 int rows = K->T, cols = 0;
436 int err = 0;
437
438 if (i == K_E) {
439 cols = K->n;
440 } else if (i == K_V) {
441 cols = (K->n * K->n + K->n) / 2;
442 } else if (i == K_BIG_S) {
443 cols = K->r;
444 } else if (i == K_BIG_P) {
445 cols = (K->r * K->r + K->r) / 2;
446 } else if (i == K_LL) {
447 cols = 1;
448 } else if (i == K_K) {
449 cols = K->r * K->n;
450 } else {
451 err = E_DATA;
452 }
453
454 if (!err && (m->rows != rows || m->cols != cols)) {
455 err = gretl_matrix_realloc(m, rows, cols);
456 if (!err) {
457 gretl_matrix_zero(m);
458 }
459 }
460
461 return err;
462 }
463
missing_matrix_error(const char * name)464 static int missing_matrix_error (const char *name)
465 {
466 if (name == NULL) {
467 gretl_errmsg_set(_("kalman: a required matrix is missing"));
468 } else {
469 gretl_errmsg_sprintf(_("kalman: required matrix %s is missing"),
470 name);
471 }
472 return E_DATA;
473 }
474
kalman_check_dimensions(kalman * K)475 static int kalman_check_dimensions (kalman *K)
476 {
477 int err = 0;
478
479 if (K->r < 1 || K->n < 1 || K->T < 2) {
480 /* the state and observation vectors must have at least one
481 element, and there must be at least two observations */
482 err = E_DATA;
483 }
484
485 /* F is mandatory, should be r x r */
486 if (!err) {
487 err = check_matrix_dims(K, K->F, K_F);
488 }
489
490 /* H is mandatory, should be r x n */
491 if (!err) {
492 err = check_matrix_dims(K, K->H, K_H);
493 }
494
495 if (K->B != NULL) {
496 /* cross-correlated disturbances */
497 if (K->C == NULL) {
498 err = missing_matrix_error("obsymat");
499 } else if (K->B->rows != K->r || K->B->cols != K->p ||
500 K->C->rows != K->n || K->C->cols != K->p) {
501 err = E_NONCONF;
502 }
503 } else {
504 /* Q is mandatory, should be r x r and symmetric */
505 if (!err) {
506 err = check_matrix_dims(K, K->Q, K_Q);
507 }
508
509 /* R should be n x n and symmetric, if present */
510 if (!err && K->R != NULL) {
511 err = check_matrix_dims(K, K->R, K_R);
512 }
513 }
514
515 /* initial S should be r x 1, if present */
516 if (!err && K->Sini != NULL) {
517 err = check_matrix_dims(K, K->Sini, K_S);
518 }
519
520 /* initial P should be r x r, if present */
521 if (!err && K->Pini != NULL) {
522 err = check_matrix_dims(K, K->Pini, K_P);
523 }
524
525 /* A should be k x n, if present (A->rows defines k; n is defined by y) */
526 if (!err) {
527 if (K->A != NULL) {
528 err = check_matrix_dims(K, K->A, K_A);
529 } else if (K->x != NULL) {
530 /* A is NULL => can't have a non-NULL x */
531 err = E_NONCONF;
532 }
533 }
534
535 /* mu should be r x 1, if present */
536 if (!err && K->mu != NULL) {
537 err = check_matrix_dims(K, K->mu, K_m);
538 }
539
540 if (err) {
541 goto bailout;
542 }
543
544 K->ifc = 0;
545
546 /* x should have T rows to match y; and it should have either k or k - 1
547 columns (the latter case indicating an implicit const) */
548 if (K->x != NULL) {
549 if (K->x->rows < K->T) {
550 fprintf(stderr, "kalman: %s has %d rows, should have %d\n",
551 kalman_matrix_name(K_x), K->x->rows, K->T);
552 return E_NONCONF;
553 } else if (K->x->cols != K->k && K->x->cols != K->k - 1) {
554 fprintf(stderr, "kalman: %s has %d columns, should have %d or %d\n",
555 kalman_matrix_name(K_x), K->x->cols, K->k, K->k - 1);
556 return E_NONCONF;
557 } else if (K->x->cols == K->k - 1) {
558 /* register the implicit const */
559 K->ifc = 1;
560 }
561 } else if (K->k == 1) {
562 /* A has one row but there's no x => implicit const */
563 K->ifc = 1;
564 } else if (K->k > 1) {
565 /* A has more than one row but there's no x => error */
566 return missing_matrix_error("obsxmat");
567 }
568
569 /* Below we have the optional "export" matrices for shipping out
570 results. If these are present but not sized correctly we'll
571 try to fix them up -- but note that they are not used in a
572 simulation run.
573 */
574
575 if (kalman_simulating(K)) {
576 return err;
577 }
578
579 /* E should be T x n */
580 if (K->E != NULL) {
581 err = maybe_resize_export_matrix(K, K->E, K_E);
582 }
583
584 /* V should be T x nn */
585 if (!err && K->V != NULL) {
586 err = maybe_resize_export_matrix(K, K->V, K_V);
587 }
588
589 /* big S should be T x r */
590 if (!err && K->S != NULL) {
591 err = maybe_resize_export_matrix(K, K->S, K_BIG_S);
592 }
593
594 /* big P should be T x nr */
595 if (!err && K->P != NULL) {
596 err = maybe_resize_export_matrix(K, K->P, K_BIG_P);
597 }
598
599 /* LL should be T x 1 */
600 if (!err && K->LL != NULL) {
601 err = maybe_resize_export_matrix(K, K->LL, K_LL);
602 }
603
604 /* K (gain) should be T x (r * n) */
605 if (!err && K->K != NULL) {
606 err = maybe_resize_export_matrix(K, K->K, K_K);
607 }
608
609 bailout:
610
611 if (err) {
612 fprintf(stderr, "kalman_check_dimensions: err = %d\n", err);
613 }
614
615 return err;
616 }
617
618 /* Write the vech of @src into row @t of @targ */
619
load_to_vech(gretl_matrix * targ,const gretl_matrix * src,int n,int t)620 static void load_to_vech (gretl_matrix *targ,
621 const gretl_matrix *src,
622 int n, int t)
623 {
624 int i, j, m = 0;
625 double x;
626
627 for (i=0; i<n; i++) {
628 for (j=i; j<n; j++) {
629 x = gretl_matrix_get(src, i, j);
630 gretl_matrix_set(targ, t, m++, x);
631 }
632 }
633 }
634
635 /* Write the vec of @src into row @t of @targ */
636
load_to_vec(gretl_matrix * targ,const gretl_matrix * src,int t)637 static void load_to_vec (gretl_matrix *targ,
638 const gretl_matrix *src,
639 int t)
640 {
641 int i;
642
643 for (i=0; i<targ->cols; i++) {
644 gretl_matrix_set(targ, t, i, src->val[i]);
645 }
646 }
647
648 /* Write the square root of diagonal of square matrix @src
649 into row @t of @targ, starting at column offset @j
650 */
651
load_to_diag(gretl_matrix * targ,const gretl_matrix * src,int t,int j)652 static void load_to_diag (gretl_matrix *targ,
653 const gretl_matrix *src,
654 int t, int j)
655 {
656 int i, n = gretl_vector_get_length(src);
657 double x;
658
659 for (i=0; i<n; i++) {
660 x = gretl_matrix_get(src, i, i);
661 if (x <= 0.0) {
662 gretl_matrix_set(targ, t, i+j, 0.0);
663 } else {
664 gretl_matrix_set(targ, t, i+j, sqrt(x));
665 }
666 }
667 }
668
669 /* copy from vector @src into row @t of @targ */
670
load_to_row(gretl_matrix * targ,const gretl_vector * src,int t)671 static void load_to_row (gretl_matrix *targ,
672 const gretl_vector *src,
673 int t)
674 {
675 double x;
676 int i;
677
678 for (i=0; i<targ->cols; i++) {
679 x = gretl_vector_get(src, i);
680 gretl_matrix_set(targ, t, i, x);
681 }
682 }
683
684 /* copy from vector @src into row @t of @targ,
685 starting at column offset @j in @targ */
686
load_to_row_offset(gretl_matrix * targ,const gretl_vector * src,int t,int j)687 static void load_to_row_offset (gretl_matrix *targ,
688 const gretl_vector *src,
689 int t, int j)
690 {
691 int i, n = gretl_vector_get_length(src);
692
693 for (i=0; i<n; i++) {
694 gretl_matrix_set(targ, t, i+j, src->val[i]);
695 }
696 }
697
set_row(gretl_matrix * targ,int t,double x)698 static void set_row (gretl_matrix *targ, int t, double x)
699 {
700 int i;
701
702 for (i=0; i<targ->cols; i++) {
703 gretl_matrix_set(targ, t, i, x);
704 }
705 }
706
707 /* checks if row has a 1 in position n and 0s elsewhere */
708
ok_companion_row(const gretl_matrix * F,int cols,int row,int n)709 static int ok_companion_row (const gretl_matrix *F,
710 int cols, int row,
711 int n)
712 {
713 double x;
714 int j;
715
716 for (j=0; j<cols; j++) {
717 x = gretl_matrix_get(F, row, j);
718 if ((j == n && x != 1.0) || (j != n && x != 0.0)) {
719 return 0;
720 }
721 }
722
723 return 1;
724 }
725
count_nonshifts(const gretl_matrix * F)726 static int count_nonshifts (const gretl_matrix *F)
727 {
728 int rmax = gretl_matrix_rows(F);
729 double x;
730 int i, j, n = 0;
731 int ret = rmax;
732
733 /* examine bottom row */
734 for (j=0; j<rmax; j++) {
735 x = gretl_matrix_get(F, rmax - 1, j);
736 if (x != 0.0) {
737 if (n == 0 && j > 0 && x == 1.0) {
738 n = j;
739 } else {
740 return rmax;
741 }
742 }
743 }
744
745 /* at this point n holds the candidate; now check n rows from
746 the (r - n - 1)-th */
747 if (n > 0) {
748 ret = n;
749 j = 0;
750 for (i=rmax-n; i<rmax; i++) {
751 if (!ok_companion_row(F, rmax, i, j++)) {
752 return rmax;
753 }
754 }
755 }
756
757 return ret;
758 }
759
kalman_init(kalman * K)760 static int kalman_init (kalman *K)
761 {
762 int err = 0;
763
764 K->SSRw = NADBL;
765 K->sumldet = NADBL;
766 K->loglik = NADBL;
767 K->s2 = NADBL;
768
769 clear_gretl_matrix_err();
770
771 if (K->Sini != NULL) {
772 K->S0 = gretl_matrix_copy(K->Sini);
773 K->S1 = gretl_matrix_copy(K->Sini);
774 } else {
775 K->S0 = gretl_zero_matrix_new(K->r, 1);
776 K->S1 = gretl_zero_matrix_new(K->r, 1);
777 }
778
779 if (K->Pini != NULL) {
780 K->P0 = gretl_matrix_copy(K->Pini);
781 K->P1 = gretl_matrix_copy(K->Pini);
782 } else {
783 K->P0 = gretl_zero_matrix_new(K->r, K->r);
784 K->P1 = gretl_zero_matrix_new(K->r, K->r);
785 }
786
787 /* forecast error vector, per observation */
788 K->e = gretl_matrix_alloc(K->n, 1);
789
790 err = get_gretl_matrix_err();
791 if (err) {
792 return err;
793 }
794
795 K->nonshift = -1;
796
797 K->Blk = gretl_matrix_block_new(&K->PH, K->r, K->n, /* P*H */
798 &K->FPH, K->r, K->n, /* F*P*H */
799 &K->HPH, K->n, K->n, /* H'*P*H */
800 &K->Vt, K->n, K->n, /* (H'*P*H + R)^{-1} */
801 &K->Ve, K->n, 1, /* (H'*P*H + R)^{-1} * e */
802 &K->PHV, K->r, K->n, /* P*H*V */
803 &K->Ax, K->n, 1, /* A'*x at obs t */
804 &K->Kt, K->r, K->n, /* gain at t */
805 &K->Tmpnn, K->n, K->n,
806 &K->Tmprr, K->r, K->r,
807 &K->Tmprr_2a, K->r, K->r,
808 &K->Tmprr_2b, K->r, K->r,
809 &K->Tmpr1, K->r, 1,
810 NULL);
811
812 if (K->Blk == NULL) {
813 err = E_ALLOC;
814 }
815
816 if (!err && K->Pini == NULL && !(K->flags & KALMAN_USER)) {
817 /* in the "user" case we do this later */
818 err = construct_Pini(K);
819 }
820
821 return err;
822 }
823
kalman_set_dimensions(kalman * K)824 static void kalman_set_dimensions (kalman *K)
825 {
826 K->r = gretl_matrix_rows(K->F); /* F->rows defines r */
827 K->k = gretl_matrix_rows(K->A); /* A->rows defines k */
828 K->n = gretl_matrix_cols(K->y); /* y->cols defines n */
829
830 if (!kalman_simulating(K)) {
831 /* y->rows defines T, except when simulating */
832 K->T = gretl_matrix_rows(K->y);
833 }
834
835 K->okT = K->T;
836
837 /* K->p is non-zero only under cross-correlation; in that case the
838 matrix given as 'Q' in Kalman set-up in fact represents B (as
839 in v_t = B \varepsilon_t) and it must be r x p, where p is the
840 number of elements in the "combined" disturbance vector
841 \varepsilon_t.
842 */
843 K->p = (K->B != NULL)? gretl_matrix_cols(K->B): 0;
844 }
845
846 /**
847 * kalman_new:
848 * @S: r x 1 initial state vector.
849 * @P: r x r initial precision matrix.
850 * @F: r x r state transition matrix.
851 * @A: n x k matrix of coefficients on exogenous variables in the
852 * observation equation.
853 * @H: n x r matrix of coefficients on the state variables in the
854 * observation equation.
855 * @Q: r x r contemporaneous covariance matrix for the errors in the
856 * state equation.
857 * @R: n x n contemporaneous covariance matrix for the errors in the
858 * observation equation (or NULL if this is not applicable).
859 * @y: T x n matrix of dependent variable(s).
860 * @x: T x k matrix of exogenous variable(s). May be NULL if there
861 * are no exogenous variables, or if there's only a constant.
862 * @m: r x 1 vector of constants in the state transition, or NULL.
863 * @E: T x n matrix in which to record forecast errors (or NULL if
864 * this is not required).
865 * @err: location to receive error code.
866 *
867 * Allocates and initializes a Kalman struct, which can subsequently
868 * be used for forecasting with kalman_forecast(). The nomenclature
869 * for the various required matrices is that in Hamilton's Time
870 * Series Analysis (1994, chapter 13), except that "S" is used in
871 * place of Hamilton's \xi for the state vector.
872 *
873 * Returns: pointer to allocated struct, or NULL on failure, in
874 * which case @err will receive a non-zero code.
875 */
876
kalman_new(gretl_matrix * S,gretl_matrix * P,gretl_matrix * F,gretl_matrix * A,gretl_matrix * H,gretl_matrix * Q,gretl_matrix * R,gretl_matrix * y,gretl_matrix * x,gretl_matrix * m,gretl_matrix * E,int * err)877 kalman *kalman_new (gretl_matrix *S, gretl_matrix *P,
878 gretl_matrix *F, gretl_matrix *A,
879 gretl_matrix *H, gretl_matrix *Q,
880 gretl_matrix *R, gretl_matrix *y,
881 gretl_matrix *x, gretl_matrix *m,
882 gretl_matrix *E, int *err)
883 {
884 kalman *K;
885
886 *err = 0;
887
888 if (y == NULL || F == NULL || H == NULL || Q == NULL) {
889 fprintf(stderr, "kalman_new: y=%p, F=%p, H=%p, Q=%p\n",
890 (void *) y, (void *) F, (void *) H, (void *) Q);
891 *err = missing_matrix_error(NULL);
892 return NULL;
893 }
894
895 K = kalman_new_empty(0);
896 if (K == NULL) {
897 *err = E_ALLOC;
898 return NULL;
899 }
900
901 /* use pointers for input matrices, don't copy */
902 K->F = F;
903 K->A = A;
904 K->H = H;
905 K->Q = Q;
906 K->R = R;
907 K->y = y;
908 K->x = x;
909 K->Sini = S;
910 K->Pini = P;
911 K->mu = m;
912
913 /* output, but again use external pointer */
914 K->E = E;
915
916 kalman_set_dimensions(K);
917
918 *err = kalman_check_dimensions(K);
919 if (*err) {
920 free(K);
921 return NULL;
922 }
923
924 *err = kalman_init(K);
925
926 if (*err) {
927 kalman_free(K);
928 K = NULL;
929 } else {
930 gretl_matrix_zero(K->e);
931 }
932
933 return K;
934 }
935
936 /* supports hansl function for creating a named Kalman bundle */
937
kalman_new_minimal(gretl_matrix * M[],int copy[],int nmat,int * err)938 kalman *kalman_new_minimal (gretl_matrix *M[], int copy[],
939 int nmat, int *err)
940 {
941 gretl_matrix **targ[5];
942 kalman *K;
943 int i;
944
945 *err = 0;
946
947 if (M[0] == NULL || M[1] == NULL || M[2] == NULL || M[3] == NULL) {
948 fprintf(stderr, "kalman_new_minimal: nmat=%d, y=%p, H=%p, F=%p, Q=%p\n",
949 nmat, (void *) M[0], (void *) M[1], (void *) M[2], (void *) M[3]);
950 *err = missing_matrix_error(NULL);
951 return NULL;
952 }
953
954 K = kalman_new_empty(KALMAN_USER | KALMAN_BUNDLE);
955 if (K == NULL) {
956 *err = E_ALLOC;
957 return NULL;
958 }
959
960 targ[0] = &K->y;
961 targ[1] = &K->H;
962 targ[2] = &K->F;
963
964 if (nmat == 5) {
965 K->flags |= KALMAN_CROSS;
966 targ[3] = &K->B;
967 targ[4] = &K->C;
968 } else {
969 targ[3] = &K->Q;
970 }
971
972 for (i=0; i<nmat; i++) {
973 if (copy[i]) {
974 *targ[i] = gretl_matrix_copy(M[i]);
975 } else {
976 *targ[i] = M[i];
977 }
978 }
979
980 #if 0
981 gretl_matrix_print(K->y, "K->y");
982 gretl_matrix_print(K->H, "K->H");
983 gretl_matrix_print(K->F, "K->F");
984 #endif
985
986 kalman_set_dimensions(K);
987
988 if (K->p > 0) {
989 *err = kalman_revise_variance(K);
990 }
991
992 if (!*err) {
993 *err = kalman_check_dimensions(K);
994 }
995
996 if (!*err) {
997 *err = kalman_init(K);
998 }
999
1000 if (*err) {
1001 kalman_free(K);
1002 K = NULL;
1003 } else {
1004 gretl_matrix_zero(K->e);
1005 }
1006
1007 return K;
1008 }
1009
matrix_is_varying(kalman * K,int i)1010 static int matrix_is_varying (kalman *K, int i)
1011 {
1012 if (K->matcall != NULL) {
1013 if (K->varying == NULL) {
1014 check_for_matrix_updates(K, NULL);
1015 }
1016 if (K->varying != NULL) {
1017 return K->varying[i];
1018 }
1019 }
1020
1021 return 0;
1022 }
1023
1024 enum {
1025 UPDATE_INIT, /* initialization of matrices */
1026 UPDATE_STEP /* refreshing matrices per time-step */
1027 };
1028
1029 /* After reading 'Q' = B and 'R' = C from the user, either at
1030 (re-)initialization or at a given time-step in the case where either
1031 of these matrices is time-varying, record the user input in K->B
1032 and K->C and form the 'real' Q and R. But note that in the
1033 time-step case it may be that only one of Q, R needs to be treated
1034 in this way (if only one is time-varying, only one will have
1035 been redefined via a function call).
1036 */
1037
kalman_update_crossinfo(kalman * K,int mode)1038 static int kalman_update_crossinfo (kalman *K, int mode)
1039 {
1040 int err = 0;
1041
1042 /* Note that B and C may be needed as such for simulation */
1043
1044 if (mode == UPDATE_INIT || matrix_is_varying(K, K_Q)) {
1045 /* recreate BB using modified B */
1046 err = gretl_matrix_multiply_mod(K->B, GRETL_MOD_NONE,
1047 K->B, GRETL_MOD_TRANSPOSE,
1048 K->cross->BB, GRETL_MOD_NONE);
1049 }
1050
1051 if (!err && (mode == UPDATE_INIT || matrix_is_varying(K, K_R))) {
1052 /* recreate CC using modified C */
1053 err = gretl_matrix_multiply_mod(K->C, GRETL_MOD_NONE,
1054 K->C, GRETL_MOD_TRANSPOSE,
1055 K->cross->CC, GRETL_MOD_NONE);
1056 }
1057
1058 if (!err && (mode == UPDATE_INIT || matrix_is_varying(K, K_Q) ||
1059 matrix_is_varying(K, K_R))) {
1060 /* recreate BC using modified B and/or C */
1061 err = gretl_matrix_multiply_mod(K->B, GRETL_MOD_NONE,
1062 K->C, GRETL_MOD_TRANSPOSE,
1063 K->cross->BC, GRETL_MOD_NONE);
1064 }
1065
1066 /* redundant? */
1067 if (!err && mode == UPDATE_STEP) {
1068 if (matrix_is_varying(K, K_Q)) {
1069 K->Q = K->cross->BB;
1070 }
1071 if (matrix_is_varying(K, K_R)) {
1072 K->R = K->cross->CC;
1073 }
1074 }
1075
1076 return err;
1077 }
1078
kalman_add_crossinfo(kalman * K)1079 static int kalman_add_crossinfo (kalman *K)
1080 {
1081 crossinfo *c = malloc(sizeof *c);
1082
1083 if (c == NULL) {
1084 return E_ALLOC;
1085 }
1086
1087 c->BB = gretl_matrix_alloc(K->r, K->r);
1088 c->CC = gretl_matrix_alloc(K->n, K->n);
1089 c->BC = gretl_matrix_alloc(K->r, K->n);
1090
1091 if (c->BB == NULL || c->CC == NULL || c->BC == NULL) {
1092 free_crossinfo(c);
1093 return E_ALLOC;
1094 }
1095
1096 K->cross = c;
1097
1098 return 0;
1099 }
1100
1101 /* kalman_revise_variance: the user has actually given B in place of Q
1102 and C in place of R; so we have to form Q = BB', R = CC', and BC'
1103 (cross-correlated disturbances).
1104
1105 This function is called in the course of initial set-up of a
1106 filter, and also when the Kalman matrices are being re-checked at
1107 the start of filtering, smoothing or simulation.
1108 */
1109
kalman_revise_variance(kalman * K)1110 static int kalman_revise_variance (kalman *K)
1111 {
1112 int err = 0;
1113
1114 if (K->B == NULL || K->C == NULL) {
1115 return missing_matrix_error("'statevar' or 'obsvar'");
1116 }
1117
1118 #if 0
1119 fprintf(stderr, "kalman_revise_variance\n");
1120 gretl_matrix_print(K->B, "B");
1121 gretl_matrix_print(K->C, "C");
1122 #endif
1123
1124 if (K->cross == NULL) {
1125 /* not allocated yet: this should be the case only
1126 on initial set-up */
1127 err = kalman_add_crossinfo(K);
1128 }
1129
1130 if (!err) {
1131 err = kalman_update_crossinfo(K, UPDATE_INIT);
1132 }
1133
1134 if (!err) {
1135 /* establish convenience pointers */
1136 K->Q = K->cross->BB;
1137 K->R = K->cross->CC;
1138 #if 0
1139 gretl_matrix_print(K->cross->BB, "BB");
1140 gretl_matrix_print(K->cross->CC, "CC");
1141 gretl_matrix_print(K->cross->BC, "BC");
1142 #endif
1143 }
1144
1145 if (err) {
1146 fprintf(stderr, "kalman_revise_variance: err = %d\n", err);
1147 }
1148
1149 return err;
1150 }
1151
matrix_diff(const gretl_matrix * a,const gretl_matrix * b,double tol)1152 static int matrix_diff (const gretl_matrix *a,
1153 const gretl_matrix *b,
1154 double tol)
1155 {
1156 int i, n = a->rows * a->cols;
1157
1158 /* note: we ignore the 0,0 element here */
1159
1160 for (i=1; i<n; i++) {
1161 if (fabs(b->val[i] - a->val[i]) > tol) {
1162 return 1;
1163 }
1164 }
1165
1166 return 0;
1167 }
1168
1169 /* below: if postmult is non-zero, we're post-multiplying by the
1170 transpose of F */
1171
multiply_by_F(kalman * K,const gretl_matrix * A,gretl_matrix * B,int postmult)1172 static int multiply_by_F (kalman *K, const gretl_matrix *A,
1173 gretl_matrix *B, int postmult)
1174 {
1175 int ret = 0;
1176
1177 #if 0 /* uninitialized stuff here? */
1178 gretl_matrix_print(A, "A, in multiply_by_F");
1179 #endif
1180
1181 if (gretl_is_zero_matrix(A)) {
1182 gretl_matrix_zero(B);
1183 } else if (K->nonshift == K->r) {
1184 if (postmult) {
1185 gretl_matrix_multiply_mod(A, GRETL_MOD_NONE,
1186 K->F, GRETL_MOD_TRANSPOSE,
1187 B, GRETL_MOD_NONE);
1188 } else {
1189 gretl_matrix_multiply(K->F, A, B);
1190 }
1191 } else {
1192 gretl_matrix *topF = K->Tmprr_2a;
1193 gretl_matrix *top = K->Tmprr_2b;
1194 int r1 = K->nonshift;
1195 int r2 = K->r - r1;
1196 int i, j, c;
1197 double x;
1198
1199 if (postmult) {
1200 /* if post-multiplying by F', "top" actually means "left" */
1201
1202 topF->rows = K->r;
1203 topF->cols = r1;
1204 c = A->rows;
1205 top->rows = c;
1206 top->cols = r1;
1207
1208 /* copy from F to topF */
1209 for (i=0; i<r1; i++) {
1210 for (j=0; j<K->r; j++) {
1211 x = gretl_matrix_get(K->F, i, j);
1212 gretl_matrix_set(topF, j, i, x);
1213 }
1214 }
1215
1216 gretl_matrix_multiply(A, topF, top);
1217
1218 for (i=0; i<c; i++) {
1219 for (j=0; j<r1; j++) {
1220 x = gretl_matrix_get(top, i, j);
1221 gretl_matrix_set(B, i, j, x);
1222 }
1223 }
1224
1225 for (i=0; i<c; i++) {
1226 for (j=0; j<r2; j++) {
1227 x = gretl_matrix_get(A, i, j);
1228 gretl_matrix_set(B, i, j + r1, x);
1229 }
1230 }
1231
1232 } else {
1233 /* pre-multiplying by F */
1234
1235 topF->rows = r1;
1236 topF->cols = K->r;
1237 c = A->cols;
1238 top->rows = r1;
1239 top->cols = c;
1240
1241 for (i=0; i<r1; i++) {
1242 for (j=0; j<K->r; j++) {
1243 x = gretl_matrix_get(K->F, i, j);
1244 gretl_matrix_set(topF, i, j, x);
1245 }
1246 }
1247
1248 gretl_matrix_multiply(topF, A, top);
1249
1250 for (i=0; i<r1; i++) {
1251 for (j=0; j<c; j++) {
1252 x = gretl_matrix_get(top, i, j);
1253 gretl_matrix_set(B, i, j, x);
1254 }
1255 }
1256
1257 for (i=0; i<r2; i++) {
1258 for (j=0; j<c; j++) {
1259 x = gretl_matrix_get(A, i, j);
1260 gretl_matrix_set(B, i + r1, j, x);
1261 }
1262 }
1263 }
1264 }
1265
1266 return ret;
1267 }
1268
1269 /* Simplified version of the function below, for ARMA:
1270 in this case we have H = (r x 1) and S = (r x 1),
1271 although F and P are (r x r).
1272 */
1273
kalman_arma_iter_1(kalman * K,int missobs)1274 static int kalman_arma_iter_1 (kalman *K, int missobs)
1275 {
1276 double Ve;
1277 int i, err = 0;
1278
1279 /* write F*S into S+ */
1280 err += multiply_by_F(K, K->S0, K->S1, 0);
1281
1282 if (missobs) {
1283 return err;
1284 }
1285
1286 /* form e = y - A'x - H'S */
1287 K->e->val[0] -= K->Ax->val[0];
1288 for (i=0; i<K->r; i++) {
1289 K->e->val[0] -= K->H->val[i] * K->S0->val[i];
1290 }
1291
1292 /* form FPH */
1293 err += multiply_by_F(K, K->PH, K->FPH, 0);
1294
1295 /* form (H'PH + R)^{-1} * (y - Ax - H'S) = "Ve" */
1296 Ve = K->Vt->val[0] * K->e->val[0];
1297
1298 /* form FPH * (H'PH + R)^{-1} * (y - A'x - H'S),
1299 and complete calculation of S+ */
1300 for (i=0; i<K->r; i++) {
1301 K->S1->val[i] += K->FPH->val[i] * Ve;
1302 }
1303
1304 if (!err) {
1305 /* scalar quadratic form */
1306 K->SSRw += Ve * K->e->val[0];
1307 }
1308
1309 /* 2018-04-07: the following (which is now standard)
1310 was confined to the case of KALMAN_ETT being set,
1311 and I believe the formula was wrong, namely
1312 K->e->val[0] = Ve, which amounts dividing e->val[0]
1313 by its estimated variance rather than standard
1314 deviation. AC. */
1315 K->e->val[0] = sqrt(K->Vt->val[0]) * K->e->val[0];
1316
1317 return err;
1318 }
1319
1320 /* Hamilton (1994) equation [13.2.23] page 381, in simplified notation:
1321
1322 S+ = FS + FPH(H'PH + R)^{-1} * (y - A'x - H'S)
1323
1324 "S" is Hamilton's \xi (state vector)
1325 */
1326
kalman_iter_1(kalman * K,int missobs,double * llt)1327 static int kalman_iter_1 (kalman *K, int missobs, double *llt)
1328 {
1329 int err = 0;
1330
1331 /* write F*S into S+ */
1332 err += multiply_by_F(K, K->S0, K->S1, 0);
1333
1334 /* add \mu if present */
1335 if (K->mu != NULL) {
1336 gretl_matrix_add_to(K->S1, K->mu);
1337 }
1338
1339 if (missobs) {
1340 /* the observable is not in fact observed at t */
1341 *llt = 0.0;
1342 if (K->K != NULL) {
1343 set_row(K->K, K->t, 0.0);
1344 }
1345 return err;
1346 }
1347
1348 /* form e = y - A'x - H'S (e is already initialized to y) */
1349 err += gretl_matrix_subtract_from(K->e, K->Ax);
1350 err += gretl_matrix_multiply_mod(K->H, GRETL_MOD_TRANSPOSE,
1351 K->S0, GRETL_MOD_NONE,
1352 K->e, GRETL_MOD_DECREMENT);
1353
1354 if (!err) {
1355 /* contribution to log-likelihood -- see Hamilton
1356 (1994) equation [13.4.1] page 385.
1357 */
1358 double x = gretl_scalar_qform(K->e, K->Vt, &err);
1359
1360 if (err) {
1361 *llt = NADBL;
1362 } else {
1363 *llt -= 0.5 * x;
1364 K->SSRw += x;
1365 }
1366 }
1367
1368 /* form the gain, Kt = (FPH + BC') * (H'PH + R)^{-1} */
1369 err += multiply_by_F(K, K->PH, K->FPH, 0);
1370 if (K->p > 0) {
1371 /* cross-correlated case */
1372 gretl_matrix_add_to(K->FPH, K->cross->BC);
1373 }
1374 err += gretl_matrix_multiply(K->FPH, K->Vt, K->Kt);
1375
1376 /* form K_t * e_t and add to S+ */
1377 err += gretl_matrix_multiply_mod(K->Kt, GRETL_MOD_NONE,
1378 K->e, GRETL_MOD_NONE,
1379 K->S1, GRETL_MOD_CUMULATE);
1380
1381 if (!err && K->K != NULL) {
1382 /* record the gain */
1383 load_to_vec(K->K, K->Kt, K->t);
1384 }
1385
1386 return err;
1387 }
1388
1389 /* Hamilton (1994) equation [13.2.22] page 380, in simplified notation:
1390
1391 P+ = F[P - PH(H'PH + R)^{-1}H'P]F' + Q
1392
1393 Or use the alternative formulation:
1394
1395 P+ = FPF' - KVK + Q where K = gain and V = error variance
1396 */
1397
kalman_iter_2(kalman * K,int missobs)1398 static int kalman_iter_2 (kalman *K, int missobs)
1399 {
1400 int err = 0;
1401
1402 if (K->p == 0 && !missobs) {
1403 /* revise "P0" as P_{t|t} = P - PH(H'PH + R)^{-1}H'P */
1404 err = gretl_matrix_qform(K->PH, GRETL_MOD_NONE, K->Vt,
1405 K->P0, GRETL_MOD_DECREMENT);
1406 }
1407
1408 /* pre-multiply by F, post-multiply by F' */
1409 if (K->nonshift == K->r) {
1410 err += gretl_matrix_qform(K->F, GRETL_MOD_NONE, K->P0,
1411 K->P1, GRETL_MOD_NONE);
1412 } else {
1413 err += multiply_by_F(K, K->P0, K->Tmprr, 0);
1414 err += multiply_by_F(K, K->Tmprr, K->P1, 1);
1415 gretl_matrix_xtr_symmetric(K->P1);
1416 }
1417
1418 if (K->p > 0 && !missobs) {
1419 /* subtract K(H'PH + R)K' */
1420 gretl_matrix_qform(K->Kt, GRETL_MOD_NONE, K->HPH,
1421 K->P1, GRETL_MOD_DECREMENT);
1422 }
1423
1424 /* add Q */
1425 err += gretl_matrix_add_to(K->P1, K->Q);
1426
1427 return err;
1428 }
1429
1430 #if KDEBUG > 1
kalman_print_state(kalman * K)1431 static void kalman_print_state (kalman *K)
1432 {
1433 int j;
1434
1435 /* if (t > 5) return; */
1436
1437 fprintf(stderr, "t = %d:\n", K->t);
1438
1439 for (j=0; j<K->n; j++) {
1440 fprintf(stderr, "y[%d] = %.10g, err[%d] = %.10g\n", j,
1441 gretl_matrix_get(K->y, K->t, j),
1442 j, gretl_vector_get(K->e, j));
1443 }
1444
1445 gretl_matrix_print(K->S0, "K->S0");
1446 gretl_matrix_print(K->P0, "K->P0");
1447 }
1448 #endif
1449
kalman_record_state(kalman * K)1450 static void kalman_record_state (kalman *K)
1451 {
1452 if (K->S != NULL) {
1453 load_to_row(K->S, K->S0, K->t);
1454 }
1455
1456 if (K->P != NULL) {
1457 load_to_vech(K->P, K->P0, K->r, K->t);
1458 }
1459 }
1460
kalman_set_nonshift(kalman * K,int n)1461 void kalman_set_nonshift (kalman *K, int n)
1462 {
1463 K->nonshift = n;
1464 }
1465
kalman_set_options(kalman * K,int opts)1466 void kalman_set_options (kalman *K, int opts)
1467 {
1468 K->flags |= opts;
1469 }
1470
kalman_get_options(kalman * K)1471 int kalman_get_options (kalman *K)
1472 {
1473 return K->flags;
1474 }
1475
1476 /* Read from the appropriate row of x (T x k) and multiply by A' to
1477 form A'x_t. Note: the flag K->ifc is used to indicate that the
1478 observation equation has an implicit constant, with an entry in
1479 the A matrix (the first) but no explicit entry in the x matrix.
1480 */
1481
kalman_set_Ax(kalman * K,int * missobs)1482 static void kalman_set_Ax (kalman *K, int *missobs)
1483 {
1484 double xjt, axi;
1485 int i, j;
1486
1487 for (i=0; i<K->n; i++) {
1488 axi = 0.0;
1489 for (j=0; j<K->k; j++) {
1490 if (K->ifc) {
1491 xjt = (j == 0)? 1.0 : gretl_matrix_get(K->x, K->t, j - 1);
1492 } else {
1493 xjt = gretl_matrix_get(K->x, K->t, j);
1494 }
1495 if (isnan(xjt)) {
1496 #if KDEBUG
1497 fprintf(stderr, "x: got nan at obs %d\n", K->t);
1498 #endif
1499 *missobs += 1;
1500 }
1501 axi += xjt * gretl_matrix_get(K->A, j, i);
1502 }
1503 gretl_vector_set(K->Ax, i, axi);
1504 }
1505 }
1506
1507 /* read from the appropriate row of y (T x n) and transcribe to
1508 the current e (n x 1)
1509 */
1510
kalman_initialize_error(kalman * K,int * missobs)1511 static void kalman_initialize_error (kalman *K, int *missobs)
1512 {
1513 double yti;
1514 int i;
1515
1516 for (i=0; i<K->n; i++) {
1517 yti = gretl_matrix_get(K->y, K->t, i);
1518 K->e->val[i] = yti;
1519 if (isnan(yti)) {
1520 #if KDEBUG
1521 fprintf(stderr, "y: got nan at obs %d\n", K->t);
1522 #endif
1523 *missobs += 1;
1524 }
1525 }
1526 }
1527
1528 /* Given a unified function to update one or more of the
1529 potentially time-varying matrices, try to figure out
1530 which matrix or matrices are actually modified by
1531 this function.
1532 */
1533
check_for_matrix_updates(kalman * K,ufunc * uf)1534 static int check_for_matrix_updates (kalman *K, ufunc *uf)
1535 {
1536 char **lines;
1537 int nlines = 0;
1538
1539 if (K->varying != NULL) {
1540 free(K->varying);
1541 K->varying = NULL;
1542 }
1543
1544 if (uf == NULL) {
1545 uf = get_user_function_by_name(K->matcall);
1546 if (uf == NULL) {
1547 gretl_errmsg_sprintf("Couldn't find function '%s'", K->matcall);
1548 return E_DATA;
1549 }
1550 }
1551
1552 K->varying = calloc(K_N_MATCALLS, 1);
1553
1554 lines = gretl_function_retrieve_code(uf, &nlines);
1555
1556 if (lines != NULL) {
1557 const char *bname = fn_param_name(uf, 0);
1558 char test[VNAMELEN+1];
1559 const char *s;
1560 int n = strlen(bname) + 1;
1561 int i, j;
1562
1563 sprintf(test, "%s.", bname);
1564 for (i=0; i<nlines; i++) {
1565 if (!strncmp(lines[i], test, n)) {
1566 for (j=K_F; j<=K_m; j++) {
1567 s = kalman_matrix_name(j);
1568 if (!strncmp(lines[i] + n, s, strlen(s))) {
1569 fprintf(stderr, "matrix %s is varying\n", s);
1570 K->varying[j] = 1;
1571 break;
1572 }
1573 }
1574 }
1575 }
1576 free(lines);
1577 }
1578
1579 return 0;
1580 }
1581
1582 /* Function to update any time-varying matrices, for use
1583 with a kalman bundle. Bypasses the regular "genr" apparatus,
1584 passing the attached bundle directly to the given user
1585 function after is has been found by name.
1586 */
1587
kalman_update_matrices(kalman * K,PRN * prn)1588 static int kalman_update_matrices (kalman *K, PRN *prn)
1589 {
1590 ufunc *uf;
1591 fncall *fc;
1592 int err = 0;
1593
1594 uf = get_user_function_by_name(K->matcall);
1595
1596 if (uf == NULL) {
1597 gretl_errmsg_sprintf("Couldn't find function '%s'", K->matcall);
1598 return E_DATA;
1599 }
1600
1601 if (K->varying == NULL) {
1602 check_for_matrix_updates(K, uf);
1603 }
1604
1605 fc = fncall_new(uf, 0);
1606 err = push_anon_function_arg(fc, GRETL_TYPE_BUNDLE_REF, K->b);
1607
1608 if (!err) {
1609 err = gretl_function_exec(fc, GRETL_TYPE_NONE, NULL, NULL,
1610 NULL, prn);
1611 }
1612
1613 if (err) {
1614 fprintf(stderr, "kalman_update_matrices: call='%s', err=%d\n",
1615 K->matcall, err);
1616 }
1617
1618 return err;
1619 }
1620
1621 /* If we have any time-varying coefficient matrices, refresh these for
1622 the current time step. This is called on a forward filtering pass.
1623 */
1624
kalman_refresh_matrices(kalman * K,PRN * prn)1625 static int kalman_refresh_matrices (kalman *K, PRN *prn)
1626 {
1627 gretl_matrix **mptr[] = {
1628 &K->F, &K->A, &K->H, &K->Q, &K->R, &K->mu
1629 };
1630 int cross_update = 0;
1631 int i, err = 0;
1632
1633 if (kalman_xcorr(K)) {
1634 mptr[3] = &K->B;
1635 mptr[4] = &K->C;
1636 }
1637
1638 if (K->matcall != NULL) {
1639 err = kalman_update_matrices(K, prn);
1640 }
1641
1642 for (i=0; i<K_N_MATCALLS && !err; i++) {
1643 if (matrix_is_varying(K, i)) {
1644 if (kalman_xcorr(K) && (i == K_Q || i == K_R)) {
1645 /* handle revised B and/or C */
1646 cross_update = 1;
1647 } else {
1648 err = check_matrix_dims(K, *mptr[i], i);
1649 }
1650 if (err) {
1651 fprintf(stderr, "kalman_refresh_matrices: err = %d at t = %d\n",
1652 err, K->t);
1653 }
1654 }
1655 }
1656
1657 if (!err && K->step != NULL) {
1658 /* keep a record of F and/or H at the given time step */
1659 if (K->step->F != NULL) {
1660 load_to_vec(K->step->F, K->F, K->t);
1661 }
1662 if (K->step->H != NULL) {
1663 load_to_vec(K->step->H, K->H, K->t);
1664 }
1665 }
1666
1667 if (!err && cross_update) {
1668 /* cross-correlated case */
1669 err = kalman_update_crossinfo(K, UPDATE_STEP);
1670 }
1671
1672 return err;
1673 }
1674
1675 /* Variant of the above for use when Koopman-smoothing */
1676
ksmooth_refresh_matrices(kalman * K,PRN * prn)1677 static int ksmooth_refresh_matrices (kalman *K, PRN *prn)
1678 {
1679 gretl_matrix **mptr[] = {
1680 &K->Q, &K->R
1681 };
1682 int idx[] = {
1683 K_Q, K_R
1684 };
1685 int cross_update = 0;
1686 int i, ii, err = 0;
1687
1688 if (kalman_xcorr(K)) {
1689 mptr[0] = &K->B;
1690 mptr[1] = &K->C;
1691 }
1692
1693 if (K->matcall != NULL) {
1694 err = kalman_update_matrices(K, prn);
1695 }
1696
1697 for (i=0; i<2 && !err; i++) {
1698 ii = idx[i];
1699 if (matrix_is_varying(K, ii)) {
1700 if (kalman_xcorr(K) && (ii == K_Q || ii == K_R)) {
1701 /* handle revised B and/or C */
1702 cross_update = 1;
1703 } else {
1704 err = check_matrix_dims(K, *mptr[i], ii);
1705 }
1706 if (err) {
1707 fprintf(stderr, "ksmooth_refresh_matrices: err = %d at t = %d\n",
1708 err, K->t);
1709 }
1710 }
1711 }
1712
1713 if (!err && cross_update) {
1714 /* cross-correlated case */
1715 err = kalman_update_crossinfo(K, UPDATE_STEP);
1716 }
1717
1718 return err;
1719 }
1720
1721 /**
1722 * kalman_forecast:
1723 * @K: pointer to Kalman struct: see kalman_new().
1724 * @prn: printing apparatus (or NULL).
1725 *
1726 * Generates a series of one-step ahead forecasts for y, based on
1727 * information entered initially using kalman_new(), and possibly
1728 * modified using kalman_set_initial_state_vector() and/or
1729 * kalman_set_initial_MSE_matrix(). The log-likelihood is
1730 * calculated for the sequence of forecast errors on the assumption
1731 * of normality: this can be accessed using kalman_get_loglik().
1732 *
1733 * Returns: 0 on success, non-zero on error.
1734 */
1735
kalman_forecast(kalman * K,PRN * prn)1736 int kalman_forecast (kalman *K, PRN *prn)
1737 {
1738 double ldet;
1739 int smoothing, update_P = 1;
1740 int Tmiss = 0;
1741 int i, err = 0;
1742
1743 #if KDEBUG
1744 fprintf(stderr, "kalman_forecast: T = %d\n", K->T);
1745 #endif
1746
1747 smoothing = (K->flags & KALMAN_SMOOTH)? 1 : 0;
1748
1749 if (K->nonshift < 0) {
1750 K->nonshift = count_nonshifts(K->F);
1751 }
1752
1753 K->SSRw = K->sumldet = K->loglik = 0.0;
1754 K->s2 = NADBL;
1755 K->okT = K->T;
1756
1757 if (K->x == NULL) {
1758 /* no exogenous vars */
1759 if (K->A != NULL) {
1760 /* implicit const case: A is 1 x n and A'x is n x 1 */
1761 gretl_vector_copy_values(K->Ax, K->A);
1762 } else {
1763 gretl_matrix_zero(K->Ax);
1764 }
1765 }
1766
1767 set_kalman_running(K);
1768
1769 for (K->t = 0; K->t < K->T && !err; K->t += 1) {
1770 int missobs = 0;
1771 double llt = 0.0;
1772
1773 #if KDEBUG > 1
1774 kalman_print_state(K);
1775 #endif
1776
1777 if (K->S != NULL || K->P != NULL) {
1778 kalman_record_state(K);
1779 }
1780
1781 if (filter_is_varying(K)) {
1782 /* we have time-varying coefficients */
1783 err = kalman_refresh_matrices(K, prn);
1784 if (err) {
1785 K->loglik = NADBL;
1786 break;
1787 }
1788 }
1789
1790 /* read slice from y */
1791 kalman_initialize_error(K, &missobs);
1792
1793 if (K->x != NULL) {
1794 /* read from x if applicable
1795 note 2010-09-18: this was conditional on missobs < K->n
1796 FIXME?
1797 */
1798 kalman_set_Ax(K, &missobs);
1799 }
1800
1801 if (missobs) {
1802 /* 2010-09-18: this was right after kalman_initialize_error
1803 FIXME?
1804 */
1805 Tmiss++;
1806 }
1807
1808 /* initial matrix calculations: form PH and H'PH
1809 (note that we need PH later) */
1810 gretl_matrix_multiply(K->P0, K->H, K->PH);
1811 if (K->n == 1) {
1812 /* slight speed-up for univariate observable */
1813 double x = (K->R == NULL)? 0.0 : K->R->val[0];
1814
1815 for (i=0; i<K->r; i++) {
1816 x += K->H->val[i] * K->PH->val[i];
1817 }
1818 if (x <= 0.0) {
1819 err = E_NAN;
1820 } else {
1821 K->HPH->val[0] = x;
1822 ldet = log(x);
1823 K->Vt->val[0] = 1.0 / x;
1824 }
1825 } else {
1826 gretl_matrix_qform(K->H, GRETL_MOD_TRANSPOSE,
1827 K->P0, K->HPH, GRETL_MOD_NONE);
1828 if (K->R != NULL) {
1829 gretl_matrix_add_to(K->HPH, K->R);
1830 }
1831 gretl_matrix_copy_values(K->Vt, K->HPH);
1832 err = gretl_invert_symmetric_matrix2(K->Vt, &ldet);
1833 if (err) {
1834 fprintf(stderr, "kalman_forecast: failed to invert V\n");
1835 gretl_matrix_print(K->Vt, "V");
1836 }
1837 }
1838
1839 /* likelihood bookkeeping */
1840 if (err) {
1841 K->loglik = NADBL;
1842 break;
1843 } else if (!missobs) {
1844 K->sumldet += ldet;
1845 }
1846
1847 if (K->V != NULL) {
1848 /* record MSE for estimate of observables */
1849 if (missobs) {
1850 if (smoothing) {
1851 set_row(K->V, K->t, 0.0);
1852 } else {
1853 set_row(K->V, K->t, 1.0/0.0);
1854 }
1855 } else if (smoothing) {
1856 /* record inverse */
1857 load_to_vech(K->V, K->Vt, K->n, K->t);
1858 } else {
1859 /* record uninverted matrix */
1860 load_to_vech(K->V, K->HPH, K->n, K->t);
1861 }
1862 }
1863
1864 /* first stage of dual iteration */
1865 if (arma_ll(K) && !smoothing) {
1866 err = kalman_arma_iter_1(K, missobs);
1867 } else {
1868 err = kalman_iter_1(K, missobs, &llt);
1869 if (K->LL != NULL) {
1870 if (na(llt) || missobs) {
1871 llt = NADBL;
1872 } else {
1873 llt -= 0.5 * (K->n * LN_2_PI + ldet);
1874 }
1875 gretl_vector_set(K->LL, K->t, llt);
1876 }
1877 }
1878
1879 /* record forecast errors if wanted */
1880 if (!err && K->E != NULL) {
1881 if (missobs && smoothing) {
1882 set_row(K->E, K->t, 0.0);
1883 } else {
1884 load_to_row(K->E, K->e, K->t);
1885 }
1886 }
1887
1888 /* update state vector */
1889 if (!err) {
1890 gretl_matrix_copy_values(K->S0, K->S1);
1891 }
1892
1893 if (!err && update_P) {
1894 /* second stage of dual iteration */
1895 err = kalman_iter_2(K, missobs);
1896 }
1897
1898 if (!err) {
1899 /* update MSE matrix, if needed */
1900 if (arma_ll(K) && !smoothing && update_P && K->t > 20) {
1901 if (!matrix_diff(K->P1, K->P0, 1.0e-20)) {
1902 K->P0->val[0] += 1.0;
1903 update_P = 0;
1904 }
1905 }
1906 if (update_P) {
1907 gretl_matrix_copy_values(K->P0, K->P1);
1908 }
1909 }
1910 }
1911
1912 set_kalman_stopped(K);
1913
1914 if (isnan(K->loglik) || isinf(K->loglik)) {
1915 K->loglik = NADBL;
1916 }
1917
1918 if (na(K->loglik)) {
1919 goto bailout;
1920 }
1921
1922 K->okT = K->T - Tmiss;
1923
1924 if (arma_ll(K)) {
1925 double ll1 = 1.0 + LN_2_PI + log(K->SSRw / K->okT);
1926
1927 if (K->flags & KALMAN_AVG_LL) {
1928 K->loglik = -0.5 * (ll1 + K->sumldet / K->okT);
1929 } else {
1930 K->loglik = -0.5 * (K->okT * ll1 + K->sumldet);
1931 }
1932 } else {
1933 /* For K->s2 see Koopman, Shephard and Doornik, "Statistical
1934 algorithms for models in state space using SsfPack 2.2",
1935 Econometrics Journal, 1999, vol. 2, pp. 113-166; also
1936 available at http://www.ssfpack.com/ . For the role of
1937 'd', see in addition Koopman's 1997 JASA article.
1938 */
1939 int nT = K->n * K->okT;
1940 int d = (K->flags & KALMAN_DIFFUSE)? K->r : 0;
1941
1942 if (d > 0) {
1943 K->loglik = -0.5 * ((nT - d) * LN_2_PI + K->sumldet + K->SSRw)
1944 + (d / 2.0) * log(kappa);
1945 } else {
1946 K->loglik = -0.5 * (nT * LN_2_PI + K->sumldet + K->SSRw);
1947 }
1948 K->s2 = K->SSRw / (nT - d);
1949 }
1950
1951 if (isnan(K->loglik) || isinf(K->loglik)) {
1952 K->loglik = NADBL;
1953 err = E_NAN;
1954 }
1955
1956 bailout:
1957
1958 #if KDEBUG
1959 fprintf(stderr, "kalman_forecast: err=%d, ll=%#.12g\n", err,
1960 K->loglik);
1961 #endif
1962
1963 return err;
1964 }
1965
1966 /**
1967 * kalman_get_loglik:
1968 * @K: pointer to Kalman struct.
1969 *
1970 * Retrieves the log-likelhood calculated via a run of
1971 * kalman_forecast().
1972 *
1973 * Returns: ll value, or #NADBL on failure.
1974 */
1975
kalman_get_loglik(const kalman * K)1976 double kalman_get_loglik (const kalman *K)
1977 {
1978 return K->loglik;
1979 }
1980
1981 /**
1982 * kalman_get_arma_variance:
1983 * @K: pointer to Kalman struct.
1984 *
1985 * Retrieves the estimated variance for an ARMA model
1986 * estimated using the Kalman filter.
1987 *
1988 * Returns: sigma-squared value, or #NADBL on failure.
1989 */
1990
kalman_get_arma_variance(const kalman * K)1991 double kalman_get_arma_variance (const kalman *K)
1992 {
1993 if (na(K->SSRw)) {
1994 return NADBL;
1995 } else {
1996 return K->SSRw / K->okT;
1997 }
1998 }
1999
2000 /**
2001 * kalman_set_initial_state_vector:
2002 * @K: pointer to Kalman struct.
2003 * @S: matrix of values to set.
2004 *
2005 * Resets the initial value of the state vector in a Kalman
2006 * struct, using the values from @S. See also kalman_new().
2007 *
2008 * Returns: 0 on success, non-zero on error.
2009 */
2010
kalman_set_initial_state_vector(kalman * K,const gretl_matrix * S)2011 int kalman_set_initial_state_vector (kalman *K, const gretl_matrix *S)
2012 {
2013 return gretl_matrix_copy_values(K->S0, S);
2014 }
2015
2016 /**
2017 * kalman_set_initial_MSE_matrix:
2018 * @K: pointer to Kalman struct.
2019 * @P: matrix of values to set.
2020 *
2021 * Resets the initial value of the MSE matrix in a Kalman
2022 * struct, using the values from @P. See also kalman_new().
2023 *
2024 * Returns: 0 on success, non-zero on error.
2025 */
2026
kalman_set_initial_MSE_matrix(kalman * K,const gretl_matrix * P)2027 int kalman_set_initial_MSE_matrix (kalman *K, const gretl_matrix *P)
2028 {
2029 return gretl_matrix_copy_values(K->P0, P);
2030 }
2031
2032 struct K_input_mat {
2033 int sym;
2034 const char *name;
2035 };
2036
2037 /* mapping to names used in "kalman ... end kalman" block */
2038
2039 struct K_input_mat K_input_mats[] = {
2040 { K_y, "obsy" },
2041 { K_H, "obsymat" },
2042 { K_x, "obsx" },
2043 { K_A, "obsxmat" },
2044 { K_R, "obsvar" },
2045 { K_F, "statemat" },
2046 { K_Q, "statevar" },
2047 { K_m, "stconst" },
2048 { K_S, "inistate" },
2049 { K_P, "inivar" }
2050 };
2051
obsy_check(kalman * K)2052 static int obsy_check (kalman *K)
2053 {
2054 if (K->y == NULL) {
2055 return missing_matrix_error("obsy");
2056 } else if (K->y->rows != K->T || K->y->cols != K->n) {
2057 fprintf(stderr, "obsy_check: K->y should be %d x %d, is %d x %d\n",
2058 K->T, K->n, gretl_matrix_rows(K->y),
2059 gretl_matrix_cols(K->y));
2060 return E_NONCONF;
2061 } else {
2062 return 0;
2063 }
2064 }
2065
kalman_bundle_recheck_matrices(kalman * K,PRN * prn)2066 static int kalman_bundle_recheck_matrices (kalman *K, PRN *prn)
2067 {
2068 int err = 0;
2069
2070 K->flags |= KALMAN_CHECK;
2071
2072 if (filter_is_varying(K)) {
2073 err = kalman_update_matrices(K, prn);
2074 }
2075
2076 K->flags ^= KALMAN_CHECK;
2077
2078 if (!err && (K->H == NULL || K->F == NULL || K->Q == NULL)) {
2079 fprintf(stderr, "kalman_bundle_kalman_recheck_matrices: H=%p, F=%p, Q=%p\n",
2080 K->H, K->F, K->Q);
2081 err = missing_matrix_error(NULL);
2082 }
2083
2084 if (err) {
2085 return err;
2086 }
2087
2088 /* redundant? */
2089 kalman_set_dimensions(K);
2090
2091 if (gretl_matrix_rows(K->F) != K->r ||
2092 gretl_matrix_rows(K->A) != K->k) {
2093 err = E_NONCONF;
2094 } else if (!kalman_simulating(K)) {
2095 err = obsy_check(K);
2096 }
2097
2098 if (!err && K->p > 0) {
2099 err = kalman_revise_variance(K);
2100 }
2101
2102 if (!err) {
2103 err = kalman_check_dimensions(K);
2104 }
2105
2106 if (!err) {
2107 if (K->Sini != NULL) {
2108 gretl_matrix_copy_values(K->S0, K->Sini);
2109 } else {
2110 gretl_matrix_zero(K->S0);
2111 }
2112 if (K->Pini != NULL) {
2113 gretl_matrix_copy_values(K->P0, K->Pini);
2114 } else {
2115 err = construct_Pini(K);
2116 }
2117 }
2118
2119 return err;
2120 }
2121
kalman_attach_data(kalman * K,void * data)2122 void kalman_attach_data (kalman *K, void *data)
2123 {
2124 if (K != NULL) {
2125 K->data = data;
2126 }
2127 }
2128
kalman_get_data(const kalman * K)2129 void *kalman_get_data (const kalman *K)
2130 {
2131 return (K != NULL)? K->data : NULL;
2132 }
2133
kalman_attach_printer(kalman * K,PRN * prn)2134 void kalman_attach_printer (kalman *K, PRN *prn)
2135 {
2136 if (K != NULL) {
2137 K->prn = prn;
2138 }
2139 }
2140
kalman_get_printer(const kalman * K)2141 PRN *kalman_get_printer (const kalman *K)
2142 {
2143 return (K != NULL)? K->prn : NULL;
2144 }
2145
kalman_matrix_name(int sym)2146 static const char *kalman_matrix_name (int sym)
2147 {
2148 int i;
2149
2150 for (i=0; i<K_MMAX; i++) {
2151 if (K_input_mats[i].sym == sym) {
2152 return K_input_mats[i].name;
2153 }
2154 }
2155
2156 /* failed */
2157 return "matrix";
2158 }
2159
kalman_ensure_output_matrices(kalman * K)2160 static int kalman_ensure_output_matrices (kalman *K)
2161 {
2162 int err = 0;
2163
2164 if (K->E == NULL) {
2165 K->E = gretl_null_matrix_new();
2166 }
2167 if (K->V == NULL) {
2168 K->V = gretl_null_matrix_new();
2169 }
2170 if (K->S == NULL) {
2171 K->S = gretl_null_matrix_new();
2172 }
2173 if (K->P == NULL) {
2174 K->P = gretl_null_matrix_new();
2175 }
2176 if (K->K == NULL) {
2177 K->K = gretl_null_matrix_new();
2178 }
2179
2180 if (K->E == NULL || K->V == NULL || K->S == NULL ||
2181 K->P == NULL || K->K == NULL) {
2182 err = E_ALLOC;
2183 }
2184
2185 return err;
2186 }
2187
kalman_bundle_run(gretl_bundle * b,PRN * prn,int * errp)2188 int kalman_bundle_run (gretl_bundle *b, PRN *prn, int *errp)
2189 {
2190 kalman *K = gretl_bundle_get_private_data(b);
2191 int err;
2192
2193 K->b = b; /* attach bundle pointer */
2194 err = kalman_ensure_output_matrices(K);
2195
2196 if (!err) {
2197 gretl_matrix_zero(K->e);
2198 err = kalman_bundle_recheck_matrices(K, prn);
2199 }
2200
2201 if (!err && K->LL == NULL) {
2202 K->LL = gretl_matrix_alloc(K->T, 1);
2203 if (K->LL == NULL) {
2204 err = E_ALLOC;
2205 }
2206 }
2207
2208 if (!err) {
2209 err = kalman_forecast(K, prn);
2210 }
2211
2212 if (err != E_NAN) {
2213 *errp = err;
2214 } else {
2215 /* we'll flag E_NAN with a return value of 1 but
2216 won't count it as a 'true' error */
2217 *errp = 0;
2218 }
2219
2220 return err;
2221 }
2222
2223 /* Copy row @t from @src into @targ; or add row @t of @src to
2224 @targ; or subtract row @t of @src from @targ. We allow the
2225 possibility that the length of vector @targ is less than
2226 the number of columns in @src, but not the converse.
2227 */
2228
load_from_row(gretl_vector * targ,const gretl_matrix * src,int t,GretlMatrixMod mod)2229 static int load_from_row (gretl_vector *targ,
2230 const gretl_matrix *src,
2231 int t, GretlMatrixMod mod)
2232 {
2233 int i, n = gretl_vector_get_length(targ);
2234 double x;
2235
2236 if (n > src->cols) {
2237 fprintf(stderr, "load_from_row: targ length = %d, but src "
2238 "has %d columns\n", n, src->cols);
2239 return E_NONCONF;
2240 }
2241
2242 for (i=0; i<n; i++) {
2243 x = gretl_matrix_get(src, t, i);
2244 if (mod == GRETL_MOD_CUMULATE) {
2245 targ->val[i] += x;
2246 } else if (mod == GRETL_MOD_DECREMENT) {
2247 targ->val[i] -= x;
2248 } else {
2249 targ->val[i] = x;
2250 }
2251 }
2252
2253 return 0;
2254 }
2255
2256 /* As load_from_row(), except that a column offset, @j, is
2257 supported for the reading of a row from @src, and we
2258 don't support the @mod option.
2259 */
2260
load_from_row_offset(gretl_vector * targ,const gretl_matrix * src,int t,int j)2261 static int load_from_row_offset (gretl_vector *targ,
2262 const gretl_matrix *src,
2263 int t, int j)
2264 {
2265 int i, n = gretl_vector_get_length(targ);
2266 double x;
2267
2268 if (n > src->cols - j) {
2269 return E_NONCONF;
2270 }
2271
2272 for (i=0; i<n; i++) {
2273 x = gretl_matrix_get(src, t, i + j);
2274 targ->val[i] = x;
2275 }
2276
2277 return 0;
2278 }
2279
2280 /* Row @t of @src represents the vech of an n x n matrix: extract the
2281 row and apply the inverse operation of vech to reconstitute the
2282 matrix in @targ -- or subtract the newly reconstituted matrix
2283 from @targ.
2284 */
2285
load_from_vech(gretl_matrix * targ,const gretl_matrix * src,int n,int t,int mod)2286 static void load_from_vech (gretl_matrix *targ, const gretl_matrix *src,
2287 int n, int t, int mod)
2288 {
2289 int i, j, m = 0;
2290 double x;
2291
2292 for (i=0; i<n; i++) {
2293 for (j=i; j<n; j++) {
2294 x = gretl_matrix_get(src, t, m++);
2295 if (mod == GRETL_MOD_DECREMENT) {
2296 x = gretl_matrix_get(targ, i, j) - x;
2297 }
2298 gretl_matrix_set(targ, i, j, x);
2299 if (i != j) {
2300 gretl_matrix_set(targ, j, i, x);
2301 }
2302 }
2303 }
2304 }
2305
2306 /* Row @t of @src represents the vec of a certain matrix: extract the
2307 row and reconstitute the matrix in @targ.
2308 */
2309
load_from_vec(gretl_matrix * targ,const gretl_matrix * src,int t)2310 static int load_from_vec (gretl_matrix *targ, const gretl_matrix *src,
2311 int t)
2312 {
2313 int i, k = targ->rows * targ->cols;
2314
2315 for (i=0; i<k; i++) {
2316 targ->val[i] = gretl_matrix_get(src, t, i);
2317 }
2318
2319 return 0;
2320 }
2321
2322 /* For disturbance smoothing: ensure we have on hand
2323 matrices that are correctly sized to hold estimates
2324 of the variance of the disturbance(s) in the state
2325 and (if applicable) observation equations.
2326 */
2327
maybe_resize_dist_mse(kalman * K,gretl_matrix ** vvt,gretl_matrix ** vwt)2328 static int maybe_resize_dist_mse (kalman *K,
2329 gretl_matrix **vvt,
2330 gretl_matrix **vwt)
2331 {
2332 int n = K->R == NULL ? 0 : K->n;
2333 int k, err = 0;
2334
2335 /* combined results: how many columns do we need? */
2336 k = K->r + n;
2337
2338 if (K->Vsd == NULL) {
2339 K->Vsd = gretl_matrix_alloc(K->T, k);
2340 if (K->Vsd == NULL) {
2341 err = E_ALLOC;
2342 }
2343 } else if (K->Vsd->rows != K->T || K->Vsd->cols != k) {
2344 err = gretl_matrix_realloc(K->Vsd, K->T, k);
2345 }
2346
2347 if (!err) {
2348 /* step-t square state matrix */
2349 *vvt = gretl_matrix_alloc(K->r, K->r);
2350 if (*vvt == NULL) {
2351 err = E_ALLOC;
2352 }
2353 }
2354
2355 if (!err && n > 0) {
2356 /* step-t square obs matrix */
2357 *vwt = gretl_matrix_alloc(K->n, K->n);
2358 if (*vwt == NULL) {
2359 err = E_ALLOC;
2360 }
2361 }
2362
2363 return err;
2364 }
2365
retrieve_Ft(kalman * K,int t)2366 static int retrieve_Ft (kalman *K, int t)
2367 {
2368 if (K->step == NULL || K->step->F == NULL) {
2369 return E_DATA;
2370 } else {
2371 return load_from_vec((gretl_matrix *) K->F, K->step->F, t);
2372 }
2373 }
2374
retrieve_Ht(kalman * K,int t)2375 static int retrieve_Ht (kalman *K, int t)
2376 {
2377 if (K->step == NULL || K->step->H == NULL) {
2378 return E_DATA;
2379 } else {
2380 return load_from_vec((gretl_matrix *) K->H, K->step->H, t);
2381 }
2382 }
2383
2384 /* Calculate the variance of the smoothed disturbances
2385 for the cross-correlated case. See Koopman, Shephard
2386 and Doornik (1998), page 19, var(\varepsilon_t|Y_n).
2387 */
2388
combined_dist_variance(kalman * K,gretl_matrix * D,gretl_matrix * N,gretl_matrix * vv,gretl_matrix * vw,gretl_matrix_block * BX,int dkstyle)2389 static int combined_dist_variance (kalman *K,
2390 gretl_matrix *D,
2391 gretl_matrix *N,
2392 gretl_matrix *vv,
2393 gretl_matrix *vw,
2394 gretl_matrix_block *BX,
2395 int dkstyle)
2396 {
2397 gretl_matrix *DC, *KN, *Veps, *NB, *NK;
2398
2399 DC = gretl_matrix_block_get_matrix(BX, 0);
2400 KN = gretl_matrix_block_get_matrix(BX, 1);
2401 Veps = gretl_matrix_block_get_matrix(BX, 2);
2402 NB = gretl_matrix_block_get_matrix(BX, 3);
2403
2404 KN = gretl_matrix_reuse(KN, K->n, K->r);
2405
2406 /* First chunk of Veps:
2407 Koopman's notation: G_t'(D_t*G_t - K_t'*N_t*H_t)
2408 In gretl: C_t'(D_t*C_t - K_t'*N_t*B_t)
2409 */
2410
2411 gretl_matrix_multiply(D, K->C, DC);
2412 gretl_matrix_multiply_mod(K->Kt, GRETL_MOD_TRANSPOSE,
2413 N, GRETL_MOD_NONE,
2414 KN, GRETL_MOD_NONE);
2415 gretl_matrix_multiply_mod(KN, GRETL_MOD_NONE,
2416 K->B, GRETL_MOD_NONE,
2417 DC, GRETL_MOD_DECREMENT);
2418 gretl_matrix_multiply_mod(K->C, GRETL_MOD_TRANSPOSE,
2419 DC, GRETL_MOD_NONE,
2420 Veps, GRETL_MOD_NONE);
2421
2422 NK = gretl_matrix_reuse(KN, K->r, K->n);
2423
2424 /* Second chunk of Veps:
2425 Koopman's notation: H_t'(N_t*H_t - N_t*K_t*G_t)
2426 In gretl: B_t'(N_t*B_t - N_t*K_t*C_t),
2427 and add to the first chunk calculated above.
2428 */
2429
2430 gretl_matrix_multiply(N, K->B, NB);
2431 gretl_matrix_multiply_mod(N, GRETL_MOD_TRANSPOSE,
2432 K->Kt, GRETL_MOD_NONE,
2433 NK, GRETL_MOD_NONE);
2434 gretl_matrix_multiply_mod(NK, GRETL_MOD_NONE,
2435 K->C, GRETL_MOD_NONE,
2436 NB, GRETL_MOD_DECREMENT);
2437 gretl_matrix_multiply_mod(K->B, GRETL_MOD_TRANSPOSE,
2438 NB, GRETL_MOD_NONE,
2439 Veps, GRETL_MOD_CUMULATE);
2440
2441 if (dkstyle) {
2442 /* Veps = I_p - Veps */
2443 double vii;
2444 int i;
2445
2446 gretl_matrix_multiply_by_scalar(Veps, -1.0);
2447
2448 for (i=0; i<K->p; i++) {
2449 vii = gretl_matrix_get(Veps, i, i);
2450 gretl_matrix_set(Veps, i, i, 1.0 + vii);
2451 }
2452 }
2453
2454 /* Veps (p x p) holds the variance of \epsilon_t
2455 conditional on Y_n: now form the per-equation
2456 disturbance variance matrices, @vv and @vw, for
2457 this time-step.
2458 */
2459
2460 gretl_matrix_qform(K->B, GRETL_MOD_NONE, Veps,
2461 vv, GRETL_MOD_NONE);
2462 gretl_matrix_qform(K->C, GRETL_MOD_NONE, Veps,
2463 vw, GRETL_MOD_NONE);
2464
2465 return 0;
2466 }
2467
2468 /* See Koopman, Shephard and Doornik, section 4.4 */
2469
koopman_smooth(kalman * K,int dkstyle)2470 static int koopman_smooth (kalman *K, int dkstyle)
2471 {
2472 gretl_matrix_block *B, *BX = NULL;
2473 gretl_matrix *u, *D, *L, *R;
2474 gretl_matrix *r0, *r1, *r2, *N1, *N2;
2475 gretl_matrix *n1 = NULL;
2476 gretl_matrix *Vvt = NULL;
2477 gretl_matrix *Vwt = NULL;
2478 gretl_matrix *DC = NULL;
2479 gretl_matrix *KN = NULL;
2480 gretl_matrix *RHS = NULL;
2481 gretl_matrix *NB = NULL;
2482 gretl_matrix *Ut = NULL;
2483 double x;
2484 int i, t, err = 0;
2485
2486 B = gretl_matrix_block_new(&u, K->n, 1,
2487 &D, K->n, K->n,
2488 &L, K->r, K->r,
2489 &R, K->T, K->r,
2490 &r0, K->r, 1,
2491 &r1, K->r, 1,
2492 &r2, K->r, 1,
2493 &N1, K->r, K->r,
2494 &N2, K->r, K->r,
2495 NULL);
2496
2497 if (B == NULL) {
2498 return E_ALLOC;
2499 }
2500
2501 if (K->b != NULL) {
2502 /* for variance of smoothed disturbances */
2503 err = maybe_resize_dist_mse(K, &Vvt, &Vwt);
2504 }
2505
2506 if (K->b != NULL && K->p > 0) {
2507 BX = gretl_matrix_block_new(&DC, K->n, K->p,
2508 &KN, K->n, K->r,
2509 &RHS, K->p, K->p,
2510 &NB, K->r, K->p,
2511 &Ut, K->p, 1,
2512 NULL);
2513 if (BX == NULL) {
2514 err = E_ALLOC;
2515 }
2516 }
2517
2518 if (err) {
2519 gretl_matrix_block_destroy(B);
2520 gretl_matrix_block_destroy(BX);
2521 gretl_matrix_free(Vvt);
2522 gretl_matrix_free(Vwt);
2523 return err;
2524 }
2525
2526 gretl_matrix_zero(r1);
2527 gretl_matrix_zero(N1);
2528
2529 /* The backward recursion */
2530
2531 for (t=K->T-1; t>=0 && !err; t--) {
2532 /* load et, Vt and Kt for time t */
2533 load_from_row(K->e, K->E, t, GRETL_MOD_NONE);
2534 load_from_vech(K->Vt, K->V, K->n, t, GRETL_MOD_NONE);
2535 load_from_vec(K->Kt, K->K, t);
2536
2537 /* get F_t and/or H_t if need be */
2538 if (matrix_is_varying(K, K_F)) {
2539 err = retrieve_Ft(K, t);
2540 }
2541 if (!err && matrix_is_varying(K, K_H)) {
2542 err = retrieve_Ht(K, t);
2543 }
2544 if (err) {
2545 break;
2546 }
2547
2548 if (filter_is_varying(K)) {
2549 /* K->Q and/or K->R may be time-varying */
2550 K->t = t;
2551 ksmooth_refresh_matrices(K, NULL);
2552 }
2553
2554 /* u_t = V_t^{-1} e_t - K_t' r_t */
2555 gretl_matrix_multiply(K->Vt, K->e, u);
2556 if (t < K->T - 1) {
2557 gretl_matrix_multiply_mod(K->Kt, GRETL_MOD_TRANSPOSE,
2558 r1, GRETL_MOD_NONE,
2559 u, GRETL_MOD_DECREMENT);
2560 }
2561
2562 /* save u_t values in E */
2563 load_to_row(K->E, u, t);
2564
2565 if (K->Vsd != NULL && K->p == 0) {
2566 /* variance of state disturbance */
2567 if (dkstyle) {
2568 /* Q_t - Q_t N_t Q_t */
2569 gretl_matrix_copy_values(Vvt, K->Q);
2570 gretl_matrix_qform(K->Q, GRETL_MOD_TRANSPOSE,
2571 N1, Vvt, GRETL_MOD_DECREMENT);
2572
2573 } else {
2574 /* Q_t N_t Q_t */
2575 gretl_matrix_qform(K->Q, GRETL_MOD_TRANSPOSE,
2576 N1, Vvt, GRETL_MOD_NONE);
2577 }
2578 load_to_diag(K->Vsd, Vvt, t, 0);
2579 }
2580
2581 /* D_t = V_t^{-1} + K_t' N_t K_t */
2582 gretl_matrix_copy_values(D, K->Vt);
2583 if (t < K->T - 1) {
2584 gretl_matrix_qform(K->Kt, GRETL_MOD_TRANSPOSE,
2585 N1, D, GRETL_MOD_CUMULATE);
2586 }
2587
2588 if (K->R != NULL && K->Vsd != NULL && K->p == 0) {
2589 /* variance of obs disturbance */
2590 if (dkstyle) {
2591 /* R_t - R_t D_t R_t */
2592 gretl_matrix_copy_values(Vwt, K->R);
2593 gretl_matrix_qform(K->R, GRETL_MOD_TRANSPOSE,
2594 D, Vwt, GRETL_MOD_DECREMENT);
2595
2596 } else {
2597 /* R_t D_t R_t */
2598 gretl_matrix_qform(K->R, GRETL_MOD_TRANSPOSE,
2599 D, Vwt, GRETL_MOD_NONE);
2600 }
2601 load_to_diag(K->Vsd, Vwt, t, K->r);
2602 }
2603
2604 if (K->Vsd != NULL && K->p > 0) {
2605 /* variance of combined disturbance */
2606 err = combined_dist_variance(K, D, N1, Vvt, Vwt, BX,
2607 dkstyle);
2608 if (err) {
2609 break;
2610 } else {
2611 load_to_diag(K->Vsd, Vvt, t, 0);
2612 load_to_diag(K->Vsd, Vwt, t, K->r);
2613 }
2614 }
2615
2616 /* L_t = F - KH' */
2617 gretl_matrix_copy_values(L, K->F);
2618 gretl_matrix_multiply_mod(K->Kt, GRETL_MOD_NONE,
2619 K->H, GRETL_MOD_TRANSPOSE,
2620 L, GRETL_MOD_DECREMENT);
2621
2622 /* save r_t values in R */
2623 load_to_row(R, r1, t);
2624
2625 /* r_{t-1} = H V_t^{-1}e_t + L_t' r_t */
2626 gretl_matrix_multiply(K->H, K->Vt, K->Tmpr1);
2627 gretl_matrix_multiply(K->Tmpr1, K->e, r2);
2628 if (t < K->T - 1) {
2629 gretl_matrix_multiply_mod(L, GRETL_MOD_TRANSPOSE,
2630 r1, GRETL_MOD_NONE,
2631 r2, GRETL_MOD_CUMULATE);
2632 }
2633 /* transcribe for next step */
2634 gretl_matrix_copy_values(r1, r2);
2635
2636 /* preserve r_0 for smoothing of state */
2637 if (t == 0) {
2638 gretl_matrix_copy_values(r0, r2);
2639 }
2640
2641 /* N_{t-1} = H V_t^{-1}H' + L' N L */
2642 gretl_matrix_qform(K->H, GRETL_MOD_NONE,
2643 K->Vt, N2, GRETL_MOD_NONE);
2644 if (t < K->T - 1) {
2645 gretl_matrix_qform(L, GRETL_MOD_TRANSPOSE,
2646 N1, N2, GRETL_MOD_CUMULATE);
2647 }
2648 /* transcribe for next step */
2649 gretl_matrix_copy_values(N1, N2);
2650 }
2651
2652 #if 0
2653 gretl_matrix_print(K->E, "u_t, all t");
2654 gretl_matrix_print(R, "r_t, all t");
2655 #endif
2656
2657 if (K->R != NULL && K->U != NULL) {
2658 /* we need an n x 1 temp matrix */
2659 n1 = gretl_matrix_reuse(K->Tmpnn, K->n, 1);
2660 }
2661
2662 /* Smoothed disturbances, all time steps */
2663
2664 if (K->p > 0) {
2665 /* eps = B' r_t + C' e_t */
2666 for (t=0; t<K->T; t++) {
2667 if (filter_is_varying(K)) {
2668 K->t = t;
2669 ksmooth_refresh_matrices(K, NULL);
2670 }
2671 load_from_row(r1, R, t, GRETL_MOD_NONE);
2672 gretl_matrix_multiply_mod(K->B, GRETL_MOD_TRANSPOSE,
2673 r1, GRETL_MOD_NONE,
2674 Ut, GRETL_MOD_NONE);
2675 load_from_row(K->e, K->E, t, GRETL_MOD_NONE);
2676 gretl_matrix_multiply_mod(K->C, GRETL_MOD_TRANSPOSE,
2677 K->e, GRETL_MOD_NONE,
2678 Ut, GRETL_MOD_CUMULATE);
2679 gretl_matrix_multiply(K->B, Ut, r2);
2680 load_to_row(R, r2, t);
2681 gretl_matrix_multiply(K->C, Ut, n1);
2682 for (i=0; i<K->r; i++) {
2683 x = gretl_vector_get(r2, i);
2684 gretl_matrix_set(K->U, t, i, x);
2685 }
2686 for (i=0; i<K->n; i++) {
2687 x = gretl_vector_get(n1, i);
2688 gretl_matrix_set(K->U, t, K->r + i, x);
2689 }
2690 }
2691 } else {
2692 /* the independent case */
2693 for (t=0; t<K->T; t++) {
2694 if (filter_is_varying(K)) {
2695 K->t = t;
2696 ksmooth_refresh_matrices(K, NULL);
2697 }
2698 load_from_row(r1, R, t, GRETL_MOD_NONE);
2699 gretl_matrix_multiply(K->Q, r1, r2);
2700 load_to_row(R, r2, t);
2701 if (K->U != NULL) {
2702 for (i=0; i<K->r; i++) {
2703 x = gretl_vector_get(r2, i);
2704 gretl_matrix_set(K->U, t, i, x);
2705 }
2706 }
2707 if (n1 != NULL) {
2708 load_from_row(K->e, K->E, t, GRETL_MOD_NONE);
2709 gretl_matrix_multiply(K->R, K->e, n1);
2710 for (i=0; i<K->n; i++) {
2711 x = gretl_vector_get(n1, i);
2712 gretl_matrix_set(K->U, t, K->r + i, x);
2713 }
2714 }
2715 }
2716 }
2717
2718 if (K->R != NULL && K->U != NULL) {
2719 /* reset to full size */
2720 gretl_matrix_reuse(K->Tmpnn, K->n, K->n);
2721 }
2722
2723 /* Write initial smoothed state into first row of S */
2724
2725 if (K->Pini != NULL) {
2726 gretl_matrix_multiply(K->Pini, r0, K->S0);
2727 } else {
2728 construct_Pini(K);
2729 gretl_matrix_multiply(K->P0, r0, K->S0);
2730 }
2731 if (K->Sini != NULL) {
2732 gretl_matrix_add_to(K->S0, K->Sini);
2733 }
2734 load_to_row(K->S, K->S0, 0);
2735
2736 /* Smoothed state, remaining time steps */
2737
2738 for (t=1; t<K->T; t++) {
2739 /* S_{t+1} = FS_t + v_t (or + B*eps_t) */
2740 load_from_row(K->S0, K->S, t-1, GRETL_MOD_NONE);
2741 gretl_matrix_multiply(K->F, K->S0, K->S1);
2742 load_from_row(K->S1, R, t-1, GRETL_MOD_CUMULATE);
2743 load_to_row(K->S, K->S1, t);
2744 }
2745
2746 gretl_matrix_block_destroy(B);
2747 gretl_matrix_block_destroy(BX);
2748 gretl_matrix_free(Vvt);
2749 gretl_matrix_free(Vwt);
2750
2751 return err;
2752 }
2753
2754 /* Anderson-Moore Kalman smoothing: see Iskander Karibzhanov's
2755 exposition at http://karibzhanov.com/help/kalcvs.htm
2756 This is much the clearest account I have seen (AC 2009-04-14,
2757 URL updated 2016-03-24).
2758
2759 This method uses S_{t|t-1} and P_{t|t-1} for all t, but we can
2760 overwrite these with the smoothed values as we go. We also need
2761 stored values for the prediction error, its MSE, and the gain at
2762 each time step. Note that u_t and U_t are set to zero for
2763 t = T - 1.
2764 */
2765
anderson_moore_smooth(kalman * K)2766 static int anderson_moore_smooth (kalman *K)
2767 {
2768 gretl_matrix_block *B;
2769 gretl_matrix *L = K->Tmprr;
2770 gretl_matrix *u, *u1, *U, *U1;
2771 gretl_matrix *StT, *PtT;
2772 int t, err = 0;
2773
2774 B = gretl_matrix_block_new(&StT, K->r, 1,
2775 &PtT, K->r, K->r,
2776 &u, K->r, 1,
2777 &u1, K->r, 1,
2778 &U, K->r, K->r,
2779 &U1, K->r, K->r,
2780 NULL);
2781 if (B == NULL) {
2782 return E_ALLOC;
2783 }
2784
2785 gretl_matrix_zero(u);
2786 gretl_matrix_zero(U);
2787
2788 for (t=K->T-1; t>=0 && !err; t--) {
2789 /* get F_t and/or H_t if need be */
2790 if (matrix_is_varying(K, K_F)) {
2791 err = retrieve_Ft(K, t);
2792 }
2793 if (!err && matrix_is_varying(K, K_H)) {
2794 err = retrieve_Ht(K, t);
2795 }
2796 if (err) {
2797 break;
2798 }
2799
2800 /* L_t = F_t - K_t H_t' */
2801 gretl_matrix_copy_values(L, K->F);
2802 load_from_vec(K->Kt, K->K, t);
2803 gretl_matrix_multiply_mod(K->Kt, GRETL_MOD_NONE,
2804 K->H, GRETL_MOD_TRANSPOSE,
2805 L, GRETL_MOD_DECREMENT);
2806
2807 /* u_{t-1} = H_t V_t e_t + L_t' u_t */
2808 load_from_vech(K->Vt, K->V, K->n, t, GRETL_MOD_NONE);
2809 load_from_row(K->e, K->E, t, GRETL_MOD_NONE);
2810 gretl_matrix_multiply(K->Vt, K->e, K->Ve);
2811 gretl_matrix_multiply(K->H, K->Ve, u1);
2812 if (t == K->T - 1) {
2813 gretl_matrix_multiply(K->H, K->Ve, u);
2814 } else {
2815 gretl_matrix_multiply_mod(L, GRETL_MOD_TRANSPOSE,
2816 u, GRETL_MOD_NONE,
2817 u1, GRETL_MOD_CUMULATE);
2818 gretl_matrix_copy_values(u, u1);
2819 }
2820
2821 /* U_{t-1} = H_t V_t H_t' + L_t' U_t L_t */
2822 if (t == K->T - 1) {
2823 gretl_matrix_qform(K->H, GRETL_MOD_NONE,
2824 K->Vt, U, GRETL_MOD_NONE);
2825 } else {
2826 gretl_matrix_qform(K->H, GRETL_MOD_NONE,
2827 K->Vt, U1, GRETL_MOD_NONE);
2828 gretl_matrix_qform(L, GRETL_MOD_TRANSPOSE,
2829 U, U1, GRETL_MOD_CUMULATE);
2830 gretl_matrix_copy_values(U, U1);
2831 }
2832
2833 /* S_{t|T} = S_{t|t-1} + P_{t|t-1} u_{t-1} */
2834 load_from_row(StT, K->S, t, GRETL_MOD_NONE);
2835 load_from_vech(K->P0, K->P, K->r, t, GRETL_MOD_NONE);
2836 gretl_matrix_multiply_mod(K->P0, GRETL_MOD_NONE,
2837 u, GRETL_MOD_NONE,
2838 StT, GRETL_MOD_CUMULATE);
2839 load_to_row(K->S, StT, t);
2840
2841 /* P_{t|T} = P_{t|t-1} - P_{t|t-1} U_{t-1} P_{t|t-1} */
2842 gretl_matrix_copy_values(PtT, K->P0);
2843 gretl_matrix_qform(K->P0, GRETL_MOD_NONE,
2844 U, PtT, GRETL_MOD_DECREMENT);
2845 load_to_vech(K->P, PtT, K->r, t);
2846 }
2847
2848 gretl_matrix_block_destroy(B);
2849
2850 return err;
2851 }
2852
2853 /* If we're doing smoothing for a system that has time-varying
2854 coefficients in K->F or K->H we'll record the vec of the
2855 coefficient matrices for each time-step on the forward pass.
2856 Here we allocate the required storage.
2857 */
2858
kalman_add_stepinfo(kalman * K)2859 static int kalman_add_stepinfo (kalman *K)
2860 {
2861 int err = 0;
2862
2863 K->step = malloc(sizeof *K->step);
2864
2865 if (K->step == NULL) {
2866 return E_ALLOC;
2867 }
2868
2869 K->step->F = K->step->H = NULL;
2870
2871 if (matrix_is_varying(K, K_F)) {
2872 K->step->F = gretl_matrix_alloc(K->T, K->r * K->r);
2873 if (K->step->F == NULL) {
2874 err = E_ALLOC;
2875 }
2876 }
2877
2878 if (!err && matrix_is_varying(K, K_H)) {
2879 K->step->H = gretl_matrix_alloc(K->T, K->r * K->n);
2880 if (K->step->H == NULL) {
2881 err = E_ALLOC;
2882 }
2883 }
2884
2885 if (err) {
2886 free_stepinfo(K);
2887 }
2888
2889 return err;
2890 }
2891
2892 /**
2893 * kalman_arma_smooth:
2894 * @K: pointer to Kalman struct.
2895 * @err: location to receive error code.
2896 *
2897 * Runs a filtering pass followed by a smoothing pass.
2898 *
2899 * Returns: matrix containing the smoothed estimate of the
2900 * dependent variable, or NULL on error.
2901 */
2902
kalman_arma_smooth(kalman * K,int * err)2903 gretl_matrix *kalman_arma_smooth (kalman *K, int *err)
2904 {
2905 int nr = (K->r * K->r + K->r) / 2;
2906 int nn = (K->n * K->n + K->n) / 2;
2907 gretl_vector *ys = NULL;
2908
2909 /* Set up the matrices we need to store computed results from all
2910 time steps on the forward pass: state S_{t|t-1},(inverse) error
2911 variance, gain, and MSE of state, P_{t|t-1}.
2912 */
2913
2914 K->S = gretl_matrix_alloc(K->T, K->r);
2915 K->V = gretl_matrix_alloc(K->T, nn);
2916 K->K = gretl_matrix_alloc(K->T, K->r * K->n);
2917 K->P = gretl_matrix_alloc(K->T, nr);
2918
2919 if (K->S == NULL || K->V == NULL || K->K == NULL || K->P == NULL) {
2920 *err = E_ALLOC;
2921 } else {
2922 double yst, Sti;
2923 int i, t, miss = 0;
2924
2925 K->flags |= KALMAN_SMOOTH;
2926 *err = kalman_forecast(K, NULL);
2927 K->flags &= ~KALMAN_SMOOTH;
2928 K->t = 0;
2929
2930 if (!*err) {
2931 *err = anderson_moore_smooth(K);
2932 }
2933
2934 if (!*err) {
2935 ys = gretl_column_vector_alloc(K->T);
2936 if (ys == NULL) {
2937 *err = E_ALLOC;
2938 } else {
2939 for (t=0; t<K->T; t++) {
2940 yst = 0.0;
2941 for (i=0; i<K->r; i++) {
2942 Sti = gretl_matrix_get(K->S, t, i);
2943 yst += K->H->val[i] * Sti;
2944 }
2945 if (K->Ax != NULL) {
2946 K->t = t;
2947 kalman_set_Ax(K, &miss);
2948 for (i=0; i<K->n; i++) {
2949 yst += gretl_vector_get(K->Ax, i);
2950 }
2951 }
2952 gretl_vector_set(ys, t, yst);
2953 }
2954 K->t = 0;
2955 }
2956 }
2957
2958 gretl_matrix_replace(&K->S, NULL);
2959 gretl_matrix_replace(&K->V, NULL);
2960 gretl_matrix_replace(&K->K, NULL);
2961 gretl_matrix_replace(&K->P, NULL);
2962 }
2963
2964 if (*err && ys != NULL) {
2965 gretl_matrix_free(ys);
2966 ys = NULL;
2967 }
2968
2969 return ys;
2970 }
2971
2972 /* optional matrix to hold smoothed disturbances a la Koopman */
2973
ensure_U_matrix(kalman * K)2974 static int ensure_U_matrix (kalman *K)
2975 {
2976 int Ucols = K->r;
2977 int Urows = K->T;
2978 int err = 0;
2979
2980 if (K->R != NULL) {
2981 Ucols += K->n;
2982 }
2983
2984 if (K->U == NULL) {
2985 K->U = gretl_matrix_alloc(Urows, Ucols);
2986 if (K->U == NULL) {
2987 err = E_ALLOC;
2988 }
2989 } else if (K->U->rows != Urows || K->U->cols != Ucols) {
2990 err = gretl_matrix_realloc(K->U, Urows, Ucols);
2991 }
2992
2993 return err;
2994 }
2995
2996 /**
2997 * kalman_smooth:
2998 * @K: pointer to kalman struct.
2999 * @pP: pointer to matrix in which to retrieve the MSE of the
3000 * smoothed state (or NULL if this is not required).
3001 * @pU: pointer to matrix in which to retrieve the smoothed
3002 * disturbances (or NULL if this is not required).
3003 * @err: location to receive error code.
3004 *
3005 * Runs a filtering pass followed by a backward, smoothing pass.
3006 * At present the @pU argument is experimental and a bodge: it will
3007 * not actually do anything unless @pP is left NULL.
3008 *
3009 * Returns: matrix containing the smoothed estimate of the
3010 * state, or NULL on error.
3011 */
3012
kalman_smooth(kalman * K,gretl_matrix ** pP,gretl_matrix ** pU,int * err)3013 gretl_matrix *kalman_smooth (kalman *K,
3014 gretl_matrix **pP,
3015 gretl_matrix **pU,
3016 int *err)
3017 {
3018 gretl_matrix *E, *S, *P = NULL;
3019 gretl_matrix *G, *V;
3020 int nr, nn;
3021
3022 if (pP == NULL && pU != NULL) {
3023 /* optional accessor for smoothed disturbances a la Koopman:
3024 experimental, and avilable only if @pP is not given
3025 */
3026 *err = ensure_U_matrix(K);
3027 if (*err) {
3028 return NULL;
3029 }
3030 }
3031
3032 nr = (K->r * K->r + K->r) / 2;
3033 nn = (K->n * K->n + K->n) / 2;
3034
3035 /* Set up the matrices we need to store computed results from all
3036 time steps on the forward pass: prediction error, (inverse)
3037 error variance, gain, state S_{t|t-1} and MSE of state,
3038 P_{t|t-1}.
3039 */
3040
3041 E = gretl_matrix_alloc(K->T, K->n);
3042 V = gretl_matrix_alloc(K->T, nn);
3043 G = gretl_matrix_alloc(K->T, K->r * K->n);
3044 S = gretl_matrix_alloc(K->T, K->r);
3045 P = gretl_matrix_alloc(K->T, nr);
3046
3047 if (E == NULL || V == NULL || G == NULL ||
3048 S == NULL || P == NULL) {
3049 *err = E_ALLOC;
3050 goto bailout;
3051 }
3052
3053 if (matrix_is_varying(K, K_F) || matrix_is_varying(K, K_H)) {
3054 /* add recorder for F_t and/or H_t */
3055 *err = kalman_add_stepinfo(K);
3056 if (*err) {
3057 goto bailout;
3058 }
3059 }
3060
3061 /* attach all export matrices to Kalman */
3062 K->E = E;
3063 K->V = V;
3064 K->K = G;
3065 K->S = S;
3066 K->P = P;
3067
3068 #if 0
3069 /* and recheck dimensions */
3070 *err = user_kalman_recheck_matrices(K, prn);
3071 #endif
3072
3073 if (!*err) {
3074 /* forward pass */
3075 K->flags |= KALMAN_SMOOTH;
3076 *err = kalman_forecast(K, NULL);
3077 K->flags &= ~KALMAN_SMOOTH;
3078 }
3079
3080 K->t = 0;
3081
3082 if (!*err) {
3083 /* bodge */
3084 if (K->U != NULL) {
3085 *err = koopman_smooth(K, 0);
3086 } else {
3087 *err = anderson_moore_smooth(K);
3088 }
3089 }
3090
3091 /* detach matrices */
3092 K->E = NULL;
3093 K->V = NULL;
3094 K->K = NULL;
3095 K->S = NULL;
3096 K->P = NULL;
3097
3098 /* and trash the "stepinfo" storage */
3099 free_stepinfo(K);
3100
3101 bailout:
3102
3103 gretl_matrix_free(E);
3104 gretl_matrix_free(V);
3105 gretl_matrix_free(G);
3106
3107 if (!*err && pP != NULL) {
3108 *pP = P;
3109 } else {
3110 gretl_matrix_free(P);
3111 }
3112
3113 if (!*err && pU != NULL) {
3114 *pU = K->U;
3115 } else {
3116 gretl_matrix_free(K->U);
3117 }
3118
3119 K->U = NULL;
3120
3121 if (*err) {
3122 gretl_matrix_free(S);
3123 S = NULL;
3124 }
3125
3126 return S;
3127 }
3128
kalman_bundle_smooth(gretl_bundle * b,int dist,PRN * prn)3129 int kalman_bundle_smooth (gretl_bundle *b, int dist, PRN *prn)
3130 {
3131 kalman *K = gretl_bundle_get_private_data(b);
3132 int err;
3133
3134 if (K == NULL) {
3135 fprintf(stderr, "kalman_bundle_smooth: K is NULL\n");
3136 return E_DATA;
3137 }
3138
3139 K->b = b; /* attach bundle pointer */
3140
3141 err = kalman_ensure_output_matrices(K);
3142
3143 if (!err && dist) {
3144 err = ensure_U_matrix(K);
3145 }
3146
3147 if (err) {
3148 return err;
3149 }
3150
3151 if (matrix_is_varying(K, K_F) || matrix_is_varying(K, K_H)) {
3152 /* add recorder for F_t and/or H_t */
3153 err = kalman_add_stepinfo(K);
3154 if (err) {
3155 goto bailout;
3156 }
3157 }
3158
3159 if (!err) {
3160 err = kalman_bundle_recheck_matrices(K, prn);
3161 }
3162
3163 if (!err) {
3164 /* forward pass */
3165 K->flags |= KALMAN_SMOOTH;
3166 err = kalman_forecast(K, NULL);
3167 K->flags &= ~KALMAN_SMOOTH;
3168 }
3169
3170 K->t = 0;
3171
3172 if (!err) {
3173 if (dist > 1) {
3174 err = koopman_smooth(K, 1);
3175 } else if (dist == 1) {
3176 err = koopman_smooth(K, 0);
3177 } else {
3178 err = anderson_moore_smooth(K);
3179 }
3180 }
3181
3182 bailout:
3183
3184 /* trash the "stepinfo" storage */
3185 free_stepinfo(K);
3186
3187 return err;
3188 }
3189
extract_Q(kalman * K,const gretl_matrix * Sim0)3190 static gretl_matrix *extract_Q (kalman *K,
3191 const gretl_matrix *Sim0)
3192 {
3193 gretl_matrix *Q;
3194 double x;
3195 int i, j;
3196
3197 Q = gretl_matrix_alloc(K->r, K->r);
3198
3199 if (Q != NULL) {
3200 for (i=0; i<K->r; i++) {
3201 for (j=0; j<K->r; j++) {
3202 x = gretl_matrix_get(Sim0, i, j);
3203 gretl_matrix_set(Q, i, j, x);
3204 }
3205 }
3206 }
3207
3208 return Q;
3209 }
3210
3211 /* See the account in Koopman, Shephard and Doornik, Econometrics
3212 Journal, 1999 (volume 2, pp. 113-166), section 4.2, regarding
3213 the initialization of the state under simulation.
3214 */
3215
sim_state_0(kalman * K,const gretl_matrix * U,const gretl_matrix * Sim0)3216 static int sim_state_0 (kalman *K, const gretl_matrix *U,
3217 const gretl_matrix *Sim0)
3218 {
3219 gretl_matrix *Q, *v0 = NULL, *bv = NULL;
3220 int getroot = 1;
3221 int i, err = 0;
3222
3223 if (!kalman_ssfsim(K)) {
3224 if (Sim0 != NULL) {
3225 /* Sim0 contains the state for t = 1 */
3226 err = gretl_matrix_copy_values(K->S0, Sim0);
3227 }
3228 /* error or not, we're done */
3229 return err;
3230 }
3231
3232 /* now we're in the "ssfsim" case, emulating ssfpack */
3233
3234 if (Sim0 != NULL) {
3235 /* Sim0 contains state variance factor
3236 plus the state for t = 0
3237 */
3238 Q = extract_Q(K, Sim0);
3239 getroot = 0;
3240 } else {
3241 Q = gretl_matrix_copy(K->P0);
3242 }
3243
3244 if (Q == NULL) {
3245 err = E_ALLOC;
3246 } else if (getroot) {
3247 err = gretl_matrix_psd_root(Q, 0);
3248 }
3249
3250 if (!err) {
3251 int vlen = K->p > 0 ? K->p : K->r;
3252
3253 v0 = gretl_matrix_alloc(vlen, 1);
3254 if (v0 == NULL) {
3255 err = E_ALLOC;
3256 }
3257 }
3258
3259 if (K->p > 0) {
3260 bv = gretl_matrix_alloc(K->r, 1);
3261 if (bv == NULL) {
3262 err = E_ALLOC;
3263 }
3264 }
3265
3266 if (!err && Sim0 != NULL) {
3267 /* set S0 from last row of Sim0 */
3268 for (i=0; i<K->r; i++) {
3269 K->S0->val[i] = gretl_matrix_get(Sim0, K->r, i);
3270 }
3271 }
3272
3273 if (!err) {
3274 /* handle the t = 0 disturbance */
3275 load_from_row(v0, U, 0, GRETL_MOD_NONE);
3276 if (K->p > 0) {
3277 /* cross-correlated */
3278 gretl_matrix_multiply(K->B, v0, bv);
3279 gretl_matrix_multiply_mod(Q, GRETL_MOD_NONE,
3280 bv, GRETL_MOD_NONE,
3281 K->S0, GRETL_MOD_CUMULATE);
3282 } else {
3283 gretl_matrix_multiply_mod(Q, GRETL_MOD_NONE,
3284 v0, GRETL_MOD_NONE,
3285 K->S0, GRETL_MOD_CUMULATE);
3286 }
3287 }
3288
3289 gretl_matrix_free(Q);
3290 gretl_matrix_free(v0);
3291 gretl_matrix_free(bv);
3292
3293 return err;
3294 }
3295
3296 /* note: it's OK for @S to be NULL (if the simulated
3297 state is not wanted), so watch out for that!
3298 */
3299
kalman_simulate(kalman * K,const gretl_matrix * U,const gretl_matrix * Sim0,gretl_matrix * Y,gretl_matrix * S,PRN * prn)3300 static int kalman_simulate (kalman *K,
3301 const gretl_matrix *U,
3302 const gretl_matrix *Sim0,
3303 gretl_matrix *Y,
3304 gretl_matrix *S,
3305 PRN *prn)
3306 {
3307 gretl_matrix *yt, *et = NULL;
3308 int obs_offset = 0;
3309 int obsdist = 0;
3310 int tmin = 0;
3311 int err = 0;
3312
3313 yt = gretl_zero_matrix_new(K->n, 1);
3314 if (yt == NULL) {
3315 return E_ALLOC;
3316 }
3317
3318 if (K->p > 0) {
3319 et = gretl_matrix_alloc(K->p, 1);
3320 if (et == NULL) {
3321 gretl_matrix_free(yt);
3322 return E_ALLOC;
3323 }
3324 }
3325
3326 if (Y->cols == K->r + K->n) {
3327 /* combined (state, obs) in @Y */
3328 S = Y;
3329 obs_offset = K->r;
3330 }
3331
3332 err = sim_state_0(K, U, Sim0);
3333
3334 if (!err && kalman_ssfsim(K)) {
3335 if (S != NULL) {
3336 load_to_row_offset(S, K->S0, 0, 0);
3337 }
3338 load_to_row_offset(Y, yt, 0, obs_offset);
3339 /* the first row of output is handled */
3340 tmin = 1;
3341 }
3342
3343 if (!err && K->x == NULL) {
3344 /* no exogenous vars */
3345 if (K->A != NULL) {
3346 /* implicit const case: A is 1 x n and A'x is n x 1 */
3347 gretl_vector_copy_values(K->Ax, K->A);
3348 } else {
3349 gretl_matrix_zero(K->Ax);
3350 }
3351 }
3352
3353 if (K->p == 0 && K->R != NULL) {
3354 /* we want to read observation disturbances */
3355 obsdist = 1;
3356 }
3357
3358 for (K->t = tmin; K->t < K->T && !err; K->t += 1) {
3359 int missobs = 0;
3360
3361 if (filter_is_varying(K)) {
3362 err = kalman_refresh_matrices(K, prn);
3363 if (err) {
3364 break;
3365 }
3366 }
3367
3368 /* y_t = A'*x_t + H'*S_t + w_t */
3369 gretl_matrix_multiply_mod(K->H, GRETL_MOD_TRANSPOSE,
3370 K->S0, GRETL_MOD_NONE,
3371 yt, GRETL_MOD_NONE);
3372 if (K->x != NULL) {
3373 kalman_set_Ax(K, &missobs);
3374 }
3375 if (K->A != NULL) {
3376 gretl_matrix_add_to(yt, K->Ax);
3377 }
3378 if (K->p > 0) {
3379 /* C \varepsilon_t */
3380 load_from_row(et, U, K->t, GRETL_MOD_NONE);
3381 gretl_matrix_multiply(K->C, et, K->e);
3382 } else if (obsdist) {
3383 load_from_row_offset(K->e, U, K->t, K->r);
3384 }
3385 gretl_matrix_add_to(yt, K->e);
3386
3387 /* record the t-dated observables */
3388 load_to_row_offset(Y, yt, K->t, obs_offset);
3389
3390 /* record the t-dated state? */
3391 if (S != NULL && tmin == 0) {
3392 load_to_row_offset(S, K->S0, K->t, 0);
3393 }
3394
3395 /* S_{t+1} = F*S_t + v_t */
3396 gretl_matrix_multiply(K->F, K->S0, K->S1);
3397 if (K->p > 0) {
3398 /* B \varepsilon_t */
3399 gretl_matrix_multiply_mod(K->B, GRETL_MOD_NONE,
3400 et, GRETL_MOD_NONE,
3401 K->S1, GRETL_MOD_CUMULATE);
3402 } else {
3403 load_from_row(K->S1, U, K->t, GRETL_MOD_CUMULATE);
3404 }
3405
3406 if (K->mu != NULL) {
3407 gretl_matrix_add_to(K->S1, K->mu);
3408 }
3409
3410 /* record the (t+1)-dated state? */
3411 if (S != NULL && tmin == 1) {
3412 load_to_row_offset(S, K->S1, K->t, 0);
3413 }
3414
3415 gretl_matrix_copy_values(K->S0, K->S1);
3416 }
3417
3418 gretl_matrix_free(yt);
3419 gretl_matrix_free(et);
3420
3421 return err;
3422 }
3423
check_simul_inputs(kalman * K,const gretl_matrix * U,const gretl_matrix * Sim0,const gretl_matrix * SimX,int ssfsim,PRN * prn)3424 static int check_simul_inputs (kalman *K,
3425 const gretl_matrix *U,
3426 const gretl_matrix *Sim0,
3427 const gretl_matrix *SimX,
3428 int ssfsim,
3429 PRN *prn)
3430 {
3431 int err = 0;
3432
3433 if (U == NULL) {
3434 err = missing_matrix_error("U");
3435 } else {
3436 int ncols;
3437
3438 if (K->p > 0) {
3439 /* cross-correlated */
3440 ncols = K->p;
3441 } else {
3442 ncols = K->R == NULL ? K->r : K->r + K->n;
3443 }
3444
3445 if (U->cols != ncols) {
3446 pprintf(prn, "U should have %d columns but has %d\n",
3447 ncols, U->cols);
3448 err = E_NONCONF;
3449 }
3450 }
3451
3452 if (!err && Sim0 != NULL) {
3453 int r = ssfsim ? K->r + 1 : K->r;
3454 int c = ssfsim ? K->r : 1;
3455
3456 if (Sim0->rows != r || Sim0->cols != c) {
3457 pprintf(prn, "simstart should be %d x %d, is %d x %d\n",
3458 r, c, Sim0->rows, Sim0->cols);
3459 }
3460 }
3461
3462 if (!err && K->x != NULL) {
3463 /* do we have enough "obsx" data? */
3464 const gretl_matrix *X = SimX != NULL ? SimX : K->x;
3465
3466 if (X->rows < U->rows) {
3467 pprintf(prn, "obsx should have %d rows but has %d\n",
3468 U->rows, X->rows);
3469 err = E_NONCONF;
3470 }
3471 }
3472
3473 return err;
3474 }
3475
kalman_bundle_simulate(gretl_bundle * b,const gretl_matrix * U,int get_state,PRN * prn,int * err)3476 gretl_matrix *kalman_bundle_simulate (gretl_bundle *b,
3477 const gretl_matrix *U,
3478 int get_state,
3479 PRN *prn, int *err)
3480 {
3481 kalman *K = gretl_bundle_get_private_data(b);
3482 const gretl_matrix *Sim0 = NULL;
3483 const gretl_matrix *SimX = NULL;
3484 gretl_matrix *Ret = NULL;
3485 gretl_matrix *savex = NULL;
3486 double ssfx;
3487 int ssfsim = 0;
3488 int saveT;
3489
3490 if (K == NULL) {
3491 *err = E_DATA;
3492 return NULL;
3493 }
3494
3495 /* try accessing auxiliary info from the bundle */
3496 Sim0 = gretl_bundle_get_matrix(b, "simstart", NULL);
3497 ssfx = gretl_bundle_get_scalar(b, "ssfsim", NULL);
3498 if (K->x != NULL) {
3499 SimX = gretl_bundle_get_matrix(b, "simx", NULL);
3500 }
3501
3502 ssfsim = !na(ssfx) && ssfx != 0;
3503
3504 *err = check_simul_inputs(K, U, Sim0, SimX, ssfsim, prn);
3505 if (*err) {
3506 return NULL;
3507 }
3508
3509 K->b = b; /* attach bundle pointer */
3510
3511 saveT = K->T;
3512 savex = K->x;
3513
3514 /* we let U temporarily define the sample length */
3515 K->T = U->rows;
3516
3517 /* and we allow temporary replacement of K->x */
3518 if (SimX != NULL) {
3519 K->x = (gretl_matrix *) SimX;
3520 }
3521
3522 /* set state */
3523 if (ssfsim) {
3524 K->flags |= (KALMAN_SIM | KALMAN_SSFSIM);
3525 } else {
3526 K->flags |= KALMAN_SIM;
3527 }
3528
3529 /* now, are the other needed matrices in place? */
3530 *err = kalman_bundle_recheck_matrices(K, prn);
3531
3532 /* matrix to hold simulated observables, and state
3533 if wanted */
3534 if (!*err) {
3535 int ncols = get_state ? (K->r + K->n) : K->n;
3536
3537 Ret = gretl_matrix_alloc(K->T, ncols);
3538 if (Ret == NULL) {
3539 *err = E_ALLOC;
3540 }
3541 }
3542
3543 if (!*err) {
3544 *err = kalman_simulate(K, U, Sim0, Ret, NULL, prn);
3545 }
3546
3547 if (*err) {
3548 gretl_matrix_free(Ret);
3549 Ret = NULL;
3550 }
3551
3552 /* restore state */
3553 K->flags &= ~KALMAN_SIM;
3554 K->flags &= ~KALMAN_SSFSIM;
3555 K->T = saveT;
3556 K->x = savex;
3557
3558 return Ret;
3559 }
3560
matrix_is_diagonal(const gretl_matrix * m)3561 static int matrix_is_diagonal (const gretl_matrix *m)
3562 {
3563 double x;
3564 int i, j;
3565
3566 for (j=0; j<m->cols; j++) {
3567 for (i=0; i<m->rows; i++) {
3568 x = gretl_matrix_get(m, i, j);
3569 if (i != j && x != 0.0) return 0;
3570 }
3571 }
3572
3573 return 1;
3574 }
3575
simdata_refresh_QR(kalman * K,PRN * prn)3576 static int simdata_refresh_QR (kalman *K, PRN *prn)
3577 {
3578 int err = 0;
3579
3580 if (matrix_is_varying(K, K_Q) || matrix_is_varying(K, K_R)) {
3581 err = kalman_update_matrices(K, prn);
3582 }
3583
3584 return err;
3585 }
3586
3587 /* Return a matrix in which the standard normal variates
3588 in @U are scaled according to K->Q, and K->R if present.
3589 */
3590
kalman_bundle_simdata(gretl_bundle * b,const gretl_matrix * U,PRN * prn,int * err)3591 gretl_matrix *kalman_bundle_simdata (gretl_bundle *b,
3592 const gretl_matrix *U,
3593 PRN *prn, int *err)
3594 {
3595 kalman *K = gretl_bundle_get_private_data(b);
3596 gretl_matrix *E = NULL;
3597
3598 if (K == NULL || U == NULL) {
3599 *err = E_DATA;
3600 return NULL;
3601 }
3602
3603 if (K->p > 0) {
3604 if (U->cols != K->p) {
3605 *err = E_DATA;
3606 return NULL;
3607 } else {
3608 /* nothing to be done, really */
3609 E = gretl_matrix_copy(U);
3610 }
3611 } else {
3612 int n = K->R == NULL ? 0 : K->n;
3613 int t, j, rn = K->r + n;
3614 int T = U->rows;
3615 int varying = 0;
3616 double vjj, utj;
3617
3618 if (U->cols != rn) {
3619 *err = E_DATA;
3620 return NULL;
3621 }
3622
3623 E = gretl_matrix_alloc(U->rows, rn);
3624 if (E == NULL) {
3625 *err = E_ALLOC;
3626 return NULL;
3627 }
3628
3629 if (matrix_is_varying(K, K_Q) || matrix_is_varying(K, K_R)) {
3630 varying = 1;
3631 }
3632
3633 K->b = b;
3634 set_kalman_running(K);
3635
3636 if (matrix_is_diagonal(K->Q) &&
3637 (K->R == NULL || matrix_is_diagonal(K->R))) {
3638 for (t=0; t<T && !*err; t++) {
3639 if (varying) {
3640 K->t = t;
3641 *err = simdata_refresh_QR(K, prn);
3642 }
3643 for (j=0; j<rn && !*err; j++) {
3644 if (j < K->r) {
3645 vjj = gretl_matrix_get(K->Q, j, j);
3646 } else {
3647 vjj = gretl_matrix_get(K->R, j-K->r, j-K->r);
3648 }
3649 utj = gretl_matrix_get(U, t, j);
3650 gretl_matrix_set(E, t, j, sqrt(vjj) * utj);
3651 }
3652 }
3653 } else {
3654 gretl_matrix *V = gretl_zero_matrix_new(rn, rn);
3655 gretl_matrix *Ut = NULL;
3656 gretl_matrix *Et = NULL;
3657
3658 if (V == NULL) {
3659 *err = E_ALLOC;
3660 goto bailout;
3661 }
3662
3663 if (varying) {
3664 Ut = gretl_matrix_alloc(1, rn);
3665 Et = gretl_matrix_alloc(1, rn);
3666 if (Ut == NULL || Et == NULL) {
3667 gretl_matrix_free(V);
3668 *err = E_ALLOC;
3669 goto bailout;
3670 }
3671 }
3672
3673 if (varying) {
3674 for (t=0; t<T && !*err; t++) {
3675 K->t = t;
3676 *err = simdata_refresh_QR(K, prn);
3677 if (!*err) {
3678 gretl_matrix_inscribe_matrix(V, K->Q, 0, 0,
3679 GRETL_MOD_NONE);
3680 if (n > 0) {
3681 gretl_matrix_inscribe_matrix(V, K->R, K->r, K->r,
3682 GRETL_MOD_NONE);
3683 }
3684 *err = gretl_matrix_psd_root(V, 0);
3685 if (*err) {
3686 gretl_errmsg_set("Failed to compute factor of Omega_t");
3687 } else {
3688 load_from_row(Ut, U, t, GRETL_MOD_NONE);
3689 gretl_matrix_multiply_mod(Ut, GRETL_MOD_NONE,
3690 V, GRETL_MOD_TRANSPOSE,
3691 Et, GRETL_MOD_NONE);
3692 load_to_row(E, Et, t);
3693 }
3694 }
3695 }
3696
3697 gretl_matrix_free(Ut);
3698 gretl_matrix_free(Et);
3699 } else {
3700 gretl_matrix_inscribe_matrix(V, K->Q, 0, 0,
3701 GRETL_MOD_NONE);
3702 if (n > 0) {
3703 gretl_matrix_inscribe_matrix(V, K->R, K->r, K->r,
3704 GRETL_MOD_NONE);
3705 }
3706 *err = gretl_matrix_psd_root(V, 0);
3707 if (*err) {
3708 gretl_errmsg_set("Failed to compute factor of Omega");
3709 } else {
3710 gretl_matrix_multiply_mod(U, GRETL_MOD_NONE,
3711 V, GRETL_MOD_TRANSPOSE,
3712 E, GRETL_MOD_NONE);
3713 }
3714 }
3715
3716 gretl_matrix_free(V);
3717 }
3718 }
3719
3720 bailout:
3721
3722 K->t = 0;
3723 set_kalman_stopped(K);
3724
3725 if (E == NULL && !*err) {
3726 *err = E_ALLOC;
3727 } else if (E != NULL && *err) {
3728 gretl_matrix_free(E);
3729 E = NULL;
3730 }
3731
3732 return E;
3733 }
3734
3735 /* below: functions to support the "wrapping" of a Kalman
3736 struct in a gretl bundle */
3737
check_replacement_dims(const gretl_matrix * orig,const gretl_matrix * repl,int id)3738 static int check_replacement_dims (const gretl_matrix *orig,
3739 const gretl_matrix *repl,
3740 int id)
3741 {
3742 if (id == K_H || id == K_A || id == K_R || id == K_F ||
3743 id == K_Q || id == K_m || id == K_S || id == K_P) {
3744 if (repl->rows != orig->rows ||
3745 repl->cols != orig->cols) {
3746 gretl_errmsg_set("You cannot resize a state-space "
3747 "system matrix");
3748 return E_DATA;
3749 }
3750 } else if (id == K_y || id == K_x) {
3751 if (repl->cols != orig->cols) {
3752 gretl_errmsg_set("You cannot resize a state-space "
3753 "system matrix");
3754 return E_DATA;
3755 }
3756 }
3757
3758 return 0;
3759 }
3760
add_or_replace_k_matrix(gretl_matrix ** targ,gretl_matrix * src,int id,int copy)3761 static int add_or_replace_k_matrix (gretl_matrix **targ,
3762 gretl_matrix *src,
3763 int id, int copy)
3764 {
3765 int err = 0;
3766
3767 if (*targ != src) {
3768 if (*targ != NULL) {
3769 err = check_replacement_dims(*targ, src, id);
3770 if (err) {
3771 return err;
3772 }
3773 /* destroy old Kalman-owned matrix */
3774 gretl_matrix_free(*targ);
3775 }
3776 if (copy) {
3777 *targ = gretl_matrix_copy(src);
3778 if (*targ == NULL) {
3779 err = E_ALLOC;
3780 }
3781 } else {
3782 *targ = src;
3783 }
3784 }
3785
3786 return err;
3787 }
3788
3789 static gretl_matrix **
get_input_matrix_target_by_id(kalman * K,int i)3790 get_input_matrix_target_by_id (kalman *K, int i)
3791 {
3792 gretl_matrix **targ = NULL;
3793
3794 if (i == K_F) {
3795 targ = &K->F;
3796 } else if (i == K_A) {
3797 targ = &K->A;
3798 } else if (i == K_H) {
3799 targ = &K->H;
3800 } else if (i == K_Q) {
3801 if (kalman_xcorr(K)) {
3802 targ = &K->B;
3803 } else {
3804 targ = &K->Q;
3805 }
3806 } else if (i == K_R) {
3807 if (kalman_xcorr(K)) {
3808 targ = &K->C;
3809 } else {
3810 targ = &K->R;
3811 }
3812 } else if (i == K_m) {
3813 targ = &K->mu;
3814 } else if (i == K_y) {
3815 targ = &K->y;
3816 } else if (i == K_x) {
3817 targ = &K->x;
3818 } else if (i == K_S) {
3819 targ = &K->Sini;
3820 } else if (i == K_P) {
3821 targ = &K->Pini;
3822 }
3823
3824 return targ;
3825 }
3826
3827 /* try attaching a matrix to a Kalman bundle */
3828
3829 static int
kalman_bundle_try_set_matrix(kalman * K,void * data,GretlType vtype,int i,int copy)3830 kalman_bundle_try_set_matrix (kalman *K, void *data,
3831 GretlType vtype, int i,
3832 int copy)
3833 {
3834 gretl_matrix **targ;
3835 int err = 0;
3836
3837 targ = get_input_matrix_target_by_id(K, i);
3838
3839 if (targ == NULL) {
3840 err = E_DATA;
3841 } else {
3842 gretl_matrix *m;
3843
3844 if (vtype == GRETL_TYPE_MATRIX) {
3845 m = data;
3846 err = add_or_replace_k_matrix(targ, m, i, copy);
3847 } else if (vtype == GRETL_TYPE_DOUBLE) {
3848 m = gretl_matrix_from_scalar(*(double *) data);
3849 if (m == NULL) {
3850 err = E_ALLOC;
3851 } else {
3852 err = add_or_replace_k_matrix(targ, m, i, 0);
3853 }
3854 } else {
3855 err = E_TYPES;
3856 }
3857 }
3858
3859 return err;
3860 }
3861
3862 static int
kalman_bundle_set_matrix(kalman * K,gretl_matrix * m,int i,int copy)3863 kalman_bundle_set_matrix (kalman *K, gretl_matrix *m, int i, int copy)
3864 {
3865 gretl_matrix **targ;
3866 int err = 0;
3867
3868 targ = get_input_matrix_target_by_id(K, i);
3869
3870 if (targ == NULL) {
3871 err = E_DATA;
3872 } else {
3873 err = add_or_replace_k_matrix(targ, m, i, copy);
3874 }
3875
3876 return err;
3877 }
3878
kalman_output_matrix(kalman * K,const char * key)3879 static gretl_matrix **kalman_output_matrix (kalman *K,
3880 const char *key)
3881 {
3882 gretl_matrix **pm = NULL;
3883
3884 if (!strcmp(key, "prederr")) {
3885 pm = &K->E;
3886 } else if (!strcmp(key, "pevar")) {
3887 pm = &K->V;
3888 } else if (!strcmp(key, "state")) {
3889 pm = &K->S;
3890 } else if (!strcmp(key, "stvar")) {
3891 pm = &K->P;
3892 } else if (!strcmp(key, "gain")) {
3893 pm = &K->K;
3894 } else if (!strcmp(key, "llt")) {
3895 pm = &K->LL;
3896 } else if (!strcmp(key, "smdist")) {
3897 pm = &K->U;
3898 } else if (!strcmp(key, "smdisterr")) {
3899 pm = &K->Vsd;
3900 } else if (!strcmp(key, "uhat")) {
3901 pm = &K->e;
3902 }
3903
3904 return pm;
3905 }
3906
3907 #define K_N_OUTPUTS 9
3908
3909 static const char *kalman_output_matrix_names[K_N_OUTPUTS] = {
3910 "prederr",
3911 "pevar",
3912 "state",
3913 "stvar",
3914 "gain",
3915 "llt",
3916 "smdist",
3917 "smdisterr",
3918 "uhat"
3919 };
3920
output_matrix_slot(const char * s)3921 static int output_matrix_slot (const char *s)
3922 {
3923 int i;
3924
3925 for (i=0; i<K_N_OUTPUTS; i++) {
3926 if (!strcmp(s, kalman_output_matrix_names[i])) {
3927 return i;
3928 }
3929 }
3930
3931 return -1;
3932 }
3933
3934 #define K_N_SCALARS 9
3935
3936 enum {
3937 Ks_t = 0,
3938 Ks_DIFFUSE,
3939 Ks_CROSS,
3940 Ks_S2,
3941 Ks_LNL,
3942 Ks_r,
3943 Ks_n,
3944 Ks_T,
3945 Ks_p,
3946 };
3947
3948 static const char *kalman_output_scalar_names[K_N_SCALARS] = {
3949 "t",
3950 "diffuse",
3951 "cross",
3952 "s2",
3953 "lnl",
3954 "r",
3955 "n",
3956 "T",
3957 "p"
3958 };
3959
kalman_output_scalar(kalman * K,const char * key)3960 static double *kalman_output_scalar (kalman *K,
3961 const char *key)
3962 {
3963 /* static storage for on-the-fly scalars */
3964 static double retval[K_N_SCALARS];
3965 int i, idx = -1;
3966
3967 for (i=0; i<K_N_SCALARS; i++) {
3968 if (!strcmp(key, kalman_output_scalar_names[i])) {
3969 idx = i;
3970 break;
3971 }
3972 }
3973
3974 if (idx < 0) {
3975 return NULL;
3976 }
3977
3978 switch (idx) {
3979 case Ks_t:
3980 if (kalman_is_running(K)) {
3981 retval[idx] = K->t + 1;
3982 } else {
3983 retval[idx] = kalman_checking(K) ? 1 : 0;
3984 }
3985 break;
3986 case Ks_DIFFUSE:
3987 retval[idx] = (K->flags & KALMAN_DIFFUSE)? 1 : 0;
3988 break;
3989 case Ks_CROSS:
3990 retval[idx] = (K->flags & KALMAN_CROSS)? 1 : 0;
3991 break;
3992 case Ks_S2:
3993 retval[idx] = K->s2;
3994 break;
3995 case Ks_LNL:
3996 retval[idx] = K->loglik;
3997 break;
3998 case Ks_r:
3999 retval[idx] = K->r;
4000 break;
4001 case Ks_n:
4002 retval[idx] = K->n;
4003 break;
4004 case Ks_T:
4005 retval[idx] = K->T;
4006 break;
4007 case Ks_p:
4008 retval[idx] = K->p;
4009 break;
4010 default:
4011 break;
4012 }
4013
4014 return &retval[idx];
4015 }
4016
k_input_matrix_by_id(kalman * K,int i)4017 static const gretl_matrix *k_input_matrix_by_id (kalman *K, int i)
4018 {
4019 const gretl_matrix *m = NULL;
4020
4021 if (i == K_F) {
4022 m = K->F;
4023 } else if (i == K_A) {
4024 m = K->A;
4025 } else if (i == K_H) {
4026 m = K->H;
4027 } else if (i == K_Q) {
4028 if (kalman_xcorr(K)) {
4029 m = K->B;
4030 } else {
4031 m = K->Q;
4032 }
4033 } else if (i == K_R) {
4034 if (kalman_xcorr(K)) {
4035 m = K->C;
4036 } else {
4037 m = K->R;
4038 }
4039 } else if (i == K_m) {
4040 m = K->mu;
4041 } else if (i == K_y) {
4042 m = K->y;
4043 } else if (i == K_x) {
4044 m = K->x;
4045 } else if (i == K_S) {
4046 m = K->Sini;
4047 } else if (i == K_P) {
4048 m = K->Pini;
4049 }
4050
4051 return m;
4052 }
4053
input_matrix_slot(const char * s)4054 static int input_matrix_slot (const char *s)
4055 {
4056 int i;
4057
4058 for (i=0; i<K_MMAX; i++) {
4059 if (!strcmp(s, K_input_mats[i].name)) {
4060 return K_input_mats[i].sym;
4061 }
4062 }
4063
4064 return -1;
4065 };
4066
kalman_extra_type(const char * key)4067 static GretlType kalman_extra_type (const char *key)
4068 {
4069 if (!strcmp(key, "ssfsim")) {
4070 return GRETL_TYPE_DOUBLE;
4071 } else if (!strcmp(key, "simstart") ||
4072 !strcmp(key, "simx")) {
4073 return GRETL_TYPE_MATRIX;
4074 } else {
4075 return GRETL_TYPE_NONE;
4076 }
4077 }
4078
maybe_set_kalman_element(void * kptr,const char * key,void * vptr,GretlType vtype,int copy,int * err)4079 int maybe_set_kalman_element (void *kptr,
4080 const char *key,
4081 void *vptr,
4082 GretlType vtype,
4083 int copy,
4084 int *err)
4085 {
4086 GretlType targ;
4087 kalman *K = kptr;
4088 int fncall = 0;
4089 int i, id = -1;
4090 int Kflag = 0;
4091 int done = 0;
4092
4093 if (K == NULL) {
4094 *err = E_DATA;
4095 return 0;
4096 }
4097
4098 /* check for optional "extra" kalman items that
4099 live outside of the kalman struct itself
4100 */
4101 targ = kalman_extra_type(key);
4102 if (targ != GRETL_TYPE_NONE) {
4103 if (vtype != targ) {
4104 *err = E_TYPES;
4105 }
4106 return 0;
4107 }
4108
4109 if (!strcmp(key, "diffuse")) {
4110 Kflag = KALMAN_DIFFUSE;
4111 }
4112
4113 if (Kflag) {
4114 if (vtype == GRETL_TYPE_DOUBLE) {
4115 double x = *(double *) vptr;
4116
4117 if (na(x)) {
4118 *err = E_DATA;
4119 return 0;
4120 } else if (x == 0) {
4121 K->flags &= ~Kflag;
4122 } else {
4123 K->flags |= Kflag;
4124 }
4125 return 1;
4126 } else {
4127 *err = E_TYPES;
4128 return 0;
4129 }
4130 }
4131
4132 if (!strcmp(key, "timevar_call")) {
4133 /* try for a function call specifier (string) */
4134 if (vtype == GRETL_TYPE_STRING) {
4135 fncall = 1;
4136 } else {
4137 *err = E_TYPES;
4138 }
4139 } else {
4140 /* try for a matrix specifier */
4141 i = input_matrix_slot(key);
4142 if (i >= 0) {
4143 if (vtype == GRETL_TYPE_MATRIX ||
4144 vtype == GRETL_TYPE_DOUBLE) {
4145 id = i;
4146 } else {
4147 *err = E_TYPES;
4148 }
4149 }
4150 }
4151
4152 if (*err) {
4153 return 0;
4154 } else if (fncall) {
4155 if (copy) {
4156 K->matcall = gretl_strdup((char *) vptr);
4157 } else {
4158 K->matcall = (char *) vptr;
4159 }
4160 /* re-evaluate what's actually varying */
4161 *err = check_for_matrix_updates(K, NULL);
4162 if (!*err) {
4163 done = 1;
4164 }
4165 } else if (id < 0) {
4166 if (kalman_output_matrix(K, key) != NULL ||
4167 kalman_output_scalar(K, key) != NULL) {
4168 *err = E_DATA;
4169 gretl_errmsg_sprintf("The member %s is read-only", key);
4170 }
4171 } else {
4172 *err = kalman_bundle_try_set_matrix(K, vptr, vtype, id, copy);
4173 done = (*err == 0);
4174 }
4175
4176 return done;
4177 }
4178
maybe_delete_kalman_element(void * kptr,const char * key,int * err)4179 int maybe_delete_kalman_element (void *kptr,
4180 const char *key,
4181 int *err)
4182 {
4183 kalman *K = kptr;
4184 gretl_matrix **pm;
4185 int done = 0;
4186
4187 if (K == NULL) {
4188 return 0;
4189 }
4190
4191 if (kalman_output_scalar(K, key) != NULL ||
4192 input_matrix_slot(key) >= 0 || !strcmp(key, "uhat")) {
4193 /* note: the matrix under the key "uhat" is part of
4194 the internal kalman apparatus */
4195 gretl_errmsg_sprintf("%s: cannot be deleted", key);
4196 *err = E_DATA;
4197 } else if ((pm = kalman_output_matrix(K, key)) != NULL) {
4198 /* OK to delete a user-output matrix */
4199 gretl_matrix_free(*pm);
4200 *pm = NULL;
4201 } else if (!strcmp(key, "timevar_call")) {
4202 /* OK to delete time-variation call */
4203 if (K->matcall != NULL) {
4204 free(K->matcall);
4205 K->matcall = NULL;
4206 free(K->varying);
4207 K->varying = NULL;
4208 done = 1;
4209 } else {
4210 *err = E_DATA;
4211 }
4212 }
4213
4214 return done;
4215 }
4216
maybe_retrieve_kalman_element(void * kptr,const char * key,GretlType * type,int * reserved,int * err)4217 void *maybe_retrieve_kalman_element (void *kptr,
4218 const char *key,
4219 GretlType *type,
4220 int *reserved,
4221 int *err)
4222 {
4223 kalman *K = kptr;
4224 void *ret = NULL;
4225 int i, id = -1;
4226
4227 *type = GRETL_TYPE_NONE;
4228
4229 if (K == NULL) {
4230 *err = E_DATA;
4231 return NULL;
4232 }
4233
4234 if (!strcmp(key, "timevar_call")) {
4235 /* function call specifier? */
4236 *reserved = 1;
4237 if (K->matcall != NULL) {
4238 ret = K->matcall;
4239 *type = GRETL_TYPE_STRING;
4240 }
4241 } else {
4242 /* try for an input matrix specifier */
4243 for (i=0; i<K_MMAX; i++) {
4244 if (!strcmp(key, K_input_mats[i].name)) {
4245 id = K_input_mats[i].sym;
4246 ret = (gretl_matrix *) k_input_matrix_by_id(K, id);
4247 if (ret != NULL) {
4248 *type = GRETL_TYPE_MATRIX;
4249 }
4250 break;
4251 }
4252 }
4253 if (id < 0) {
4254 /* try for an output matrix */
4255 gretl_matrix **pm = kalman_output_matrix(K, key);
4256
4257 if (pm != NULL) {
4258 *reserved = 1;
4259 if (*pm != NULL) {
4260 ret = *pm;
4261 *type = GRETL_TYPE_MATRIX;
4262 }
4263 }
4264 }
4265 if (id < 0 && *reserved == 0) {
4266 /* nothing matched yet: try scalar member */
4267 ret = kalman_output_scalar(K, key);
4268 if (ret != NULL) {
4269 *type = GRETL_TYPE_DOUBLE;
4270 }
4271 }
4272 }
4273
4274 if (id >= 0 && *reserved == 0) {
4275 /* flag the fact that @key was a kalman-reserved
4276 identifier */
4277 *reserved = 1;
4278 }
4279
4280 if (*reserved && ret == NULL) {
4281 gretl_errmsg_sprintf("\"%s\": %s", key, _("no such item"));
4282 *err = E_DATA;
4283 }
4284
4285 return ret;
4286 }
4287
output_matrix_count(kalman * K)4288 static int output_matrix_count (kalman *K)
4289 {
4290 gretl_matrix **pm;
4291 int i, n = 0;
4292
4293 for (i=0; i<K_N_OUTPUTS; i++) {
4294 pm = kalman_output_matrix(K, kalman_output_matrix_names[i]);
4295 n += (pm != NULL && *pm != NULL);
4296 }
4297
4298 return n;
4299 }
4300
print_kalman_bundle_info(void * kptr,PRN * prn)4301 int print_kalman_bundle_info (void *kptr, PRN *prn)
4302 {
4303 kalman *K = kptr;
4304 int err = 0;
4305
4306 if (K == NULL) {
4307 pputs(prn, "Kalman struct: empty\n");
4308 err = E_DATA;
4309 } else {
4310 const gretl_matrix *m;
4311 gretl_matrix **pm;
4312 double *px;
4313 const char *name;
4314 int i, id;
4315
4316 pputs(prn, "\nKalman input matrices\n");
4317
4318 for (i=0; i<K_MMAX; i++) {
4319 id = K_input_mats[i].sym;
4320 m = k_input_matrix_by_id(K, id);
4321 if (m != NULL) {
4322 pprintf(prn, " %s: ", K_input_mats[i].name);
4323 pprintf(prn, "%d x %d\n", m->rows, m->cols);
4324 }
4325 }
4326
4327 if (output_matrix_count(K) > 0) {
4328 pputs(prn, "\nKalman output matrices\n");
4329 for (i=0; i<K_N_OUTPUTS; i++) {
4330 name = kalman_output_matrix_names[i];
4331 pm = kalman_output_matrix(K, name);
4332 if (pm != NULL && *pm != NULL) {
4333 m = *pm;
4334 pprintf(prn, " %s: ", name);
4335 pprintf(prn, "%d x %d\n", m->rows, m->cols);
4336 }
4337 }
4338 }
4339
4340 pputs(prn, "\nKalman scalars\n");
4341
4342 for (i=0; i<K_N_SCALARS; i++) {
4343 name = kalman_output_scalar_names[i];
4344 pprintf(prn, " %s: ", name);
4345 px = kalman_output_scalar(K, name);
4346 if (px == NULL || na(*px)) {
4347 pputs(prn, "NA\n");
4348 } else {
4349 pprintf(prn, "%g\n", *px);
4350 }
4351 }
4352
4353 if (K->matcall != NULL) {
4354 pputs(prn, "\nKalman strings\n");
4355 pprintf(prn, " timevar_call: %s\n", K->matcall);
4356 }
4357 }
4358
4359 return err;
4360 }
4361
4362 /* for use in context of a kalman bundle: serialize the
4363 information in the kalman struct to XML
4364 */
4365
kalman_serialize(void * kptr,PRN * prn)4366 int kalman_serialize (void *kptr, PRN *prn)
4367 {
4368 kalman *K = kptr;
4369 const gretl_matrix *m;
4370 gretl_matrix **pm;
4371 double *px;
4372 const char *name;
4373 int i, err = 0;
4374
4375 if (K == NULL) {
4376 fputs("kalman_serialize: got NULL\n", stderr);
4377 return E_DATA;
4378 }
4379
4380 pputs(prn, "<gretl-kalman>\n");
4381
4382 for (i=0; i<K_MMAX; i++) {
4383 m = k_input_matrix_by_id(K, K_input_mats[i].sym);
4384 if (m != NULL) {
4385 gretl_matrix_serialize(m, K_input_mats[i].name, prn);
4386 }
4387 }
4388
4389 for (i=0; i<K_N_OUTPUTS; i++) {
4390 name = kalman_output_matrix_names[i];
4391 pm = kalman_output_matrix(K, name);
4392 if (pm != NULL && *pm != NULL) {
4393 gretl_matrix_serialize(*pm, name, prn);
4394 }
4395 }
4396
4397 for (i=0; i<=Ks_LNL; i++) {
4398 name = kalman_output_scalar_names[i];
4399 px = kalman_output_scalar(K, name);
4400 if (px != NULL && !na(*px)) {
4401 gretl_finite_scalar_serialize(*px, name, prn);
4402 }
4403 }
4404
4405 if (K->matcall != NULL) {
4406 gretl_string_serialize(K->matcall, "timevar_call", prn);
4407 }
4408
4409 pputs(prn, "</gretl-kalman>\n");
4410
4411 return err;
4412 }
4413
required_matrix_slot(const char * s)4414 static int required_matrix_slot (const char *s)
4415 {
4416 if (!strcmp(s, "obsy")) return 0;
4417 if (!strcmp(s, "obsymat")) return 1;
4418 if (!strcmp(s, "statemat")) return 2;
4419 if (!strcmp(s, "statevar")) return 3;
4420 if (!strcmp(s, "obsvar")) return 4;
4421 return -1;
4422 };
4423
4424 /* for use in context of a kalman bundle: deserialize the
4425 kalman struct from XML
4426 */
4427
kalman_deserialize(void * p1,void * p2,int * err)4428 gretl_bundle *kalman_deserialize (void *p1, void *p2, int *err)
4429 {
4430 xmlNodePtr cur, node = p1;
4431 xmlDocPtr doc = p2;
4432 gretl_matrix *Mreq[5] = {NULL};
4433 gretl_matrix *Mopt[K_MMAX] = {NULL};
4434 gretl_matrix *Mout[K_N_OUTPUTS] = {NULL};
4435 char *tvcall = NULL;
4436 double s2 = NADBL;
4437 double lnl = NADBL;
4438 int copy[5] = {0};
4439 int i, nmats = 0;
4440 int Kflags = 0;
4441 gretl_matrix *m;
4442 double x;
4443 char *key, *strv;
4444 gretl_bundle *b = NULL;
4445
4446 while (node != NULL && !*err) {
4447 if (!xmlStrcmp(node->name, (XUC) "gretl-kalman")) {
4448 cur = node->xmlChildrenNode;
4449 while (cur != NULL && !*err) {
4450 key = (char *) xmlGetProp(cur, (XUC) "name");
4451 if (!xmlStrcmp(cur->name, (XUC) "gretl-matrix")) {
4452 /* pick up kalman matrices */
4453 m = gretl_xml_get_matrix(cur, doc, err);
4454 if ((i = required_matrix_slot(key)) >= 0) {
4455 nmats++;
4456 Mreq[i] = m;
4457 } else if ((i = input_matrix_slot(key)) >= 0) {
4458 Mopt[i] = m;
4459 } else if ((i = output_matrix_slot(key)) >= 0) {
4460 Mout[i] = m;
4461 }
4462 } else if (!xmlStrcmp(cur->name, (XUC) "scalar")) {
4463 /* pick up kalman scalars */
4464 if (gretl_xml_get_prop_as_double(cur, "value", &x)) {
4465 if (!strcmp(key, "diffuse") && x > 0) {
4466 Kflags |= KALMAN_DIFFUSE;
4467 } else if (!strcmp(key, "cross") && x > 0) {
4468 Kflags |= KALMAN_CROSS;
4469 } else if (!strcmp(key, "s2")) {
4470 s2 = x;
4471 } else if (!strcmp(key, "lnl")) {
4472 lnl = x;
4473 }
4474 }
4475 } else if (!xmlStrcmp(cur->name, (XUC) "string")) {
4476 /* pick up kalman strings */
4477 if (!strcmp(key, "timevar_call") &&
4478 gretl_xml_get_prop_as_string(cur, "value", &strv)) {
4479 tvcall = strv;
4480 }
4481 }
4482 free(key);
4483 cur = cur->next;
4484 }
4485 break;
4486 }
4487 node = node->next;
4488 }
4489
4490 if (nmats == 5 && !(Kflags & KALMAN_CROSS)) {
4491 /* drop obsvar from initialization */
4492 Mopt[K_R] = Mreq[4];
4493 Mreq[4] = NULL;
4494 nmats--;
4495 }
4496
4497 if (((Kflags & KALMAN_CROSS) && nmats != 5) ||
4498 (!(Kflags & KALMAN_CROSS) && nmats != 4)) {
4499 *err = E_DATA;
4500 } else {
4501 b = kalman_bundle_new(Mreq, copy, nmats, err);
4502 if (b != NULL) {
4503 kalman *K = gretl_bundle_get_private_data(b);
4504 gretl_matrix **pm;
4505 const char *name;
4506
4507 K->flags = Kflags;
4508 K->s2 = s2;
4509 K->loglik = lnl;
4510
4511 for (i=0; i<K_MMAX; i++) {
4512 if (Mopt[i] != NULL) {
4513 kalman_bundle_set_matrix(K, Mopt[i], i, 0);
4514 }
4515 }
4516 for (i=0; i<K_N_OUTPUTS; i++) {
4517 if (Mout[i] != NULL) {
4518 name = kalman_output_matrix_names[i];
4519 pm = kalman_output_matrix(K, name);
4520 *pm = Mout[i];
4521 }
4522 }
4523 K->matcall = tvcall;
4524 }
4525 }
4526
4527 if (*err) {
4528 /* clean up */
4529 for (i=0; i<5; i++) {
4530 gretl_matrix_free(Mreq[i]);
4531 }
4532 for (i=0; i<K_MMAX; i++) {
4533 gretl_matrix_free(Mopt[i]);
4534 }
4535 for (i=0; i<K_N_OUTPUTS; i++) {
4536 gretl_matrix_free(Mout[i]);
4537 }
4538 free(tvcall);
4539 }
4540
4541 return b;
4542 }
4543
4544 /* Called from gretl_bundle.c to meet the case where the user calls
4545 for a kalman bundle to be copied: here we create a new kalman
4546 struct and copy across the required elements (since they are
4547 not regular bundle members).
4548 */
4549
kalman_bundle_copy(const gretl_bundle * src,int * err)4550 gretl_bundle *kalman_bundle_copy (const gretl_bundle *src, int *err)
4551 {
4552 kalman *K, *Knew;
4553 gretl_bundle *b = NULL;
4554 gretl_matrix *M[5] = {NULL};
4555 gretl_matrix *m, **pm, **pm1;
4556 const char *name;
4557 int copy[5] = {1, 1, 1, 1, 1};
4558 int i, id, k = 4;
4559
4560 K = gretl_bundle_get_private_data((gretl_bundle *) src);
4561
4562 if (K == NULL) {
4563 *err = E_DATA;
4564 return NULL;
4565 }
4566
4567 M[0] = K->y;
4568 M[1] = K->H;
4569 M[2] = K->F;
4570
4571 if (kalman_xcorr(K)) {
4572 M[3] = K->B;
4573 M[4] = K->C;
4574 k = 5;
4575 } else {
4576 M[3] = K->Q;
4577 }
4578
4579 b = kalman_bundle_new(M, copy, k, err);
4580
4581 if (*err) {
4582 return b;
4583 }
4584
4585 Knew = gretl_bundle_get_private_data(b);
4586 Knew->flags = K->flags;
4587
4588 for (i=k; i<K_MMAX && !*err; i++) {
4589 id = K_input_mats[i].sym;
4590 m = (gretl_matrix *) k_input_matrix_by_id(K, id);
4591 if (m != NULL) {
4592 *err = kalman_bundle_set_matrix(Knew, m, i, 1);
4593 }
4594 }
4595
4596 for (i=0; i<K_N_OUTPUTS && !*err; i++) {
4597 name = kalman_output_matrix_names[i];
4598 pm = kalman_output_matrix(K, name);
4599 if (pm != NULL && *pm != NULL) {
4600 pm1 = kalman_output_matrix(Knew, name);
4601 *pm1 = gretl_matrix_copy(*pm);
4602 }
4603 }
4604
4605 Knew->ifc = K->ifc;
4606 Knew->s2 = K->s2;
4607 Knew->loglik = K->loglik;
4608
4609 if (K->flags & KALMAN_CROSS) {
4610 Knew->flags |= KALMAN_CROSS;
4611 }
4612
4613 if (K->flags & KALMAN_DIFFUSE) {
4614 Knew->flags |= KALMAN_DIFFUSE;
4615 }
4616
4617 if (K->matcall != NULL) {
4618 Knew->matcall = gretl_strdup(K->matcall);
4619 }
4620
4621 return b;
4622 }
4623
4624 /* for use in constructing GUI bundle save menu */
4625
kalman_bundle_get_matrix_names(kalman * K,int * ns)4626 char **kalman_bundle_get_matrix_names (kalman *K, int *ns)
4627 {
4628 char **S = NULL;
4629 gretl_matrix **pm;
4630 const char *name;
4631 int i, id, err = 0;
4632
4633 *ns = 0;
4634
4635 for (i=0; i<K_MMAX && !err; i++) {
4636 id = K_input_mats[i].sym;
4637 if (k_input_matrix_by_id(K, id) != NULL) {
4638 err = strings_array_add(&S, ns, K_input_mats[i].name);
4639 }
4640 }
4641
4642 for (i=0; i<K_N_OUTPUTS && !err; i++) {
4643 name = kalman_output_matrix_names[i];
4644 pm = kalman_output_matrix(K, name);
4645 if (pm != NULL && *pm != NULL) {
4646 err = strings_array_add(&S, ns, name);
4647 }
4648 }
4649
4650 return S;
4651 }
4652
4653 /* also for use in constructing GUI bundle save menu */
4654
kalman_bundle_get_scalar_names(kalman * K,int * ns)4655 char **kalman_bundle_get_scalar_names (kalman *K, int *ns)
4656 {
4657 char **S;
4658
4659 *ns = K_N_SCALARS -1 - na(K->s2) - na(K->loglik);
4660 S = strings_array_new(*ns);
4661
4662 if (S != NULL) {
4663 int i = 0;
4664
4665 /* flags */
4666 S[i++] = gretl_strdup("cross");
4667 S[i++] = gretl_strdup("diffuse");
4668
4669 /* actual numerical outputs */
4670 if (!na(K->s2)) {
4671 S[i++] = gretl_strdup("s2");
4672 }
4673 if (!na(K->loglik)) {
4674 S[i++] = gretl_strdup("lnl");
4675 }
4676
4677 /* system dimensions */
4678 S[i++] = gretl_strdup("r");
4679 S[i++] = gretl_strdup("n");
4680 S[i++] = gretl_strdup("T");
4681 S[i++] = gretl_strdup("p");
4682 }
4683
4684 return S;
4685 }
4686
4687 /* to support the nelem() function for kalman bundles */
4688
kalman_bundle_n_members(gretl_bundle * b)4689 int kalman_bundle_n_members (gretl_bundle *b)
4690 {
4691 kalman *K = gretl_bundle_get_private_data(b);
4692 int n = 0;
4693
4694 if (K != NULL) {
4695 int i, id;
4696
4697 n = K_N_SCALARS;
4698
4699 for (i=0; i<K_MMAX; i++) {
4700 id = K_input_mats[i].sym;
4701 n += (k_input_matrix_by_id(K, id) != NULL);
4702 }
4703
4704 n += output_matrix_count(K);
4705 }
4706
4707 return n;
4708 }
4709