2 src/math/ts/innovations.c
4 Copyright (C) 2006 Free Software Foundation, Inc. Written by Jason H. Stover.
6 This program is free software; you can redistribute it and/or modify it under
7 the terms of the GNU General Public License as published by the Free
8 Software Foundation; either version 2 of the License, or (at your option)
11 This program is distributed in the hope that it will be useful, but WITHOUT
12 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
16 You should have received a copy of the GNU General Public License along with
17 this program; if not, write to the Free Software Foundation, Inc., 51
18 Franklin Street, Fifth Floor, Boston, MA 02111-1307, USA.
21 Find preliminary ARMA coefficients via the innovations algorithm.
22 Also compute the sample mean and covariance matrix for each series.
26 P. J. Brockwell and R. A. Davis. Time Series: Theory and
27 Methods. Second edition. Springer. New York. 1991. ISBN
28 0-387-97429-6. Sections 5.2, 8.3 and 8.4.
32 #include <gsl/gsl_matrix.h>
33 #include <gsl/gsl_vector.h>
34 #include <gsl/gsl_math.h>
36 #include <libpspp/alloc.h>
37 #include <libpspp/compiler.h>
38 #include <math/coefficient.h>
39 #include <math/ts/innovations.h>
42 get_mean (const gsl_matrix *data,
43 struct innovations_estimate **est)
51 for (n = 0; n < data->size2; n++)
56 for (i = 0; i < data->size1; i++)
58 for (n = 0; n < data->size2; n++)
60 tmp = gsl_matrix_get (data, i, n);
64 d = (tmp - est[n]->mean) / est[n]->n_obs;
71 update_cov (struct innovations_estimate **est, gsl_vector_const_view x,
72 gsl_vector_const_view y, size_t lag)
78 for (j = 0; j < x.vector.size; j++)
80 xj = gsl_vector_get (&x.vector, j);
81 yj = gsl_vector_get (&y.vector, j);
88 *(est[j]->cov + lag) += xj * yj;
94 get_covariance (const gsl_matrix *data,
95 struct innovations_estimate **est, size_t max_lag)
102 assert (data != NULL);
103 assert (est != NULL);
105 for (j = 0; j < data->size2; j++)
107 for (lag = 0; lag <= max_lag; lag++)
109 *(est[j]->cov + lag) = 0.0;
113 The rows are in the outer loop because a gsl_matrix is stored in
116 for (i = 0; i < data->size1; i++)
118 for (lag = 0; lag <= max_lag && lag < data->size1 - i; lag++)
120 update_cov (est, gsl_matrix_const_row (data, i),
121 gsl_matrix_const_row (data, i + lag), lag);
124 for (j = 0; j < data->size2; j++)
126 for (lag = 0; lag <= max_lag; lag++)
128 *(est[j]->cov + lag) /= est[j]->n_obs;
136 innovations_convolve (double *x, double *y, struct innovations_estimate *est,
142 assert (x != NULL && y != NULL);
143 assert (est != NULL);
144 assert (est->scale != NULL);
146 for (k = 0; k < i; k++)
148 result += x[k] * y[k] * est->scale[i-k-1];
153 innovations_update_scale (struct innovations_estimate *est, double *theta,
160 if (i < (size_t) est->max_lag)
162 result = est->cov[0];
163 for (j = 0; j < i; j++)
166 result -= theta[k] * theta[k] * est->scale[j];
168 est->scale[i] = result;
172 init_theta (double **theta, size_t max_lag)
177 for (i = 0; i < max_lag; i++)
179 for (j = 0; j <= i; j++)
186 innovations_update_coeff (double **theta, struct innovations_estimate *est,
193 for (i = 0; i < max_lag; i++)
195 theta[i][i] = est->cov[i+1] / est->scale[0];
196 for (j = 1; j <= i; j++)
199 theta[i][k] = (est->cov[k+1] -
200 innovations_convolve (theta[i] + k + 1, theta[j - 1], est, j))
203 innovations_update_scale (est, theta[i], i + 1);
207 get_coef (const gsl_matrix *data,
208 struct innovations_estimate **est, size_t max_lag)
214 theta = xnmalloc (max_lag, sizeof (*theta));
215 for (i = 0; i < max_lag; i++)
217 theta[i] = xnmalloc (max_lag, sizeof (**(theta + i)));
220 for (n = 0; n < data->size2; n++)
222 init_theta (theta, max_lag);
223 innovations_update_scale (est[n], theta[0], 0);
224 innovations_update_coeff (theta, est[n], max_lag);
225 /* Copy the final row of coefficients into EST->COEFF.*/
226 for (i = 0; i < max_lag; i++)
229 The order of storage here means that the best predicted value
230 for the time series is computed as follows:
232 Let X[m], X[m-1],... denote the original series.
233 Let X_hat[0] denote the best predicted value of X[0],
234 X_hat[1] denote the projection of X[1] onto the subspace
235 spanned by {X[0] - X_hat[0]}. Let X_hat[m] denote the
236 projection of X[m] onto the subspace spanned by {X[m-1] - X_hat[m-1],
237 X[m-2] - X_hat[m-2],...,X[0] - X_hat[0]}.
239 Then X_hat[m] = est->coeff[m-1] * (X[m-1] - X_hat[m-1])
240 + est->coeff[m-1] * (X[m-2] - X_hat[m-2])
242 + est->coeff[m-max_lag] * (X[m - max_lag] - X_hat[m - max_lag])
244 pspp_coeff_set_estimate (est[n]->coeff[i], theta[max_lag - 1][i]);
248 for (i = 0; i < max_lag; i++)
256 innovations_struct_init (struct innovations_estimate *est,
257 const struct design_matrix *dm,
263 /* COV[0] stores the lag 0 covariance (i.e., the variance), COV[1]
264 holds the lag-1 covariance, etc.
266 est->cov = xnmalloc (lag + 1, sizeof (*est->cov));
267 est->scale = xnmalloc (lag + 1, sizeof (*est->scale));
268 est->coeff = xnmalloc (lag, sizeof (*est->coeff)); /* No intercept. */
271 The loop below is an unusual use of PSPP_COEFF_INIT(). In a
272 typical model, one column of a DESIGN_MATRIX has one
273 coefficient. But in a time-series model, one column has many
276 for (j = 0; j < lag; j++)
278 pspp_coeff_init (est->coeff + j, dm);
280 est->max_lag = (double) lag;
283 The mean is subtracted from the original data before computing the
284 coefficients. The mean is NOT added back, so if you want to predict
285 a new value, you must add the mean to X_hat[m] to get the correct
289 subtract_mean (gsl_matrix *m, struct innovations_estimate **est)
295 for (i = 0; i < m->size1; i++)
297 for (j = 0; j < m->size2; j++)
299 tmp = gsl_matrix_get (m, i, j) - est[j]->mean;
300 gsl_matrix_set (m, i, j, tmp);
304 struct innovations_estimate **
305 pspp_innovations (const struct design_matrix *dm, size_t lag)
307 struct innovations_estimate **est;
310 est = xnmalloc (dm->m->size2, sizeof *est);
311 for (i = 0; i < dm->m->size2; i++)
313 est[i] = xmalloc (sizeof *est[i]);
314 /* est[i]->variable = vars[i]; */
315 innovations_struct_init (est[i], dm, lag);
318 get_mean (dm->m, est);
319 subtract_mean (dm->m, est);
320 get_covariance (dm->m, est, lag);
321 get_coef (dm->m, est, lag);
327 pspp_innovations_free_one (struct innovations_estimate *est)
331 assert (est != NULL);
332 for (i = 0; i < (size_t) est->max_lag; i++)
334 pspp_coeff_free (est->coeff[i]);
341 void pspp_innovations_free (struct innovations_estimate **est, size_t n)
345 assert (est != NULL);
346 for (i = 0; i < n; i++)
348 pspp_innovations_free_one (est[i]);