| /* |
| * Copyright (C) 1999 Martyn Plummer |
| * Copyright (C) 1999-2016 The R Core Team |
| * |
| * This program is free software; you can redistribute it and/or modify |
| * it under the terms of the GNU General Public License as published by |
| * the Free Software Foundation; either version 2 of the License, or |
| * (at your option) any later version. |
| * |
| * This program is distributed in the hope that it will be useful, |
| * but WITHOUT ANY WARRANTY; without even the implied warranty of |
| * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| * GNU General Public License for more details. |
| * |
| * You should have received a copy of the GNU General Public License |
| * along with this program; if not, a copy is available at |
| * https://www.R-project.org/Licenses/. |
| */ |
| |
| #include <math.h> |
| #include <string.h> |
| #include <R.h> |
| #include <R_ext/Applic.h> /* Fortran routines */ |
| #include "ts.h" |
| #include "stats.h" |
| |
| |
| #define MAX_DIM_LENGTH 4 |
| |
| #define VECTOR(x) (x.vec) |
| #define MATRIX(x) (x.mat) |
| #define ARRAY1(x) (x.vec) |
| #define ARRAY2(x) (x.mat) |
| #define ARRAY3(x) (x.arr3) |
| #define ARRAY4(x) (x.arr4) |
| #define DIM(x) (x.dim) |
| #define NROW(x) (x.dim[0]) |
| #define NCOL(x) (x.dim[1]) |
| #define DIM_LENGTH(x) (x.ndim) |
| |
| |
| typedef struct array { |
| double *vec; |
| double **mat; |
| double ***arr3; |
| double ****arr4; |
| int dim[MAX_DIM_LENGTH]; |
| int ndim; |
| } Array; |
| |
| static Array make_array(double vec[], int dim[], int ndim); |
| static Array make_zero_array(int dim[], int ndim); |
| static Array make_matrix(double vec[], int nrow, int ncol); |
| static Array make_zero_matrix(int nrow, int ncol); |
| static Array make_identity_matrix(int n); |
| |
| static Array subarray(Array a, int index); |
| |
| static int vector_length(Array a); |
| |
| static void set_array_to_zero(Array arr); |
| static void copy_array (Array orig, Array ans); |
| static void array_op(Array arr1, Array arr2, char op, Array ans); |
| static void scalar_op(Array arr, double s, char op, Array ans); |
| |
| static void transpose_matrix(Array mat, Array ans); |
| static void matrix_prod(Array mat1, Array mat2, int trans1, int trans2, |
| Array ans); |
| |
| |
| /* Functions for dynamically allocating arrays |
| |
| The Array structure contains pointers to arrays which are allocated |
| using the R_alloc function. Although the .C() interface cleans up |
| all memory assigned with R_alloc, judicious use of vmaxget() vmaxset() |
| to free this memory is probably wise. See memory.c in R core. |
| |
| */ |
| |
| static void assert(int bool) |
| { |
| if(!bool) |
| error(("assert failed in src/library/ts/src/carray.c")); |
| } |
| |
| static Array init_array(void) |
| { |
| int i; |
| Array a; |
| |
| /* Initialize everything to zero. Useful for debugging */ |
| ARRAY1(a) = (double *) '\0'; |
| ARRAY2(a) = (double **) '\0'; |
| ARRAY3(a) = (double ***) '\0'; |
| ARRAY4(a) = (double ****) '\0'; |
| for (i = 0; i < MAX_DIM_LENGTH; i++) |
| DIM(a)[i] = 0; |
| DIM_LENGTH(a) = 0; |
| |
| return a; |
| } |
| |
| static int vector_length(Array a) |
| { |
| int i, len; |
| |
| for (i = 0, len = 1; i < DIM_LENGTH(a); i++) { |
| len *= DIM(a)[i]; |
| } |
| |
| return len; |
| } |
| |
| |
| static Array make_array(double vec[], int dim[], int ndim) |
| { |
| int d, i, j; |
| int len[MAX_DIM_LENGTH + 1]; |
| Array a; |
| |
| assert(ndim <= MAX_DIM_LENGTH); |
| |
| a = init_array(); |
| |
| len[ndim] = 1; |
| for (d = ndim; d >= 1; d--) { |
| len[d-1] = len[d] * dim[ndim - d]; |
| } |
| |
| for (d = 1; d <= ndim; d++) { |
| switch(d) { |
| case 1: |
| VECTOR(a) = vec; |
| break; |
| case 2: |
| ARRAY2(a) = (double**) R_alloc(len[2 - 1],sizeof(double*)); |
| for(i = 0, j = 0; i < len[2 - 1]; i++, j+=dim[ndim - 2 + 1]) { |
| ARRAY2(a)[i] = ARRAY1(a) + j; |
| } |
| break; |
| case 3: |
| ARRAY3(a) = (double***) R_alloc(len[3 - 1],sizeof(double**)); |
| for(i = 0, j = 0; i < len[3 - 1]; i++, j+=dim[ndim - 3 + 1]) { |
| ARRAY3(a)[i] = ARRAY2(a) + j; |
| } |
| break; |
| case 4: |
| ARRAY4(a) = (double****) R_alloc(len[4 - 1],sizeof(double***)); |
| for(i = 0, j = 0; i < len[4 - 1]; i++, j+=dim[ndim - 4 + 1]) { |
| ARRAY4(a)[i] = ARRAY3(a) + j; |
| } |
| break; |
| default: |
| break; |
| } |
| } |
| |
| for (i = 0; i < ndim; i++) { |
| DIM(a)[i] = dim[i]; |
| } |
| DIM_LENGTH(a) = ndim; |
| |
| return a; |
| } |
| |
| static Array make_zero_array(int dim[], int ndim) |
| { |
| int i; |
| int len; |
| double *vec; |
| |
| for (i = 0, len = 1; i < ndim; i++) { |
| len *= dim[i]; |
| } |
| |
| vec = (double *) R_alloc(len, sizeof(double)); |
| for (i = 0; i < len; i++) { |
| vec[i] = 0.0; |
| } |
| |
| return make_array(vec, dim, ndim); |
| |
| } |
| |
| static Array make_matrix(double vec[], int nrow, int ncol) |
| { |
| int dim[2]; |
| |
| dim[0] = nrow; |
| dim[1] = ncol; |
| return make_array(vec, dim, 2); |
| } |
| |
| static Array make_zero_matrix(int nrow, int ncol) |
| { |
| int dim[2]; |
| Array a; |
| |
| dim[0] = nrow; |
| dim[1] = ncol; |
| a = make_zero_array(dim, 2); |
| return a; |
| } |
| |
| static Array subarray(Array a, int index) |
| /* Return subarray of array a in the form of an Array |
| structure so it can be manipulated by other functions |
| NB The data are not copied, so any changes made to the |
| subarray will affect the original array. |
| */ |
| { |
| int i, offset; |
| Array b; |
| |
| b = init_array(); |
| |
| /* is index in range? */ |
| assert( index >= 0 && index < DIM(a)[0] ); |
| |
| offset = index; |
| switch(DIM_LENGTH(a)) { |
| /* NB Falling through here */ |
| case 4: |
| offset *= DIM(a)[DIM_LENGTH(a) - 4 + 1]; |
| ARRAY3(b) = ARRAY3(a) + offset; |
| case 3: |
| offset *= DIM(a)[DIM_LENGTH(a) - 3 + 1]; |
| ARRAY2(b) = ARRAY2(a) + offset; |
| case 2: |
| offset *= DIM(a)[DIM_LENGTH(a) - 2 + 1]; |
| ARRAY1(b) = ARRAY1(a) + offset; |
| break; |
| default: |
| break; |
| } |
| |
| |
| DIM_LENGTH(b) = DIM_LENGTH(a) - 1; |
| |
| for (i = 0; i < DIM_LENGTH(b); i++) |
| DIM(b)[i] = DIM(a)[i+1]; |
| |
| return b; |
| |
| } |
| |
| static int test_array_conform(Array a1, Array a2) |
| { |
| int i, ans = FALSE; |
| |
| if (DIM_LENGTH(a1) != DIM_LENGTH(a2)) |
| return FALSE; |
| else |
| for (i = 0; i < DIM_LENGTH(a1); i++) { |
| if (DIM(a1)[i] == DIM(a2)[i]) |
| ans = TRUE; |
| else |
| return FALSE; |
| } |
| |
| return ans; |
| } |
| |
| static void copy_array (Array orig, Array ans) |
| /* copy matrix orig to ans */ |
| { |
| int i; |
| |
| assert (test_array_conform(orig, ans)); |
| |
| for(i = 0; i < vector_length(orig); i++) |
| VECTOR(ans)[i] = VECTOR(orig)[i]; |
| } |
| |
| static void transpose_matrix(Array mat, Array ans) |
| { |
| int i,j; |
| const void *vmax; |
| Array tmp; |
| |
| tmp = init_array(); |
| |
| assert(DIM_LENGTH(mat) == 2 && DIM_LENGTH(ans) == 2); |
| assert(NCOL(mat) == NROW(ans)); |
| assert(NROW(mat) == NCOL(ans)); |
| |
| vmax = vmaxget(); |
| |
| tmp = make_zero_matrix(NROW(ans), NCOL(ans)); |
| for(i = 0; i < NROW(mat); i++) |
| for(j = 0; j < NCOL(mat); j++) |
| MATRIX(tmp)[j][i] = MATRIX(mat)[i][j]; |
| copy_array(tmp, ans); |
| |
| vmaxset(vmax); |
| } |
| |
| static void array_op(Array arr1, Array arr2, char op, Array ans) |
| /* Element-wise array operations */ |
| { |
| int i; |
| |
| assert (test_array_conform(arr1, arr2)); |
| assert (test_array_conform(arr2, ans)); |
| |
| switch (op) { |
| case '*': |
| for (i = 0; i < vector_length(ans); i++) |
| VECTOR(ans)[i] = VECTOR(arr1)[i] * VECTOR(arr2)[i]; |
| break; |
| case '+': |
| for (i = 0; i < vector_length(ans); i++) |
| VECTOR(ans)[i] = VECTOR(arr1)[i] + VECTOR(arr2)[i]; |
| break; |
| case '/': |
| for (i = 0; i < vector_length(ans); i++) |
| VECTOR(ans)[i] = VECTOR(arr1)[i] / VECTOR(arr2)[i]; |
| break; |
| case '-': |
| for (i = 0; i < vector_length(ans); i++) |
| VECTOR(ans)[i] = VECTOR(arr1)[i] - VECTOR(arr2)[i]; |
| break; |
| default: |
| printf("Unknown op in array_op"); |
| } |
| } |
| |
| |
| static void scalar_op(Array arr, double s, char op, Array ans) |
| /* Elementwise scalar operations */ |
| { |
| int i; |
| |
| assert (test_array_conform(arr, ans)); |
| |
| switch (op) { |
| case '*': |
| for (i = 0; i < vector_length(ans); i++) |
| VECTOR(ans)[i] = VECTOR(arr)[i] * s; |
| break; |
| case '+': |
| for (i = 0; i < vector_length(ans); i++) |
| VECTOR(ans)[i] = VECTOR(arr)[i] + s; |
| break; |
| case '/': |
| for (i = 0; i < vector_length(ans); i++) |
| VECTOR(ans)[i] = VECTOR(arr)[i] / s; |
| break; |
| case '-': |
| for (i = 0; i < vector_length(ans); i++) |
| VECTOR(ans)[i] = VECTOR(arr)[i] - s; |
| break; |
| default: |
| printf("Unknown op in array_op"); |
| } |
| } |
| |
| static void matrix_prod(Array mat1, Array mat2, int trans1, int trans2, Array ans) |
| /* |
| General matrix product between mat1 and mat2. Put answer in ans. |
| trans1 and trans2 are logical flags which indicate if the matrix is |
| to be transposed. Normal matrix multiplication has trans1 = trans2 = 0. |
| */ |
| { |
| int i,j,k,K1,K2; |
| const void *vmax; |
| double m1, m2; |
| Array tmp; |
| |
| /* Test whether everything is a matrix */ |
| assert(DIM_LENGTH(mat1) == 2 && |
| DIM_LENGTH(mat2) == 2 && DIM_LENGTH(ans) == 2); |
| |
| /* Test whether matrices conform. K is the dimension that is |
| lost by multiplication */ |
| if (trans1) { |
| assert ( NCOL(mat1) == NROW(ans) ); |
| K1 = NROW(mat1); |
| } |
| else { |
| assert ( NROW(mat1) == NROW(ans) ); |
| K1 = NCOL(mat1); |
| } |
| if (trans2) { |
| assert ( NROW(mat2) == NCOL(ans) ); |
| K2 = NCOL(mat2); |
| } |
| else { |
| assert ( NCOL(mat2) == NCOL(ans) ); |
| K2 = NROW(mat2); |
| } |
| assert (K1 == K2); |
| |
| tmp = init_array(); |
| |
| /* In case ans is the same as mat1 or mat2, we create a temporary |
| matrix to hold the answer, then copy it to ans |
| */ |
| vmax = vmaxget(); |
| |
| tmp = make_zero_matrix(NROW(ans), NCOL(ans)); |
| for (i = 0; i < NROW(tmp); i++) { |
| for (j = 0; j < NCOL(tmp); j++) { |
| for(k = 0; k < K1; k++) { |
| m1 = (trans1) ? MATRIX(mat1)[k][i] : MATRIX(mat1)[i][k]; |
| m2 = (trans2) ? MATRIX(mat2)[j][k] : MATRIX(mat2)[k][j]; |
| MATRIX(tmp)[i][j] += m1 * m2; |
| } |
| } |
| } |
| copy_array(tmp, ans); |
| |
| vmaxset(vmax); |
| } |
| |
| static void set_array_to_zero(Array arr) |
| { |
| int i; |
| |
| for (i = 0; i < vector_length(arr); i++) |
| VECTOR(arr)[i] = 0.0; |
| } |
| |
| static Array make_identity_matrix(int n) |
| { |
| int i; |
| Array a; |
| |
| a = make_zero_matrix(n,n); |
| for(i = 0; i < n; i++) |
| MATRIX(a)[i][i] = 1.0; |
| |
| return a; |
| } |
| |
| static void qr_solve(Array x, Array y, Array coef) |
| /* Translation of the R function qr.solve into pure C |
| NB We have to transpose the matrices since the ordering of an array is different in Fortran |
| NB2 We have to copy x to avoid it being overwritten. |
| */ |
| { |
| int i, info = 0, rank, *pivot, n, p; |
| const void *vmax; |
| double tol = 1.0E-7, *qraux, *work; |
| Array xt, yt, coeft; |
| |
| assert(NROW(x) == NROW(y)); |
| assert(NCOL(coef) == NCOL(y)); |
| assert(NCOL(x) == NROW(coef)); |
| |
| vmax = vmaxget(); |
| |
| qraux = (double *) R_alloc(NCOL(x), sizeof(double)); |
| pivot = (int *) R_alloc(NCOL(x), sizeof(int)); |
| work = (double *) R_alloc(2*NCOL(x), sizeof(double)); |
| |
| for(i = 0; i < NCOL(x); i++) |
| pivot[i] = i+1; |
| |
| xt = make_zero_matrix(NCOL(x), NROW(x)); |
| transpose_matrix(x,xt); |
| |
| n = NROW(x); |
| p = NCOL(x); |
| |
| F77_CALL(dqrdc2)(VECTOR(xt), &n, &n, &p, &tol, &rank, |
| qraux, pivot, work); |
| |
| if (rank != p) |
| error(_("Singular matrix in qr_solve")); |
| |
| yt = make_zero_matrix(NCOL(y), NROW(y)); |
| coeft = make_zero_matrix(NCOL(coef), NROW(coef)); |
| transpose_matrix(y, yt); |
| |
| F77_CALL(dqrcf)(VECTOR(xt), &NROW(x), &rank, qraux, |
| yt.vec, &NCOL(y), coeft.vec, &info); |
| |
| transpose_matrix(coeft,coef); |
| |
| vmaxset(vmax); |
| } |
| |
| static double ldet(Array x) |
| /* Log determinant of square matrix */ |
| { |
| int i, rank, *pivot, n, p; |
| const void *vmax; |
| double ll, tol = 1.0E-7, *qraux, *work; |
| Array xtmp; |
| |
| assert(DIM_LENGTH(x) == 2); /* is x a matrix? */ |
| assert(NROW(x) == NCOL(x)); /* is x square? */ |
| |
| vmax = vmaxget(); |
| |
| qraux = (double *) R_alloc(NCOL(x), sizeof(double)); |
| pivot = (int *) R_alloc(NCOL(x), sizeof(int)); |
| work = (double *) R_alloc(2*NCOL(x), sizeof(double)); |
| |
| xtmp = make_zero_matrix(NROW(x), NCOL(x)); |
| copy_array(x, xtmp); |
| |
| for(i = 0; i < NCOL(x); i++) |
| pivot[i] = i+1; |
| |
| p = n = NROW(x); |
| |
| F77_CALL(dqrdc2)(VECTOR(xtmp), &n, &n, &p, &tol, &rank, |
| qraux, pivot, work); |
| |
| if (rank != p) |
| error(_("Singular matrix in ldet")); |
| |
| for (i = 0, ll=0.0; i < rank; i++) { |
| ll += log(fabs(MATRIX(xtmp)[i][i])); |
| } |
| |
| vmaxset(vmax); |
| |
| return ll; |
| } |
| |
| |
| |
| /* Burg's algorithm for autoregression estimation |
| |
| multi_burg is the interface to R. It also handles model selection |
| using AIC |
| |
| burg implements the main part of the algorithm |
| |
| burg2 estimates the partial correlation coefficient. This |
| requires iteration in the multivariate case. |
| |
| Notation |
| |
| resid residuals (forward and backward) |
| A Estimates of autocorrelation coefficients |
| V Prediction Variance |
| K Partial correlation coefficient |
| */ |
| |
| |
| #define BURG_MAX_ITER 20 |
| #define BURG_TOL 1.0E-8 |
| |
| void multi_burg(int *pn, double *x, int *pomax, int *pnser, double *coef, |
| double *pacf, double *var, double *aic, int *porder, int *useaic, |
| int *vmethod); |
| static void burg0(int omax, Array resid_f, Array resid_b, Array *A, Array *B, |
| Array P, Array V, int vmethod); |
| static void burg2(Array ss_ff, Array ss_bb, Array ss_fb, Array E, |
| Array KA, Array KB); |
| |
| void multi_burg(int *pn, double *x, int *pomax, int *pnser, double *coef, |
| double *pacf, double *var, double *aic, int *porder, int *useaic, |
| int *vmethod) |
| { |
| int i, j, m, omax = *pomax, n = *pn, nser=*pnser, order=*porder; |
| int dim1[3]; |
| double aicmin; |
| Array xarr, resid_f, resid_b, resid_f_tmp; |
| Array *A, *B, P, V; |
| |
| dim1[0] = omax+1; dim1[1] = dim1[2] = nser; |
| A = (Array *) R_alloc(omax+1, sizeof(Array)); |
| B = (Array *) R_alloc(omax+1, sizeof(Array)); |
| for (i = 0; i <= omax; i++) { |
| A[i] = make_zero_array(dim1, 3); |
| B[i] = make_zero_array(dim1, 3); |
| } |
| P = make_array(pacf, dim1, 3); |
| V = make_array(var, dim1, 3); |
| |
| xarr = make_matrix(x, nser, n); |
| resid_f = make_zero_matrix(nser, n); |
| resid_b = make_zero_matrix(nser, n); |
| set_array_to_zero(resid_b); |
| copy_array(xarr, resid_f); |
| copy_array(xarr, resid_b); |
| resid_f_tmp = make_zero_matrix(nser, n); |
| |
| burg0(omax, resid_f, resid_b, A, B, P, V, *vmethod); |
| |
| /* Model order selection */ |
| |
| for (i = 0; i <= omax; i++) { |
| aic[i] = n * ldet(subarray(V,i)) + 2 * i * nser * nser; |
| } |
| if (*useaic) { |
| order = 0; |
| aicmin = aic[0]; |
| for (i = 1; i <= omax; i++) { |
| if (aic[i] < aicmin) { |
| aicmin = aic[i]; |
| order = i; |
| } |
| } |
| } |
| else order = omax; |
| *porder = order; |
| |
| for(i = 0; i < vector_length(A[order]); i++) |
| coef[i] = VECTOR(A[order])[i]; |
| |
| if (*useaic) { |
| /* Recalculate residuals for chosen model */ |
| set_array_to_zero(resid_f); |
| set_array_to_zero(resid_f_tmp); |
| for (m = 0; m <= order; m++) { |
| for (i = 0; i < NROW(resid_f_tmp); i++) { |
| for (j = 0; j < NCOL(resid_f_tmp) - order; j++) { |
| MATRIX(resid_f_tmp)[i][j + order] = MATRIX(xarr)[i][j + order - m]; |
| } |
| } |
| matrix_prod(subarray(A[order],m), resid_f_tmp, 0, 0, resid_f_tmp); |
| array_op(resid_f_tmp, resid_f, '+', resid_f); |
| } |
| } |
| copy_array(resid_f, xarr); |
| |
| } |
| |
| |
| static void burg0(int omax, Array resid_f, Array resid_b, Array *A, Array *B, |
| Array P, Array V, int vmethod) |
| { |
| int i, j, m, n = NCOL(resid_f), nser=NROW(resid_f); |
| Array ss_ff, ss_bb, ss_fb; |
| Array resid_f_tmp, resid_b_tmp; |
| Array KA, KB, E; |
| Array id, tmp; |
| |
| ss_ff = make_zero_matrix(nser, nser); |
| ss_fb = make_zero_matrix(nser, nser); |
| ss_bb = make_zero_matrix(nser, nser); |
| |
| resid_f_tmp = make_zero_matrix(nser, n); |
| resid_b_tmp = make_zero_matrix(nser, n); |
| |
| id = make_identity_matrix(nser); |
| |
| tmp = make_zero_matrix(nser, nser); |
| |
| E = make_zero_matrix(nser, nser); |
| KA = make_zero_matrix(nser, nser); |
| KB = make_zero_matrix(nser, nser); |
| |
| set_array_to_zero(A[0]); |
| set_array_to_zero(B[0]); |
| copy_array(id, subarray(A[0],0)); |
| copy_array(id, subarray(B[0],0)); |
| |
| matrix_prod(resid_f, resid_f, 0, 1, E); |
| scalar_op(E, n, '/', E); |
| copy_array(E, subarray(V,0)); |
| |
| for (m = 0; m < omax; m++) { |
| |
| for(i = 0; i < nser; i++) { |
| for (j = n - 1; j > m; j--) { |
| MATRIX(resid_b)[i][j] = MATRIX(resid_b)[i][j-1]; |
| } |
| MATRIX(resid_f)[i][m] = 0.0; |
| MATRIX(resid_b)[i][m] = 0.0; |
| } |
| matrix_prod(resid_f, resid_f, 0, 1, ss_ff); |
| matrix_prod(resid_b, resid_b, 0, 1, ss_bb); |
| matrix_prod(resid_f, resid_b, 0, 1, ss_fb); |
| |
| burg2(ss_ff, ss_bb, ss_fb, E, KA, KB); /* Update K */ |
| |
| for (i = 0; i <= m + 1; i++) { |
| |
| matrix_prod(KA, subarray(B[m], m + 1 - i), 0, 0, tmp); |
| array_op(subarray(A[m], i), tmp, '-', subarray(A[m+1], i)); |
| |
| matrix_prod(KB, subarray(A[m], m + 1 - i), 0, 0, tmp); |
| array_op(subarray(B[m], i), tmp, '-', subarray(B[m+1], i)); |
| |
| } |
| |
| matrix_prod(KA, resid_b, 0, 0, resid_f_tmp); |
| matrix_prod(KB, resid_f, 0, 0, resid_b_tmp); |
| array_op(resid_f, resid_f_tmp, '-', resid_f); |
| array_op(resid_b, resid_b_tmp, '-', resid_b); |
| |
| if (vmethod == 1) { |
| matrix_prod(KA, KB, 0, 0, tmp); |
| array_op(id, tmp, '-', tmp); |
| matrix_prod(tmp, E, 0, 0, E); |
| } |
| else if (vmethod == 2) { |
| matrix_prod(resid_f, resid_f, 0, 1, E); |
| matrix_prod(resid_b, resid_b, 0, 1, tmp); |
| array_op(E, tmp, '+', E); |
| scalar_op(E, 2.0*(n - m - 1), '/', E); |
| } |
| else error(_("Invalid vmethod")); |
| |
| copy_array(E, subarray(V,m+1)); |
| copy_array(KA, subarray(P,m+1)); |
| } |
| } |
| |
| |
| static void burg2(Array ss_ff, Array ss_bb, Array ss_fb, Array E, |
| Array KA, Array KB) |
| /* |
| Estimate partial correlation by minimizing (1/2)*log(det(s)) where |
| "s" is the the sum of the forward and backward prediction errors. |
| |
| In the multivariate case, the forward (KA) and backward (KB) partial |
| correlation coefficients are related by |
| |
| KA = solve(E) %*% t(KB) %*% E |
| |
| where E is the prediction variance. |
| |
| */ |
| { |
| int i, j, k, l, nser = NROW(ss_ff); |
| int iter; |
| Array ss_bf; |
| Array s, tmp, d1; |
| Array D1, D2, THETA, THETAOLD, THETADIFF, TMP; |
| Array obj; |
| Array e, f, g, h, sg, sh; |
| Array theta; |
| |
| ss_bf = make_zero_matrix(nser,nser); |
| transpose_matrix(ss_fb, ss_bf); |
| s = make_zero_matrix(nser, nser); |
| tmp = make_zero_matrix(nser, nser); |
| d1 = make_zero_matrix(nser, nser); |
| |
| e = make_zero_matrix(nser, nser); |
| f = make_zero_matrix(nser, nser); |
| g = make_zero_matrix(nser, nser); |
| h = make_zero_matrix(nser, nser); |
| sg = make_zero_matrix(nser, nser); |
| sh = make_zero_matrix(nser, nser); |
| |
| theta = make_zero_matrix(nser, nser); |
| |
| D1 = make_zero_matrix(nser*nser, 1); |
| D2 = make_zero_matrix(nser*nser, nser*nser); |
| THETA = make_zero_matrix(nser*nser, 1); /* theta in vector form */ |
| THETAOLD = make_zero_matrix(nser*nser, 1); |
| THETADIFF = make_zero_matrix(nser*nser, 1); |
| TMP = make_zero_matrix(nser*nser, 1); |
| |
| obj = make_zero_matrix(1,1); |
| |
| /* utility matrices e,f,g,h */ |
| qr_solve(E, ss_bf, e); |
| qr_solve(E, ss_fb, f); |
| qr_solve(E, ss_bb, tmp); |
| transpose_matrix(tmp, tmp); |
| qr_solve(E, tmp, g); |
| qr_solve(E, ss_ff, tmp); |
| transpose_matrix(tmp, tmp); |
| qr_solve(E, tmp, h); |
| |
| for(iter = 0; iter < BURG_MAX_ITER; iter++) |
| { |
| /* Forward and backward partial correlation coefficients */ |
| transpose_matrix(theta, tmp); |
| qr_solve(E, tmp, tmp); |
| transpose_matrix(tmp, KA); |
| |
| qr_solve(E, theta, tmp); |
| transpose_matrix(tmp, KB); |
| |
| /* Sum of forward and backward prediction errors ... */ |
| set_array_to_zero(s); |
| |
| /* Forward */ |
| array_op(s, ss_ff, '+', s); |
| matrix_prod(KA, ss_bf, 0, 0, tmp); |
| array_op(s, tmp, '-', s); |
| transpose_matrix(tmp, tmp); |
| array_op(s, tmp, '-', s); |
| matrix_prod(ss_bb, KA, 0, 1, tmp); |
| matrix_prod(KA, tmp, 0, 0, tmp); |
| array_op(s, tmp, '+', s); |
| |
| /* Backward */ |
| array_op(s, ss_bb, '+', s); |
| matrix_prod(KB, ss_fb, 0, 0, tmp); |
| array_op(s, tmp, '-', s); |
| transpose_matrix(tmp, tmp); |
| array_op(s, tmp, '-', s); |
| matrix_prod(ss_ff, KB, 0, 1, tmp); |
| matrix_prod(KB, tmp, 0, 0, tmp); |
| array_op(s, tmp, '+', s); |
| |
| matrix_prod(s, f, 0, 0, d1); |
| matrix_prod(e, s, 1, 0, tmp); |
| array_op(d1, tmp, '+', d1); |
| |
| /*matrix_prod(g,s,0,0,sg);*/ |
| matrix_prod(s,g,0,0,sg); |
| matrix_prod(s,h,0,0,sh); |
| |
| for (i = 0; i < nser; i++) { |
| for (j = 0; j < nser; j++) { |
| MATRIX(D1)[nser*i+j][0] = MATRIX(d1)[i][j]; |
| for (k = 0; k < nser; k++) |
| for (l = 0; l < nser; l++) { |
| MATRIX(D2)[nser*i+j][nser*k+l] = |
| (i == k) * MATRIX(sg)[j][l] + |
| MATRIX(sh)[i][k] * (j == l); |
| } |
| } |
| } |
| |
| copy_array(THETA, THETAOLD); |
| qr_solve(D2, D1, THETA); |
| |
| for (i = 0; i < vector_length(theta); i++) |
| VECTOR(theta)[i] = VECTOR(THETA)[i]; |
| |
| matrix_prod(D2, THETA, 0, 0, TMP); |
| |
| array_op(THETAOLD, THETA, '-', THETADIFF); |
| matrix_prod(D2, THETADIFF, 0, 0, TMP); |
| matrix_prod(THETADIFF, TMP, 1, 0, obj); |
| if (VECTOR(obj)[0] < BURG_TOL) |
| break; |
| |
| } |
| |
| if (iter == BURG_MAX_ITER) |
| error(_("Burg's algorithm failed to find partial correlation")); |
| } |
| |
| /* Whittle's algorithm for autoregression estimation |
| |
| multi_yw is the interface to R. It also handles model selection using AIC |
| |
| whittle,whittle2 implement Whittle's recursion for solving the multivariate |
| Yule-Walker equations. |
| |
| Notation |
| |
| resid residuals (forward and backward) |
| A Estimates of forward autocorrelation coefficients |
| B Estimates of backward autocorrelation coefficients |
| EA,EB Prediction Variance |
| KA,KB Partial correlation coefficient |
| */ |
| |
| void multi_yw(double *acf, int *pn, int *pomax, int *pnser, double *coef, |
| double *pacf, double *var, double *aic, int *porder, |
| int *puseaic); |
| static void whittle(Array acf, int nlag, Array *A, Array *B, Array p_forward, |
| Array v_forward, Array p_back, Array v_back); |
| static void whittle2 (Array acf, Array Aold, Array Bold, int lag, |
| char *direction, Array A, Array K, Array E); |
| |
| |
| void multi_yw(double *acf, int *pn, int *pomax, int *pnser, double *coef, |
| double *pacf, double *var, double *aic, int *porder, int *useaic) |
| { |
| int i, m; |
| int omax = *pomax, n = *pn, nser=*pnser, order=*porder; |
| double aicmin; |
| Array acf_array, p_forward, p_back, v_forward, v_back; |
| Array *A, *B; |
| int dim[3]; |
| |
| dim[0] = omax+1; dim[1] = dim[2] = nser; |
| acf_array = make_array(acf, dim, 3); |
| p_forward = make_array(pacf, dim, 3); |
| v_forward = make_array(var, dim, 3); |
| |
| /* Backward equations (discarded) */ |
| p_back= make_zero_array(dim, 3); |
| v_back= make_zero_array(dim, 3); |
| |
| A = (Array *) R_alloc(omax+2, sizeof(Array)); |
| B = (Array *) R_alloc(omax+2, sizeof(Array)); |
| for (i = 0; i <= omax; i++) { |
| A[i] = make_zero_array(dim, 3); |
| B[i] = make_zero_array(dim, 3); |
| } |
| whittle(acf_array, omax, A, B, p_forward, v_forward, p_back, v_back); |
| |
| /* Model order selection */ |
| |
| for (m = 0; m <= omax; m++) { |
| aic[m] = n * ldet(subarray(v_forward,m)) + 2 * m * nser * nser; |
| } |
| if (*useaic) { |
| order = 0; |
| aicmin = aic[0]; |
| for (m = 0; m <= omax; m++) { |
| if (aic[m] < aicmin) { |
| aicmin = aic[m]; |
| order = m; |
| } |
| } |
| } |
| else order = omax; |
| *porder = order; |
| |
| for(i = 0; i < vector_length(A[order]); i++) |
| coef[i] = VECTOR(A[order])[i]; |
| } |
| |
| static void whittle(Array acf, int nlag, Array *A, Array *B, Array p_forward, |
| Array v_forward, Array p_back, Array v_back) |
| { |
| |
| int lag, nser = DIM(acf)[1]; |
| const void *vmax; |
| Array EA, EB; /* prediction variance */ |
| Array KA, KB; /* partial correlation coefficient */ |
| Array id, tmp; |
| |
| vmax = vmaxget(); |
| |
| KA = make_zero_matrix(nser, nser); |
| EA = make_zero_matrix(nser, nser); |
| |
| KB = make_zero_matrix(nser, nser); |
| EB = make_zero_matrix(nser, nser); |
| |
| id = make_identity_matrix(nser); |
| |
| copy_array(id, subarray(A[0],0)); |
| copy_array(id, subarray(B[0],0)); |
| copy_array(id, subarray(p_forward,0)); |
| copy_array(id, subarray(p_back,0)); |
| |
| for (lag = 1; lag <= nlag; lag++) { |
| |
| whittle2(acf, A[lag-1], B[lag-1], lag, "forward", A[lag], KA, EB); |
| whittle2(acf, B[lag-1], A[lag-1], lag, "back", B[lag], KB, EA); |
| |
| copy_array(EA, subarray(v_forward,lag-1)); |
| copy_array(EB, subarray(v_back,lag-1)); |
| |
| copy_array(KA, subarray(p_forward,lag)); |
| copy_array(KB, subarray(p_back,lag)); |
| |
| } |
| |
| tmp = make_zero_matrix(nser,nser); |
| |
| matrix_prod(KB,KA, 1, 1, tmp); |
| array_op(id, tmp, '-', tmp); |
| matrix_prod(EA, tmp, 0, 0, subarray(v_forward, nlag)); |
| |
| vmaxset(vmax); |
| |
| } |
| |
| static void whittle2 (Array acf, Array Aold, Array Bold, int lag, |
| char *direction, Array A, Array K, Array E) |
| { |
| |
| int d, i, nser=DIM(acf)[1]; |
| const void *vmax; |
| Array beta, tmp, id; |
| |
| d = strcmp(direction, "forward") == 0; |
| |
| vmax = vmaxget(); |
| |
| beta = make_zero_matrix(nser,nser); |
| tmp = make_zero_matrix(nser, nser); |
| id = make_identity_matrix(nser); |
| |
| set_array_to_zero(E); |
| copy_array(id, subarray(A,0)); |
| |
| for(i = 0; i < lag; i++) { |
| matrix_prod(subarray(acf,lag - i), subarray(Aold,i), d, 1, tmp); |
| array_op(beta, tmp, '+', beta); |
| matrix_prod(subarray(acf,i), subarray(Bold,i), d, 1, tmp); |
| array_op(E, tmp, '+', E); |
| } |
| qr_solve(E, beta, K); |
| transpose_matrix(K,K); |
| for (i = 1; i <= lag; i++) { |
| matrix_prod(K, subarray(Bold,lag - i), 0, 0, tmp); |
| array_op(subarray(Aold,i), tmp, '-', subarray(A,i)); |
| } |
| |
| vmaxset(vmax); |
| } |