1 /* R : A Computer Language for Statistical Data Analysis
2 *
3 * Copyright (C) 1999-2002 The R Core Team
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 2 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, a copy is available at
17 * https://www.R-project.org/Licenses/.
18 */
19
20 #ifdef HAVE_CONFIG_H
21 # include <config.h>
22 #endif
23
24 #include <R.h>
25 #include "ts.h"
26
27
28 #ifndef max
29 #define max(a,b) ((a < b)?(b):(a))
30 #endif
31 #ifndef min
32 #define min(a,b) ((a > b)?(b):(a))
33 #endif
34
35 /* Code in this file based on Applied Statistics algorithms AS154/182
36 (C) Royal Statistical Society 1980, 1982 */
37
38 static void
inclu2(int np,double * xnext,double * xrow,double ynext,double * d,double * rbar,double * thetab)39 inclu2(int np, double *xnext, double *xrow, double ynext,
40 double *d, double *rbar, double *thetab)
41 {
42 double cbar, sbar, di, xi, xk, rbthis, dpi;
43 int i, k, ithisr;
44
45 /* This subroutine updates d, rbar, thetab by the inclusion
46 of xnext and ynext. */
47
48 for (i = 0; i < np; i++) xrow[i] = xnext[i];
49
50 for (ithisr = 0, i = 0; i < np; i++) {
51 if (xrow[i] != 0.0) {
52 xi = xrow[i];
53 di = d[i];
54 dpi = di + xi * xi;
55 d[i] = dpi;
56 cbar = di / dpi;
57 sbar = xi / dpi;
58 for (k = i + 1; k < np; k++) {
59 xk = xrow[k];
60 rbthis = rbar[ithisr];
61 xrow[k] = xk - xi * rbthis;
62 rbar[ithisr++] = cbar * rbthis + sbar * xk;
63 }
64 xk = ynext;
65 ynext = xk - xi * thetab[i];
66 thetab[i] = cbar * thetab[i] + sbar * xk;
67 if (di == 0.0) return;
68 } else ithisr = ithisr + np - i - 1;
69 }
70 }
71
starma(Starma G,int * ifault)72 void starma(Starma G, int *ifault)
73 {
74 int p = G->p, q = G->q, r = G->r, np = G->np, nrbar = G->nrbar;
75 double *phi = G->phi, *theta = G->theta, *a = G->a,
76 *P = G->P, *V = G->V, *thetab = G->thetab, *xnext = G->xnext,
77 *xrow = G->xrow, *rbar = G->rbar;
78 int indi, indj, indn;
79 double phii, phij, ynext, vj, bi;
80 int i, j, k, ithisr, ind, npr, ind1, ind2, npr1, im, jm;
81
82 /* Invoking this subroutine sets the values of v and phi, and
83 obtains the initial values of a and p. */
84
85 /* Check if ar(1) */
86
87 if (!(q > 0 || p > 1)) {
88 V[0] = 1.0;
89 a[0] = 0.0;
90 P[0] = 1.0 / (1.0 - phi[0] * phi[0]);
91 return;
92 }
93
94 /* Check for failure indication. */
95 *ifault = 0;
96 if (p < 0) *ifault = 1;
97 if (q < 0) *ifault += 2;
98 if (p == 0 && q == 0) *ifault = 4;
99 k = q + 1;
100 if (k < p) k = p;
101 if (r != k) *ifault = 5;
102 if (np != r * (r + 1) / 2) *ifault = 6;
103 if (nrbar != np * (np - 1) / 2) *ifault = 7;
104 if (r == 1) *ifault = 8;
105 if (*ifault != 0) return;
106
107 /* Now set a(0), V and phi. */
108
109 for (i = 1; i < r; i++) {
110 a[i] = 0.0;
111 if (i >= p) phi[i] = 0.0;
112 V[i] = 0.0;
113 if (i < q + 1) V[i] = theta[i - 1];
114 }
115 a[0] = 0.0;
116 if (p == 0) phi[0] = 0.0;
117 V[0] = 1.0;
118 ind = r;
119 for (j = 1; j < r; j++) {
120 vj = V[j];
121 for (i = j; i < r; i++) V[ind++] = V[i] * vj;
122 }
123
124 /* Now find p(0). */
125
126 if (p > 0) {
127 /* The set of equations s * vec(p(0)) = vec(v) is solved for
128 vec(p(0)). s is generated row by row in the array xnext. The
129 order of elements in p is changed, so as to bring more leading
130 zeros into the rows of s. */
131
132 for (i = 0; i < nrbar; i++) rbar[i] = 0.0;
133 for (i = 0; i < np; i++) {
134 P[i] = 0.0;
135 thetab[i] = 0.0;
136 xnext[i] = 0.0;
137 }
138 ind = 0;
139 ind1 = -1;
140 npr = np - r;
141 npr1 = npr + 1;
142 indj = npr;
143 ind2 = npr - 1;
144 for (j = 0; j < r; j++) {
145 phij = phi[j];
146 xnext[indj++] = 0.0;
147 indi = npr1 + j;
148 for (i = j; i < r; i++) {
149 ynext = V[ind++];
150 phii = phi[i];
151 if (j != r - 1) {
152 xnext[indj] = -phii;
153 if (i != r - 1) {
154 xnext[indi] -= phij;
155 xnext[++ind1] = -1.0;
156 }
157 }
158 xnext[npr] = -phii * phij;
159 if (++ind2 >= np) ind2 = 0;
160 xnext[ind2] += 1.0;
161 inclu2(np, xnext, xrow, ynext, P, rbar, thetab);
162 xnext[ind2] = 0.0;
163 if (i != r - 1) {
164 xnext[indi++] = 0.0;
165 xnext[ind1] = 0.0;
166 }
167 }
168 }
169
170 ithisr = nrbar - 1;
171 im = np - 1;
172 for (i = 0; i < np; i++) {
173 bi = thetab[im];
174 for (jm = np - 1, j = 0; j < i; j++)
175 bi -= rbar[ithisr--] * P[jm--];
176 P[im--] = bi;
177 }
178
179 /* now re-order p. */
180
181 ind = npr;
182 for (i = 0; i < r; i++) xnext[i] = P[ind++];
183 ind = np - 1;
184 ind1 = npr - 1;
185 for (i = 0; i < npr; i++) P[ind--] = P[ind1--];
186 for (i = 0; i < r; i++) P[i] = xnext[i];
187 } else {
188
189 /* P(0) is obtained by backsubstitution for a moving average process. */
190
191 indn = np;
192 ind = np;
193 for (i = 0; i < r; i++)
194 for (j = 0; j <= i; j++) {
195 --ind;
196 P[ind] = V[ind];
197 if (j != 0) P[ind] += P[--indn];
198 }
199 }
200 }
201
karma(Starma G,double * sumlog,double * ssq,int iupd,int * nit)202 void karma(Starma G, double *sumlog, double *ssq, int iupd, int *nit)
203 {
204 int p = G->p, q = G->q, r = G->r, n = G->n, nu = 0;
205 double *phi = G->phi, *theta = G->theta, *a = G->a, *P = G->P,
206 *V = G->V, *w = G->w, *resid = G->resid, *work = G->xnext;
207
208 int i, j, l, ii, ind, indn, indw;
209 double a1, dt, et, ft, g, ut, phij, phijdt;
210
211 /* Invoking this subroutine updates a, P, sumlog and ssq by inclusion
212 of data values w(1) to w(n). the corresponding values of resid are
213 also obtained. When ft is less than (1 + delta), quick recursions
214 are used. */
215
216 /* for non-zero values of nit, perform quick recursions. */
217
218 if (*nit == 0) {
219 for (i = 0; i < n; i++) {
220
221 /* prediction. */
222
223 if (iupd != 1 || i > 0) {
224
225 /* here dt = ft - 1.0 */
226
227 dt = (r > 1) ? P[r] : 0.0;
228 if (dt < G->delta) goto L610;
229 a1 = a[0];
230 for (j = 0; j < r - 1; j++) a[j] = a[j + 1];
231 a[r - 1] = 0.0;
232 for (j = 0; j < p; j++) a[j] += phi[j] * a1;
233 if(P[0] == 0.0) { /* last obs was available */
234 ind = -1;
235 indn = r;
236 for (j = 0; j < r; j++)
237 for (l = j; l < r; l++) {
238 ++ind;
239 P[ind] = V[ind];
240 if (l < r - 1) P[ind] += P[indn++];
241 }
242 } else {
243 for (j = 0; j < r; j++) work[j] = P[j];
244 ind = -1;
245 indn = r;
246 dt = P[0];
247 for (j = 0; j < r; j++) {
248 phij = phi[j];
249 phijdt = phij * dt;
250 for(l = j; l < r; l++) {
251 ++ind;
252 P[ind] = V[ind] + phi[l] * phijdt;
253 if (j < r - 1) P[ind] += work[j+1] * phi[l];
254 if (l < r - 1)
255 P[ind] += work[l+1] * phij + P[indn++];
256 }
257 }
258 }
259 }
260
261 /* updating. */
262
263 ft = P[0];
264 if(!ISNAN(w[i])) {
265 ut = w[i] - a[0];
266 if (r > 1)
267 for (j = 1, ind = r; j < r; j++) {
268 g = P[j] / ft;
269 a[j] += g * ut;
270 for (l = j; l < r; l++) P[ind++] -= g * P[l];
271 }
272 a[0] = w[i];
273 resid[i] = ut / sqrt(ft);
274 *ssq += ut * ut / ft;
275 *sumlog += log(ft);
276 nu++;
277 for (l = 0; l < r; l++) P[l] = 0.0;
278 } else resid[i] = NA_REAL;
279
280 }
281 *nit = n;
282
283 } else {
284
285 /* quick recursions: never used with missing values */
286
287 i = 0;
288 L610:
289 *nit = i;
290 for (ii = i; ii < n; ii++) {
291 et = w[ii];
292 indw = ii;
293 for (j = 0; j < p; j++) {
294 if (--indw < 0) break;
295 et -= phi[j] * w[indw];
296 }
297 for (j = 0; j < min(ii, q); j++)
298 et -= theta[j] * resid[ii - j - 1];
299 resid[ii] = et;
300 *ssq += et * et;
301 nu++;
302 }
303 }
304 G->nused = nu;
305 }
306
307
308 /* start of AS 182 */
309 void
forkal(Starma G,int d,int il,double * delta,double * y,double * amse,int * ifault)310 forkal(Starma G, int d, int il, double *delta, double *y, double *amse,
311 int *ifault)
312 {
313 int p = G->p, q = G->q, r = G->r, n = G->n, np = G->np;
314 double *phi = G->phi, *V = G->V, *w = G->w, *xrow = G->xrow;
315 double *a, *P, *store;
316 int rd = r + d, rz = rd*(rd + 1)/2;
317 double phii, phij, sigma2, a1, aa, dt, phijdt, ams, tmp;
318 int i, j, k, l, nu = 0;
319 int k1;
320 int i45, jj, kk, lk, ll;
321 int nt;
322 int kk1, lk1;
323 int ind, jkl, kkk;
324 int ind1, ind2;
325
326 /* Finite sample prediction from ARIMA processes. */
327
328 /* This routine will calculate the finite sample predictions
329 and their conditional mean square errors for any ARIMA process. */
330
331 /* invoking this routine will calculate the finite sample predictions */
332 /* and their conditional mean square errors for any arima process. */
333
334 store = (double *) R_alloc(rd, sizeof(double));
335 Free(G->a); G->a = a = Calloc(rd, double);
336 Free(G->P); G->P = P = Calloc(rz, double);
337
338 /* check for input faults. */
339 *ifault = 0;
340 if (p < 0) *ifault = 1;
341 if (q < 0) *ifault += 2;
342 if (p * p + q * q == 0) *ifault = 4;
343 if (r != max(p, q + 1)) *ifault = 5;
344 if (np != r * (r + 1) / 2) *ifault = 6;
345 if (d < 0) *ifault = 8;
346 if (il < 1) *ifault = 11;
347 if (*ifault != 0) return;
348
349 /* Find initial likelihood conditions. */
350
351 if (r == 1) {
352 a[0] = 0.0;
353 V[0] = 1.0;
354 P[0] = 1.0 / (1.0 - phi[0] * phi[0]);
355 } else starma(G, ifault);
356
357 /* Calculate data transformations */
358
359 nt = n - d;
360 if (d > 0) {
361 for (j = 0; j < d; j++) {
362 store[j] = w[n - j - 2];
363 if(ISNAN(store[j]))
364 error(_("missing value in last %d observations"), d);
365 }
366 for (i = 0; i < nt; i++) {
367 aa = 0.0;
368 for (k = 0; k < d; ++k) aa -= delta[k] * w[d + i - k - 1];
369 w[i] = w[i + d] + aa;
370 }
371 }
372
373 /* Evaluate likelihood to obtain final Kalman filter conditions */
374
375 {
376 double sumlog = 0.0, ssq = 0.0;
377 int nit = 0;
378 G->n = nt;
379 karma(G, &sumlog, &ssq, 1, &nit);
380 }
381
382
383 /* Calculate m.l.e. of sigma squared */
384
385 sigma2 = 0.0;
386 for (j = 0; j < nt; j++) {
387 /* macOS/gcc 3.5 didn't have isnan defined properly */
388 tmp = G->resid[j];
389 if(!ISNAN(tmp)) { nu++; sigma2 += tmp * tmp; }
390 }
391
392 sigma2 /= nu;
393
394 /* reset the initial a and P when differencing occurs */
395
396 if (d > 0) {
397 for (i = 0; i < np; i++) xrow[i] = P[i];
398 for (i = 0; i < rz; i++) P[i] = 0.0;
399 ind = 0;
400 for (j = 0; j < r; j++) {
401 k = j * (rd + 1) - j * (j + 1) / 2;
402 for (i = j; i < r; i++) P[k++] = xrow[ind++];
403 }
404 for (j = 0; j < d; j++) a[r + j] = store[j];
405 }
406
407 i45 = 2*rd + 1;
408 jkl = r * (2*d + r + 1) / 2;
409
410 for (l = 0; l < il; ++l) {
411
412 /* predict a */
413
414 a1 = a[0];
415 for (i = 0; i < r - 1; i++) a[i] = a[i + 1];
416 a[r - 1] = 0.0;
417 for (j = 0; j < p; j++) a[j] += phi[j] * a1;
418 if (d > 0) {
419 for (j = 0; j < d; j++) a1 += delta[j] * a[r + j];
420 for (i = rd - 1; i > r; i--) a[i] = a[i - 1];
421 a[r] = a1;
422 }
423
424 /* predict P */
425
426 if (d > 0) {
427 for (i = 0; i < d; i++) {
428 store[i] = 0.0;
429 for (j = 0; j < d; j++) {
430 ll = max(i, j);
431 k = min(i, j);
432 jj = jkl + (ll - k) + k * (2*d + 2 - k - 1) / 2;
433 store[i] += delta[j] * P[jj];
434 }
435 }
436 if (d > 1) {
437 for (j = 0; j < d - 1; j++) {
438 jj = d - j - 1;
439 lk = (jj - 1) * (2*d + 2 - jj) / 2 + jkl;
440 lk1 = jj * (2*d + 1 - jj) / 2 + jkl;
441 for (i = 0; i <= j; i++) P[lk1++] = P[lk++];
442 }
443 for (j = 0; j < d - 1; j++)
444 P[jkl + j + 1] = store[j] + P[r + j];
445 }
446 P[jkl] = P[0];
447 for (i = 0; i < d; i++)
448 P[jkl] += delta[i] * (store[i] + 2.0 * P[r + i]);
449 for (i = 0; i < d; i++) store[i] = P[r + i];
450 for (j = 0; j < r; j++) {
451 kk1 = (j+1) * (2*rd - j - 2) / 2 + r;
452 k1 = j * (2*rd - j - 1) / 2 + r;
453 for (i = 0; i < d; i++) {
454 kk = kk1 + i;
455 k = k1 + i;
456 P[k] = phi[j] * store[i];
457 if (j < r - 1) P[k] += P[kk];
458 }
459 }
460
461 for (j = 0; j < r; j++) {
462 store[j] = 0.0;
463 kkk = (j + 1) * (i45 - j - 1) / 2 - d;
464 for (i = 0; i < d; i++) store[j] += delta[i] * P[kkk++];
465 }
466 for (j = 0; j < r; j++) {
467 k = (j + 1) * (rd + 1) - (j + 1) * (j + 2) / 2;
468 for (i = 0; i < d - 1; i++) {
469 --k;
470 P[k] = P[k - 1];
471 }
472 }
473 for (j = 0; j < r; j++) {
474 k = j * (2*rd - j - 1) / 2 + r;
475 P[k] = store[j] + phi[j] * P[0];
476 if (j < r - 1) P[k] += P[j + 1];
477 }
478 }
479 for (i = 0; i < r; i++) store[i] = P[i];
480
481 ind = 0;
482 dt = P[0];
483 for (j = 0; j < r; j++) {
484 phij = phi[j];
485 phijdt = phij * dt;
486 ind2 = j * (2*rd - j + 1) / 2 - 1;
487 ind1 = (j + 1) * (i45 - j - 1) / 2 - 1;
488 for (i = j; i < r; i++) {
489 ++ind2;
490 phii = phi[i];
491 P[ind2] = V[ind++] + phii * phijdt;
492 if (j < r - 1) P[ind2] += store[j + 1] * phii;
493 if (i < r - 1)
494 P[ind2] += store[i + 1] * phij + P[++ind1];
495 }
496 }
497
498 /* predict y */
499
500 y[l] = a[0];
501 for (j = 0; j < d; j++) y[l] += a[r + j] * delta[j];
502
503 /* calculate m.s.e. of y */
504
505 ams = P[0];
506 if (d > 0) {
507 for (j = 0; j < d; j++) {
508 k = r * (i45 - r) / 2 + j * (2*d + 1 - j) / 2;
509 tmp = delta[j];
510 ams += 2.0 * tmp * P[r + j] + P[k] * tmp * tmp;
511 }
512 for (j = 0; j < d - 1; j++) {
513 k = r * (i45 - r) / 2 + 1 + j * (2*d + 1 - j) / 2;
514 for (i = j + 1; i < d; i++)
515 ams += 2.0 * delta[i] * delta[j] * P[k++];
516 }
517 }
518 amse[l] = ams * sigma2;
519 }
520 return;
521 }
522