src/math/ts/innovations.c

   1 /*
   2   src/math/ts/innovations.c
   3
   4   Copyright (C) 2006 Free Software Foundation, Inc. Written by Jason H. Stover.
   5
   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)
   9   any later version.
  10
  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
  14   more details.
  15
  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.
  19  */
  20 /*
  21   Find preliminary ARMA coefficients via the innovations algorithm.
  22   Also compute the sample mean and covariance matrix for each series.
  23
  24   Reference:
  25
  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.
  29  */
  30
  31 #include <gsl/gsl_matrix.h>
  32 #include <gsl/gsl_vector.h>
  33 #include <gsl/gsl_math.h>
  34 #include <stdlib.h>
  35 #include <libpspp/alloc.h>
  36 #include <libpspp/compiler.h>
  37 #include <math/coefficient.h>
  38 #include <math/ts/innovations.h>
  39
  40 static void
  41 get_mean (const gsl_matrix *data,
  42           struct innovations_estimate **est)
  43
  44 {
  45   size_t n;
  46   size_t i;
  47   double d;
  48   double tmp;
  49
  50   for (n = 0; n < data->size2; n++)
  51     {
  52       est[n]->n_obs = 0.0;
  53       est[n]->mean = 0.0;
  54     }
  55   for (i = 0; i < data->size1; i++)
  56     {
  57       for (n = 0; n < data->size2; n++)
  58         {
  59           tmp = gsl_matrix_get (data, i, n);
  60           if (!gsl_isnan (tmp))
  61             {
  62               est[n]->n_obs += 1.0;
  63               d = (tmp - est[n]->mean) / est[n]->n_obs;
  64               est[n]->mean += d;
  65             }
  66         }
  67     }
  68 }
  69 static void
  70 update_cov (struct innovations_estimate **est, gsl_vector_const_view x,
  71             gsl_vector_const_view y, size_t lag)
  72 {
  73   size_t j;
  74   double xj;
  75   double yj;
  76
  77   for (j = 0; j < x.vector.size; j++)
  78     {
  79       xj = gsl_vector_get (&x.vector, j);
  80       yj = gsl_vector_get (&y.vector, j);
  81       if (!gsl_isnan (xj))
  82         {
  83           if (!gsl_isnan (yj))
  84             {
  85               xj -= est[j]->mean;
  86               yj -= est[j]->mean;
  87               *(est[j]->cov + lag) += xj * yj;
  88             }
  89         }
  90     }
  91 }
  92 static int
  93 get_covariance (const gsl_matrix *data,
  94                 struct innovations_estimate **est, size_t max_lag)
  95 {
  96   size_t lag;
  97   size_t j;
  98   size_t i;
  99   int rc = 1;
 100
 101   assert (data != NULL);
 102   assert (est != NULL);
 103
 104   for (j = 0; j < data->size2; j++)
 105     {
 106       for (lag = 0; lag <= max_lag; lag++)
 107         {
 108           *(est[j]->cov + lag) = 0.0;
 109         }
 110     }
 111   /*
 112     The rows are in the outer loop because a gsl_matrix is stored in
 113     row-major order.
 114    */
 115   for (i = 0; i < data->size1; i++)
 116     {
 117       for (lag = 0; lag < max_lag && lag < data->size1 - i; lag++)
 118         {
 119           update_cov (est, gsl_matrix_const_row (data, i),
 120                       gsl_matrix_const_row (data, i + lag), lag);
 121         }
 122     }
 123   for (j = 0; j < data->size2; j++)
 124     {
 125       for (lag = 0; lag <= max_lag; lag++)
 126         {
 127           *(est[j]->cov + lag) /= est[j]->n_obs;
 128         }
 129     }
 130
 131   return rc;
 132 }
 133
 134 static double
 135 innovations_convolve (double **theta, struct innovations_estimate *est,
 136                       int i, int j)
 137 {
 138   int k;
 139   double result = 0.0;
 140
 141   for (k = 0; k < j; k++)
 142     {
 143       result += theta[i-1][i-k-1] * theta[j][j-k-1] * est->scale[k];
 144     }
 145   return result;
 146 }
 147 static void
 148 innovations_update_scale (struct innovations_estimate *est, double *theta,
 149                           size_t i)
 150 {
 151   double result = 0.0;
 152   size_t j;
 153   size_t k;
 154
 155   if (i < (size_t) est->max_lag)
 156     {
 157       result = est->cov[0];
 158       for (j = 0; j < i; j++)
 159         {
 160           k = i - j - 1;
 161           result -= theta[k] * theta[k] * est->scale[j];
 162         }
 163       est->scale[i] = result;
 164     }
 165 }
 166 static void
 167 init_theta (double **theta, size_t max_lag)
 168 {
 169   size_t i;
 170   size_t j;
 171
 172   for (i = 0; i < max_lag; i++)
 173     {
 174       for (j = 0; j <= i; j++)
 175         {
 176           theta[i][j] = 0.0;
 177         }
 178     }
 179 }
 180 static void
 181 innovations_update_coeff (double **theta, struct innovations_estimate *est,
 182                           size_t max_lag)
 183 {
 184   size_t i;
 185   size_t j;
 186   size_t k;
 187
 188   for (i = 0; i < max_lag; i++)
 189     {
 190       for (j = 0; j <= i; j++)
 191         {
 192           k = i - j;
 193           theta[i][k] = (est->cov[k] -
 194             innovations_convolve (theta, est, i, j))
 195             / est->scale[k];
 196         }
 197       innovations_update_scale (est, theta[i], i + 1);
 198     }
 199 }
 200 static void
 201 get_coef (const gsl_matrix *data,
 202           struct innovations_estimate **est, size_t max_lag)
 203 {
 204   size_t i;
 205   size_t n;
 206   double **theta;
 207
 208   theta = xnmalloc (max_lag, sizeof (*theta));
 209   for (i = 0; i < max_lag; i++)
 210     {
 211       theta[i] = xnmalloc (max_lag, sizeof (**(theta + i)));
 212     }
 213
 214   for (n = 0; n < data->size2; n++)
 215     {
 216       init_theta (theta, max_lag);
 217       innovations_update_scale (est[n], theta[0], 0);
 218       innovations_update_coeff (theta, est[n], max_lag);
 219       /* Copy the final row of coefficients into EST->COEFF.*/
 220       for (i = 0; i < max_lag; i++)
 221         {
 222           /*
 223             The order of storage here means that the best predicted value
 224             for the time series is computed as follows:
 225
 226             Let X[m], X[m-1],... denote the original series.
 227             Let X_hat[0] denote the best predicted value of X[0],
 228             X_hat[1] denote the projection of X[1] onto the subspace
 229             spanned by {X[0] - X_hat[0]}. Let X_hat[m] denote the
 230             projection of X[m] onto the subspace spanned by {X[m-1] - X_hat[m-1],
 231             X[m-2] - X_hat[m-2],...,X[0] - X_hat[0]}.
 232
 233             Then X_hat[m] = est->coeff[m-1] * (X[m-1] - X_hat[m-1])
 234                           + est->coeff[m-1] * (X[m-2] - X_hat[m-2])
 235                           ...
 236                           + est->coeff[m-max_lag] * (X[m - max_lag] - X_hat[m - max_lag])
 237
 238             (That is what X_hat[m] SHOULD be, anyway. These routines need
 239             to be tested.)
 240            */
 241           pspp_coeff_set_estimate (est[n]->coeff[i], theta[max_lag - 1][i]);
 242         }
 243     }
 244
 245   for (i = 0; i < max_lag; i++)
 246     {
 247       free (theta[i]);
 248     }
 249   free (theta);
 250 }
 251
 252 static void
 253 innovations_struct_init (struct innovations_estimate *est,
 254                          const struct design_matrix *dm,
 255                          size_t lag)
 256 {
 257   size_t j;
 258
 259   est->mean = 0.0;
 260   /* COV[0] stores the lag 0 covariance (i.e., the variance), COV[1]
 261      holds the lag-1 covariance, etc.
 262    */
 263   est->cov = xnmalloc (lag + 1, sizeof (*est->cov));
 264   est->scale = xnmalloc (lag + 1, sizeof (*est->scale));
 265   est->coeff = xnmalloc (lag, sizeof (*est->coeff)); /* No intercept. */
 266
 267   /*
 268     The loop below is an unusual use of PSPP_COEFF_INIT(). In a
 269     typical model, one column of a DESIGN_MATRIX has one
 270     coefficient. But in a time-series model, one column has many
 271     coefficients.
 272    */
 273   for (j = 0; j < lag; j++)
 274     {
 275       pspp_coeff_init (est->coeff + j, dm);
 276     }
 277   est->max_lag = (double) lag;
 278 }
 279
 280 struct innovations_estimate **
 281 pspp_innovations (const struct design_matrix *dm, size_t lag)
 282 {
 283   struct innovations_estimate **est;
 284   size_t i;
 285
 286   est = xnmalloc (dm->m->size2, sizeof *est);
 287   for (i = 0; i < dm->m->size2; i++)
 288     {
 289       est[i] = xmalloc (sizeof *est[i]);
 290 /*       est[i]->variable = vars[i]; */
 291       innovations_struct_init (est[i], dm, lag);
 292     }
 293
 294   get_mean (dm->m, est);
 295   get_covariance (dm->m, est, lag);
 296   get_coef (dm->m, est, lag);
 297
 298   return est;
 299 }
 300
 301 static void
 302 pspp_innovations_free_one (struct innovations_estimate *est)
 303 {
 304   size_t i;
 305
 306   assert (est != NULL);
 307   for (i = 0; i < (size_t) est->max_lag; i++)
 308     {
 309       pspp_coeff_free (est->coeff[i]);
 310     }
 311   free (est->scale);
 312   free (est->cov);
 313   free (est);
 314 }
 315
 316 void pspp_innovations_free (struct innovations_estimate **est, size_t n)
 317 {
 318   size_t i;
 319
 320   assert (est != NULL);
 321   for (i = 0; i < n; i++)
 322     {
 323       pspp_innovations_free_one (est[i]);
 324     }
 325   free (est);
 326 }