1 /*
2 * Copyright (C) 1999 Martyn Plummer
3 * Copyright (C) 1999-2016 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 #include <math.h>
21 #include <string.h>
22 #include <R.h>
23 #include <R_ext/Applic.h> /* Fortran routines */
24 #include "ts.h"
25 #include "stats.h"
26
27
28 #define MAX_DIM_LENGTH 4
29
30 #define VECTOR(x) (x.vec)
31 #define MATRIX(x) (x.mat)
32 #define ARRAY1(x) (x.vec)
33 #define ARRAY2(x) (x.mat)
34 #define ARRAY3(x) (x.arr3)
35 #define ARRAY4(x) (x.arr4)
36 #define DIM(x) (x.dim)
37 #define NROW(x) (x.dim[0])
38 #define NCOL(x) (x.dim[1])
39 #define DIM_LENGTH(x) (x.ndim)
40
41
42 typedef struct array {
43 double *vec;
44 double **mat;
45 double ***arr3;
46 double ****arr4;
47 int dim[MAX_DIM_LENGTH];
48 int ndim;
49 } Array;
50
51 static Array make_array(double vec[], int dim[], int ndim);
52 static Array make_zero_array(int dim[], int ndim);
53 static Array make_matrix(double vec[], int nrow, int ncol);
54 static Array make_zero_matrix(int nrow, int ncol);
55 static Array make_identity_matrix(int n);
56
57 static Array subarray(Array a, int index);
58
59 static int vector_length(Array a);
60
61 static void set_array_to_zero(Array arr);
62 static void copy_array (Array orig, Array ans);
63 static void array_op(Array arr1, Array arr2, char op, Array ans);
64 static void scalar_op(Array arr, double s, char op, Array ans);
65
66 static void transpose_matrix(Array mat, Array ans);
67 static void matrix_prod(Array mat1, Array mat2, int trans1, int trans2,
68 Array ans);
69
70
71 /* Functions for dynamically allocating arrays
72
73 The Array structure contains pointers to arrays which are allocated
74 using the R_alloc function. Although the .C() interface cleans up
75 all memory assigned with R_alloc, judicious use of vmaxget() vmaxset()
76 to free this memory is probably wise. See memory.c in R core.
77
78 */
79
assert(int bool)80 static void assert(int bool)
81 {
82 if(!bool)
83 error(_("assert failed in src/library/ts/src/carray.c"));
84 }
85
init_array(void)86 static Array init_array(void)
87 {
88 int i;
89 Array a;
90
91 /* Initialize everything to zero. Useful for debugging */
92 ARRAY1(a) = (double *) '\0';
93 ARRAY2(a) = (double **) '\0';
94 ARRAY3(a) = (double ***) '\0';
95 ARRAY4(a) = (double ****) '\0';
96 for (i = 0; i < MAX_DIM_LENGTH; i++)
97 DIM(a)[i] = 0;
98 DIM_LENGTH(a) = 0;
99
100 return a;
101 }
102
vector_length(Array a)103 static int vector_length(Array a)
104 {
105 int i, len;
106
107 for (i = 0, len = 1; i < DIM_LENGTH(a); i++) {
108 len *= DIM(a)[i];
109 }
110
111 return len;
112 }
113
114
make_array(double vec[],int dim[],int ndim)115 static Array make_array(double vec[], int dim[], int ndim)
116 {
117 int d, i, j;
118 int len[MAX_DIM_LENGTH + 1];
119 Array a;
120
121 assert(ndim <= MAX_DIM_LENGTH);
122
123 a = init_array();
124
125 len[ndim] = 1;
126 for (d = ndim; d >= 1; d--) {
127 len[d-1] = len[d] * dim[ndim - d];
128 }
129
130 for (d = 1; d <= ndim; d++) {
131 switch(d) {
132 case 1:
133 VECTOR(a) = vec;
134 break;
135 case 2:
136 ARRAY2(a) = (double**) R_alloc(len[2 - 1],sizeof(double*));
137 for(i = 0, j = 0; i < len[2 - 1]; i++, j+=dim[ndim - 2 + 1]) {
138 ARRAY2(a)[i] = ARRAY1(a) + j;
139 }
140 break;
141 case 3:
142 ARRAY3(a) = (double***) R_alloc(len[3 - 1],sizeof(double**));
143 for(i = 0, j = 0; i < len[3 - 1]; i++, j+=dim[ndim - 3 + 1]) {
144 ARRAY3(a)[i] = ARRAY2(a) + j;
145 }
146 break;
147 case 4:
148 ARRAY4(a) = (double****) R_alloc(len[4 - 1],sizeof(double***));
149 for(i = 0, j = 0; i < len[4 - 1]; i++, j+=dim[ndim - 4 + 1]) {
150 ARRAY4(a)[i] = ARRAY3(a) + j;
151 }
152 break;
153 default:
154 break;
155 }
156 }
157
158 for (i = 0; i < ndim; i++) {
159 DIM(a)[i] = dim[i];
160 }
161 DIM_LENGTH(a) = ndim;
162
163 return a;
164 }
165
make_zero_array(int dim[],int ndim)166 static Array make_zero_array(int dim[], int ndim)
167 {
168 int i;
169 int len;
170 double *vec;
171
172 for (i = 0, len = 1; i < ndim; i++) {
173 len *= dim[i];
174 }
175
176 vec = (double *) R_alloc(len, sizeof(double));
177 for (i = 0; i < len; i++) {
178 vec[i] = 0.0;
179 }
180
181 return make_array(vec, dim, ndim);
182
183 }
184
make_matrix(double vec[],int nrow,int ncol)185 static Array make_matrix(double vec[], int nrow, int ncol)
186 {
187 int dim[2];
188
189 dim[0] = nrow;
190 dim[1] = ncol;
191 return make_array(vec, dim, 2);
192 }
193
make_zero_matrix(int nrow,int ncol)194 static Array make_zero_matrix(int nrow, int ncol)
195 {
196 int dim[2];
197 Array a;
198
199 dim[0] = nrow;
200 dim[1] = ncol;
201 a = make_zero_array(dim, 2);
202 return a;
203 }
204
subarray(Array a,int index)205 static Array subarray(Array a, int index)
206 /* Return subarray of array a in the form of an Array
207 structure so it can be manipulated by other functions
208 NB The data are not copied, so any changes made to the
209 subarray will affect the original array.
210 */
211 {
212 int i, offset;
213 Array b;
214
215 b = init_array();
216
217 /* is index in range? */
218 assert( index >= 0 && index < DIM(a)[0] );
219
220 offset = index;
221 switch(DIM_LENGTH(a)) {
222 /* NB Falling through here */
223 case 4:
224 offset *= DIM(a)[DIM_LENGTH(a) - 4 + 1];
225 ARRAY3(b) = ARRAY3(a) + offset;
226 case 3:
227 offset *= DIM(a)[DIM_LENGTH(a) - 3 + 1];
228 ARRAY2(b) = ARRAY2(a) + offset;
229 case 2:
230 offset *= DIM(a)[DIM_LENGTH(a) - 2 + 1];
231 ARRAY1(b) = ARRAY1(a) + offset;
232 break;
233 default:
234 break;
235 }
236
237
238 DIM_LENGTH(b) = DIM_LENGTH(a) - 1;
239
240 for (i = 0; i < DIM_LENGTH(b); i++)
241 DIM(b)[i] = DIM(a)[i+1];
242
243 return b;
244
245 }
246
test_array_conform(Array a1,Array a2)247 static int test_array_conform(Array a1, Array a2)
248 {
249 int i, ans = FALSE;
250
251 if (DIM_LENGTH(a1) != DIM_LENGTH(a2))
252 return FALSE;
253 else
254 for (i = 0; i < DIM_LENGTH(a1); i++) {
255 if (DIM(a1)[i] == DIM(a2)[i])
256 ans = TRUE;
257 else
258 return FALSE;
259 }
260
261 return ans;
262 }
263
copy_array(Array orig,Array ans)264 static void copy_array (Array orig, Array ans)
265 /* copy matrix orig to ans */
266 {
267 int i;
268
269 assert (test_array_conform(orig, ans));
270
271 for(i = 0; i < vector_length(orig); i++)
272 VECTOR(ans)[i] = VECTOR(orig)[i];
273 }
274
transpose_matrix(Array mat,Array ans)275 static void transpose_matrix(Array mat, Array ans)
276 {
277 int i,j;
278 const void *vmax;
279 Array tmp;
280
281 tmp = init_array();
282
283 assert(DIM_LENGTH(mat) == 2 && DIM_LENGTH(ans) == 2);
284 assert(NCOL(mat) == NROW(ans));
285 assert(NROW(mat) == NCOL(ans));
286
287 vmax = vmaxget();
288
289 tmp = make_zero_matrix(NROW(ans), NCOL(ans));
290 for(i = 0; i < NROW(mat); i++)
291 for(j = 0; j < NCOL(mat); j++)
292 MATRIX(tmp)[j][i] = MATRIX(mat)[i][j];
293 copy_array(tmp, ans);
294
295 vmaxset(vmax);
296 }
297
array_op(Array arr1,Array arr2,char op,Array ans)298 static void array_op(Array arr1, Array arr2, char op, Array ans)
299 /* Element-wise array operations */
300 {
301 int i;
302
303 assert (test_array_conform(arr1, arr2));
304 assert (test_array_conform(arr2, ans));
305
306 switch (op) {
307 case '*':
308 for (i = 0; i < vector_length(ans); i++)
309 VECTOR(ans)[i] = VECTOR(arr1)[i] * VECTOR(arr2)[i];
310 break;
311 case '+':
312 for (i = 0; i < vector_length(ans); i++)
313 VECTOR(ans)[i] = VECTOR(arr1)[i] + VECTOR(arr2)[i];
314 break;
315 case '/':
316 for (i = 0; i < vector_length(ans); i++)
317 VECTOR(ans)[i] = VECTOR(arr1)[i] / VECTOR(arr2)[i];
318 break;
319 case '-':
320 for (i = 0; i < vector_length(ans); i++)
321 VECTOR(ans)[i] = VECTOR(arr1)[i] - VECTOR(arr2)[i];
322 break;
323 default:
324 printf("Unknown op in array_op");
325 }
326 }
327
328
scalar_op(Array arr,double s,char op,Array ans)329 static void scalar_op(Array arr, double s, char op, Array ans)
330 /* Elementwise scalar operations */
331 {
332 int i;
333
334 assert (test_array_conform(arr, ans));
335
336 switch (op) {
337 case '*':
338 for (i = 0; i < vector_length(ans); i++)
339 VECTOR(ans)[i] = VECTOR(arr)[i] * s;
340 break;
341 case '+':
342 for (i = 0; i < vector_length(ans); i++)
343 VECTOR(ans)[i] = VECTOR(arr)[i] + s;
344 break;
345 case '/':
346 for (i = 0; i < vector_length(ans); i++)
347 VECTOR(ans)[i] = VECTOR(arr)[i] / s;
348 break;
349 case '-':
350 for (i = 0; i < vector_length(ans); i++)
351 VECTOR(ans)[i] = VECTOR(arr)[i] - s;
352 break;
353 default:
354 printf("Unknown op in array_op");
355 }
356 }
357
matrix_prod(Array mat1,Array mat2,int trans1,int trans2,Array ans)358 static void matrix_prod(Array mat1, Array mat2, int trans1, int trans2, Array ans)
359 /*
360 General matrix product between mat1 and mat2. Put answer in ans.
361 trans1 and trans2 are logical flags which indicate if the matrix is
362 to be transposed. Normal matrix multiplication has trans1 = trans2 = 0.
363 */
364 {
365 int i,j,k,K1,K2;
366 const void *vmax;
367 double m1, m2;
368 Array tmp;
369
370 /* Test whether everything is a matrix */
371 assert(DIM_LENGTH(mat1) == 2 &&
372 DIM_LENGTH(mat2) == 2 && DIM_LENGTH(ans) == 2);
373
374 /* Test whether matrices conform. K is the dimension that is
375 lost by multiplication */
376 if (trans1) {
377 assert ( NCOL(mat1) == NROW(ans) );
378 K1 = NROW(mat1);
379 }
380 else {
381 assert ( NROW(mat1) == NROW(ans) );
382 K1 = NCOL(mat1);
383 }
384 if (trans2) {
385 assert ( NROW(mat2) == NCOL(ans) );
386 K2 = NCOL(mat2);
387 }
388 else {
389 assert ( NCOL(mat2) == NCOL(ans) );
390 K2 = NROW(mat2);
391 }
392 assert (K1 == K2);
393
394 tmp = init_array();
395
396 /* In case ans is the same as mat1 or mat2, we create a temporary
397 matrix to hold the answer, then copy it to ans
398 */
399 vmax = vmaxget();
400
401 tmp = make_zero_matrix(NROW(ans), NCOL(ans));
402 for (i = 0; i < NROW(tmp); i++) {
403 for (j = 0; j < NCOL(tmp); j++) {
404 for(k = 0; k < K1; k++) {
405 m1 = (trans1) ? MATRIX(mat1)[k][i] : MATRIX(mat1)[i][k];
406 m2 = (trans2) ? MATRIX(mat2)[j][k] : MATRIX(mat2)[k][j];
407 MATRIX(tmp)[i][j] += m1 * m2;
408 }
409 }
410 }
411 copy_array(tmp, ans);
412
413 vmaxset(vmax);
414 }
415
set_array_to_zero(Array arr)416 static void set_array_to_zero(Array arr)
417 {
418 int i;
419
420 for (i = 0; i < vector_length(arr); i++)
421 VECTOR(arr)[i] = 0.0;
422 }
423
make_identity_matrix(int n)424 static Array make_identity_matrix(int n)
425 {
426 int i;
427 Array a;
428
429 a = make_zero_matrix(n,n);
430 for(i = 0; i < n; i++)
431 MATRIX(a)[i][i] = 1.0;
432
433 return a;
434 }
435
qr_solve(Array x,Array y,Array coef)436 static void qr_solve(Array x, Array y, Array coef)
437 /* Translation of the R function qr.solve into pure C
438 NB We have to transpose the matrices since the ordering of an array is different in Fortran
439 NB2 We have to copy x to avoid it being overwritten.
440 */
441 {
442 int i, info = 0, rank, *pivot, n, p;
443 const void *vmax;
444 double tol = 1.0E-7, *qraux, *work;
445 Array xt, yt, coeft;
446
447 assert(NROW(x) == NROW(y));
448 assert(NCOL(coef) == NCOL(y));
449 assert(NCOL(x) == NROW(coef));
450
451 vmax = vmaxget();
452
453 qraux = (double *) R_alloc(NCOL(x), sizeof(double));
454 pivot = (int *) R_alloc(NCOL(x), sizeof(int));
455 work = (double *) R_alloc(2*NCOL(x), sizeof(double));
456
457 for(i = 0; i < NCOL(x); i++)
458 pivot[i] = i+1;
459
460 xt = make_zero_matrix(NCOL(x), NROW(x));
461 transpose_matrix(x,xt);
462
463 n = NROW(x);
464 p = NCOL(x);
465
466 F77_CALL(dqrdc2)(VECTOR(xt), &n, &n, &p, &tol, &rank,
467 qraux, pivot, work);
468
469 if (rank != p)
470 error(_("Singular matrix in qr_solve"));
471
472 yt = make_zero_matrix(NCOL(y), NROW(y));
473 coeft = make_zero_matrix(NCOL(coef), NROW(coef));
474 transpose_matrix(y, yt);
475
476 F77_CALL(dqrcf)(VECTOR(xt), &NROW(x), &rank, qraux,
477 yt.vec, &NCOL(y), coeft.vec, &info);
478
479 transpose_matrix(coeft,coef);
480
481 vmaxset(vmax);
482 }
483
ldet(Array x)484 static double ldet(Array x)
485 /* Log determinant of square matrix */
486 {
487 int i, rank, *pivot, n, p;
488 const void *vmax;
489 double ll, tol = 1.0E-7, *qraux, *work;
490 Array xtmp;
491
492 assert(DIM_LENGTH(x) == 2); /* is x a matrix? */
493 assert(NROW(x) == NCOL(x)); /* is x square? */
494
495 vmax = vmaxget();
496
497 qraux = (double *) R_alloc(NCOL(x), sizeof(double));
498 pivot = (int *) R_alloc(NCOL(x), sizeof(int));
499 work = (double *) R_alloc(2*NCOL(x), sizeof(double));
500
501 xtmp = make_zero_matrix(NROW(x), NCOL(x));
502 copy_array(x, xtmp);
503
504 for(i = 0; i < NCOL(x); i++)
505 pivot[i] = i+1;
506
507 p = n = NROW(x);
508
509 F77_CALL(dqrdc2)(VECTOR(xtmp), &n, &n, &p, &tol, &rank,
510 qraux, pivot, work);
511
512 if (rank != p)
513 error(_("Singular matrix in ldet"));
514
515 for (i = 0, ll=0.0; i < rank; i++) {
516 ll += log(fabs(MATRIX(xtmp)[i][i]));
517 }
518
519 vmaxset(vmax);
520
521 return ll;
522 }
523
524
525
526 /* Burg's algorithm for autoregression estimation
527
528 multi_burg is the interface to R. It also handles model selection
529 using AIC
530
531 burg implements the main part of the algorithm
532
533 burg2 estimates the partial correlation coefficient. This
534 requires iteration in the multivariate case.
535
536 Notation
537
538 resid residuals (forward and backward)
539 A Estimates of autocorrelation coefficients
540 V Prediction Variance
541 K Partial correlation coefficient
542 */
543
544
545 #define BURG_MAX_ITER 20
546 #define BURG_TOL 1.0E-8
547
548 void multi_burg(int *pn, double *x, int *pomax, int *pnser, double *coef,
549 double *pacf, double *var, double *aic, int *porder, int *useaic,
550 int *vmethod);
551 static void burg0(int omax, Array resid_f, Array resid_b, Array *A, Array *B,
552 Array P, Array V, int vmethod);
553 static void burg2(Array ss_ff, Array ss_bb, Array ss_fb, Array E,
554 Array KA, Array KB);
555
multi_burg(int * pn,double * x,int * pomax,int * pnser,double * coef,double * pacf,double * var,double * aic,int * porder,int * useaic,int * vmethod)556 void multi_burg(int *pn, double *x, int *pomax, int *pnser, double *coef,
557 double *pacf, double *var, double *aic, int *porder, int *useaic,
558 int *vmethod)
559 {
560 int i, j, m, omax = *pomax, n = *pn, nser=*pnser, order=*porder;
561 int dim1[3];
562 double aicmin;
563 Array xarr, resid_f, resid_b, resid_f_tmp;
564 Array *A, *B, P, V;
565
566 dim1[0] = omax+1; dim1[1] = dim1[2] = nser;
567 A = (Array *) R_alloc(omax+1, sizeof(Array));
568 B = (Array *) R_alloc(omax+1, sizeof(Array));
569 for (i = 0; i <= omax; i++) {
570 A[i] = make_zero_array(dim1, 3);
571 B[i] = make_zero_array(dim1, 3);
572 }
573 P = make_array(pacf, dim1, 3);
574 V = make_array(var, dim1, 3);
575
576 xarr = make_matrix(x, nser, n);
577 resid_f = make_zero_matrix(nser, n);
578 resid_b = make_zero_matrix(nser, n);
579 set_array_to_zero(resid_b);
580 copy_array(xarr, resid_f);
581 copy_array(xarr, resid_b);
582 resid_f_tmp = make_zero_matrix(nser, n);
583
584 burg0(omax, resid_f, resid_b, A, B, P, V, *vmethod);
585
586 /* Model order selection */
587
588 for (i = 0; i <= omax; i++) {
589 aic[i] = n * ldet(subarray(V,i)) + 2 * i * nser * nser;
590 }
591 if (*useaic) {
592 order = 0;
593 aicmin = aic[0];
594 for (i = 1; i <= omax; i++) {
595 if (aic[i] < aicmin) {
596 aicmin = aic[i];
597 order = i;
598 }
599 }
600 }
601 else order = omax;
602 *porder = order;
603
604 for(i = 0; i < vector_length(A[order]); i++)
605 coef[i] = VECTOR(A[order])[i];
606
607 if (*useaic) {
608 /* Recalculate residuals for chosen model */
609 set_array_to_zero(resid_f);
610 set_array_to_zero(resid_f_tmp);
611 for (m = 0; m <= order; m++) {
612 for (i = 0; i < NROW(resid_f_tmp); i++) {
613 for (j = 0; j < NCOL(resid_f_tmp) - order; j++) {
614 MATRIX(resid_f_tmp)[i][j + order] = MATRIX(xarr)[i][j + order - m];
615 }
616 }
617 matrix_prod(subarray(A[order],m), resid_f_tmp, 0, 0, resid_f_tmp);
618 array_op(resid_f_tmp, resid_f, '+', resid_f);
619 }
620 }
621 copy_array(resid_f, xarr);
622
623 }
624
625
burg0(int omax,Array resid_f,Array resid_b,Array * A,Array * B,Array P,Array V,int vmethod)626 static void burg0(int omax, Array resid_f, Array resid_b, Array *A, Array *B,
627 Array P, Array V, int vmethod)
628 {
629 int i, j, m, n = NCOL(resid_f), nser=NROW(resid_f);
630 Array ss_ff, ss_bb, ss_fb;
631 Array resid_f_tmp, resid_b_tmp;
632 Array KA, KB, E;
633 Array id, tmp;
634
635 ss_ff = make_zero_matrix(nser, nser);
636 ss_fb = make_zero_matrix(nser, nser);
637 ss_bb = make_zero_matrix(nser, nser);
638
639 resid_f_tmp = make_zero_matrix(nser, n);
640 resid_b_tmp = make_zero_matrix(nser, n);
641
642 id = make_identity_matrix(nser);
643
644 tmp = make_zero_matrix(nser, nser);
645
646 E = make_zero_matrix(nser, nser);
647 KA = make_zero_matrix(nser, nser);
648 KB = make_zero_matrix(nser, nser);
649
650 set_array_to_zero(A[0]);
651 set_array_to_zero(B[0]);
652 copy_array(id, subarray(A[0],0));
653 copy_array(id, subarray(B[0],0));
654
655 matrix_prod(resid_f, resid_f, 0, 1, E);
656 scalar_op(E, n, '/', E);
657 copy_array(E, subarray(V,0));
658
659 for (m = 0; m < omax; m++) {
660
661 for(i = 0; i < nser; i++) {
662 for (j = n - 1; j > m; j--) {
663 MATRIX(resid_b)[i][j] = MATRIX(resid_b)[i][j-1];
664 }
665 MATRIX(resid_f)[i][m] = 0.0;
666 MATRIX(resid_b)[i][m] = 0.0;
667 }
668 matrix_prod(resid_f, resid_f, 0, 1, ss_ff);
669 matrix_prod(resid_b, resid_b, 0, 1, ss_bb);
670 matrix_prod(resid_f, resid_b, 0, 1, ss_fb);
671
672 burg2(ss_ff, ss_bb, ss_fb, E, KA, KB); /* Update K */
673
674 for (i = 0; i <= m + 1; i++) {
675
676 matrix_prod(KA, subarray(B[m], m + 1 - i), 0, 0, tmp);
677 array_op(subarray(A[m], i), tmp, '-', subarray(A[m+1], i));
678
679 matrix_prod(KB, subarray(A[m], m + 1 - i), 0, 0, tmp);
680 array_op(subarray(B[m], i), tmp, '-', subarray(B[m+1], i));
681
682 }
683
684 matrix_prod(KA, resid_b, 0, 0, resid_f_tmp);
685 matrix_prod(KB, resid_f, 0, 0, resid_b_tmp);
686 array_op(resid_f, resid_f_tmp, '-', resid_f);
687 array_op(resid_b, resid_b_tmp, '-', resid_b);
688
689 if (vmethod == 1) {
690 matrix_prod(KA, KB, 0, 0, tmp);
691 array_op(id, tmp, '-', tmp);
692 matrix_prod(tmp, E, 0, 0, E);
693 }
694 else if (vmethod == 2) {
695 matrix_prod(resid_f, resid_f, 0, 1, E);
696 matrix_prod(resid_b, resid_b, 0, 1, tmp);
697 array_op(E, tmp, '+', E);
698 scalar_op(E, 2.0*(n - m - 1), '/', E);
699 }
700 else error(_("Invalid vmethod"));
701
702 copy_array(E, subarray(V,m+1));
703 copy_array(KA, subarray(P,m+1));
704 }
705 }
706
707
burg2(Array ss_ff,Array ss_bb,Array ss_fb,Array E,Array KA,Array KB)708 static void burg2(Array ss_ff, Array ss_bb, Array ss_fb, Array E,
709 Array KA, Array KB)
710 /*
711 Estimate partial correlation by minimizing (1/2)*log(det(s)) where
712 "s" is the the sum of the forward and backward prediction errors.
713
714 In the multivariate case, the forward (KA) and backward (KB) partial
715 correlation coefficients are related by
716
717 KA = solve(E) %*% t(KB) %*% E
718
719 where E is the prediction variance.
720
721 */
722 {
723 int i, j, k, l, nser = NROW(ss_ff);
724 int iter;
725 Array ss_bf;
726 Array s, tmp, d1;
727 Array D1, D2, THETA, THETAOLD, THETADIFF, TMP;
728 Array obj;
729 Array e, f, g, h, sg, sh;
730 Array theta;
731
732 ss_bf = make_zero_matrix(nser,nser);
733 transpose_matrix(ss_fb, ss_bf);
734 s = make_zero_matrix(nser, nser);
735 tmp = make_zero_matrix(nser, nser);
736 d1 = make_zero_matrix(nser, nser);
737
738 e = make_zero_matrix(nser, nser);
739 f = make_zero_matrix(nser, nser);
740 g = make_zero_matrix(nser, nser);
741 h = make_zero_matrix(nser, nser);
742 sg = make_zero_matrix(nser, nser);
743 sh = make_zero_matrix(nser, nser);
744
745 theta = make_zero_matrix(nser, nser);
746
747 D1 = make_zero_matrix(nser*nser, 1);
748 D2 = make_zero_matrix(nser*nser, nser*nser);
749 THETA = make_zero_matrix(nser*nser, 1); /* theta in vector form */
750 THETAOLD = make_zero_matrix(nser*nser, 1);
751 THETADIFF = make_zero_matrix(nser*nser, 1);
752 TMP = make_zero_matrix(nser*nser, 1);
753
754 obj = make_zero_matrix(1,1);
755
756 /* utility matrices e,f,g,h */
757 qr_solve(E, ss_bf, e);
758 qr_solve(E, ss_fb, f);
759 qr_solve(E, ss_bb, tmp);
760 transpose_matrix(tmp, tmp);
761 qr_solve(E, tmp, g);
762 qr_solve(E, ss_ff, tmp);
763 transpose_matrix(tmp, tmp);
764 qr_solve(E, tmp, h);
765
766 for(iter = 0; iter < BURG_MAX_ITER; iter++)
767 {
768 /* Forward and backward partial correlation coefficients */
769 transpose_matrix(theta, tmp);
770 qr_solve(E, tmp, tmp);
771 transpose_matrix(tmp, KA);
772
773 qr_solve(E, theta, tmp);
774 transpose_matrix(tmp, KB);
775
776 /* Sum of forward and backward prediction errors ... */
777 set_array_to_zero(s);
778
779 /* Forward */
780 array_op(s, ss_ff, '+', s);
781 matrix_prod(KA, ss_bf, 0, 0, tmp);
782 array_op(s, tmp, '-', s);
783 transpose_matrix(tmp, tmp);
784 array_op(s, tmp, '-', s);
785 matrix_prod(ss_bb, KA, 0, 1, tmp);
786 matrix_prod(KA, tmp, 0, 0, tmp);
787 array_op(s, tmp, '+', s);
788
789 /* Backward */
790 array_op(s, ss_bb, '+', s);
791 matrix_prod(KB, ss_fb, 0, 0, tmp);
792 array_op(s, tmp, '-', s);
793 transpose_matrix(tmp, tmp);
794 array_op(s, tmp, '-', s);
795 matrix_prod(ss_ff, KB, 0, 1, tmp);
796 matrix_prod(KB, tmp, 0, 0, tmp);
797 array_op(s, tmp, '+', s);
798
799 matrix_prod(s, f, 0, 0, d1);
800 matrix_prod(e, s, 1, 0, tmp);
801 array_op(d1, tmp, '+', d1);
802
803 /*matrix_prod(g,s,0,0,sg);*/
804 matrix_prod(s,g,0,0,sg);
805 matrix_prod(s,h,0,0,sh);
806
807 for (i = 0; i < nser; i++) {
808 for (j = 0; j < nser; j++) {
809 MATRIX(D1)[nser*i+j][0] = MATRIX(d1)[i][j];
810 for (k = 0; k < nser; k++)
811 for (l = 0; l < nser; l++) {
812 MATRIX(D2)[nser*i+j][nser*k+l] =
813 (i == k) * MATRIX(sg)[j][l] +
814 MATRIX(sh)[i][k] * (j == l);
815 }
816 }
817 }
818
819 copy_array(THETA, THETAOLD);
820 qr_solve(D2, D1, THETA);
821
822 for (i = 0; i < vector_length(theta); i++)
823 VECTOR(theta)[i] = VECTOR(THETA)[i];
824
825 matrix_prod(D2, THETA, 0, 0, TMP);
826
827 array_op(THETAOLD, THETA, '-', THETADIFF);
828 matrix_prod(D2, THETADIFF, 0, 0, TMP);
829 matrix_prod(THETADIFF, TMP, 1, 0, obj);
830 if (VECTOR(obj)[0] < BURG_TOL)
831 break;
832
833 }
834
835 if (iter == BURG_MAX_ITER)
836 error(_("Burg's algorithm failed to find partial correlation"));
837 }
838
839 /* Whittle's algorithm for autoregression estimation
840
841 multi_yw is the interface to R. It also handles model selection using AIC
842
843 whittle,whittle2 implement Whittle's recursion for solving the multivariate
844 Yule-Walker equations.
845
846 Notation
847
848 resid residuals (forward and backward)
849 A Estimates of forward autocorrelation coefficients
850 B Estimates of backward autocorrelation coefficients
851 EA,EB Prediction Variance
852 KA,KB Partial correlation coefficient
853 */
854
855 void multi_yw(double *acf, int *pn, int *pomax, int *pnser, double *coef,
856 double *pacf, double *var, double *aic, int *porder,
857 int *puseaic);
858 static void whittle(Array acf, int nlag, Array *A, Array *B, Array p_forward,
859 Array v_forward, Array p_back, Array v_back);
860 static void whittle2 (Array acf, Array Aold, Array Bold, int lag,
861 char *direction, Array A, Array K, Array E);
862
863
multi_yw(double * acf,int * pn,int * pomax,int * pnser,double * coef,double * pacf,double * var,double * aic,int * porder,int * useaic)864 void multi_yw(double *acf, int *pn, int *pomax, int *pnser, double *coef,
865 double *pacf, double *var, double *aic, int *porder, int *useaic)
866 {
867 int i, m;
868 int omax = *pomax, n = *pn, nser=*pnser, order=*porder;
869 double aicmin;
870 Array acf_array, p_forward, p_back, v_forward, v_back;
871 Array *A, *B;
872 int dim[3];
873
874 dim[0] = omax+1; dim[1] = dim[2] = nser;
875 acf_array = make_array(acf, dim, 3);
876 p_forward = make_array(pacf, dim, 3);
877 v_forward = make_array(var, dim, 3);
878
879 /* Backward equations (discarded) */
880 p_back= make_zero_array(dim, 3);
881 v_back= make_zero_array(dim, 3);
882
883 A = (Array *) R_alloc(omax+2, sizeof(Array));
884 B = (Array *) R_alloc(omax+2, sizeof(Array));
885 for (i = 0; i <= omax; i++) {
886 A[i] = make_zero_array(dim, 3);
887 B[i] = make_zero_array(dim, 3);
888 }
889 whittle(acf_array, omax, A, B, p_forward, v_forward, p_back, v_back);
890
891 /* Model order selection */
892
893 for (m = 0; m <= omax; m++) {
894 aic[m] = n * ldet(subarray(v_forward,m)) + 2 * m * nser * nser;
895 }
896 if (*useaic) {
897 order = 0;
898 aicmin = aic[0];
899 for (m = 0; m <= omax; m++) {
900 if (aic[m] < aicmin) {
901 aicmin = aic[m];
902 order = m;
903 }
904 }
905 }
906 else order = omax;
907 *porder = order;
908
909 for(i = 0; i < vector_length(A[order]); i++)
910 coef[i] = VECTOR(A[order])[i];
911 }
912
whittle(Array acf,int nlag,Array * A,Array * B,Array p_forward,Array v_forward,Array p_back,Array v_back)913 static void whittle(Array acf, int nlag, Array *A, Array *B, Array p_forward,
914 Array v_forward, Array p_back, Array v_back)
915 {
916
917 int lag, nser = DIM(acf)[1];
918 const void *vmax;
919 Array EA, EB; /* prediction variance */
920 Array KA, KB; /* partial correlation coefficient */
921 Array id, tmp;
922
923 vmax = vmaxget();
924
925 KA = make_zero_matrix(nser, nser);
926 EA = make_zero_matrix(nser, nser);
927
928 KB = make_zero_matrix(nser, nser);
929 EB = make_zero_matrix(nser, nser);
930
931 id = make_identity_matrix(nser);
932
933 copy_array(id, subarray(A[0],0));
934 copy_array(id, subarray(B[0],0));
935 copy_array(id, subarray(p_forward,0));
936 copy_array(id, subarray(p_back,0));
937
938 for (lag = 1; lag <= nlag; lag++) {
939
940 whittle2(acf, A[lag-1], B[lag-1], lag, "forward", A[lag], KA, EB);
941 whittle2(acf, B[lag-1], A[lag-1], lag, "back", B[lag], KB, EA);
942
943 copy_array(EA, subarray(v_forward,lag-1));
944 copy_array(EB, subarray(v_back,lag-1));
945
946 copy_array(KA, subarray(p_forward,lag));
947 copy_array(KB, subarray(p_back,lag));
948
949 }
950
951 tmp = make_zero_matrix(nser,nser);
952
953 matrix_prod(KB,KA, 1, 1, tmp);
954 array_op(id, tmp, '-', tmp);
955 matrix_prod(EA, tmp, 0, 0, subarray(v_forward, nlag));
956
957 vmaxset(vmax);
958
959 }
960
whittle2(Array acf,Array Aold,Array Bold,int lag,char * direction,Array A,Array K,Array E)961 static void whittle2 (Array acf, Array Aold, Array Bold, int lag,
962 char *direction, Array A, Array K, Array E)
963 {
964
965 int d, i, nser=DIM(acf)[1];
966 const void *vmax;
967 Array beta, tmp, id;
968
969 d = strcmp(direction, "forward") == 0;
970
971 vmax = vmaxget();
972
973 beta = make_zero_matrix(nser,nser);
974 tmp = make_zero_matrix(nser, nser);
975 id = make_identity_matrix(nser);
976
977 set_array_to_zero(E);
978 copy_array(id, subarray(A,0));
979
980 for(i = 0; i < lag; i++) {
981 matrix_prod(subarray(acf,lag - i), subarray(Aold,i), d, 1, tmp);
982 array_op(beta, tmp, '+', beta);
983 matrix_prod(subarray(acf,i), subarray(Bold,i), d, 1, tmp);
984 array_op(E, tmp, '+', E);
985 }
986 qr_solve(E, beta, K);
987 transpose_matrix(K,K);
988 for (i = 1; i <= lag; i++) {
989 matrix_prod(K, subarray(Bold,lag - i), 0, 0, tmp);
990 array_op(subarray(Aold,i), tmp, '-', subarray(A,i));
991 }
992
993 vmaxset(vmax);
994 }
995