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_variance (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 = 2.0;
  53       est[n]->mean = 0.0;
  54       est[n]->variance = 0.0;
  55     }
  56   for (i = 0; i < data->size1; i++)
  57     {
  58       for (n = 0; n < data->size2; n++)
  59         {
  60           tmp = gsl_matrix_get (data, i, n);
  61           if (!gsl_isnan (tmp))
  62             {
  63               d = (tmp - est[n]->mean) / est[n]->n_obs;
  64               est[n]->mean += d;
  65               est[n]->variance += est[n]->n_obs * est[n]->n_obs * d * d;
  66               est[n]->n_obs += 1.0;
  67             }
  68         }
  69     }
  70   for (n = 0; n < data->size2; n++)
  71     {
  72       /* Maximum likelihood estimate of the variance. */
  73       est[n]->variance /= est[n]->n_obs;
  74     }
  75 }
  76
  77 static int
  78 get_covariance (const gsl_matrix *data,
  79                 struct innovations_estimate **est, size_t max_lag)
  80 {
  81   size_t lag;
  82   size_t j;
  83   size_t i;
  84   double x;
  85   double y;
  86   int rc = 1;
  87
  88   assert (data != NULL);
  89   assert (est != NULL);
  90
  91   for (i = 0; i < data->size1; i++)
  92     {
  93       for (j = 0; j < data->size2; j++)
  94         {
  95           x = gsl_matrix_get (data, i, j);
  96
  97           if (!gsl_isnan (x))
  98             {
  99               x -= est[j]->mean;
 100               for (lag = 1; lag <= max_lag && lag < data->size1 - max_lag; lag++)
 101                 {
 102                   y = gsl_matrix_get (data, i + lag, j);
 103                   if (!gsl_isnan (y))
 104                     {
 105                       y -= est[j]->mean;
 106                       *(est[j]->cov + lag) += y * x;
 107                       est[i]->n_obs += 1.0;
 108                     }
 109                 }
 110             }
 111         }
 112     }
 113   for (lag = 0; lag <= max_lag && lag < data->size1 - max_lag; lag++)
 114     {
 115       for (j = 0; j < data->size2; j++)
 116         {
 117           *(est[j]->cov + lag) /= (est[j]->n_obs - lag);
 118         }
 119     }
 120   return rc;
 121 }
 122 static double
 123 innovations_convolve (double **theta, struct innovations_estimate *est,
 124                       int i, int j)
 125 {
 126   int k;
 127   double result = 0.0;
 128
 129   for (k = 0; k < i; k++)
 130     {
 131       result += theta[i-1][i-k-1] * theta[j-1][j-k-1] * est->scale[k];
 132     }
 133   return result;
 134 }
 135 static void
 136 innovations_update_scale (struct innovations_estimate *est, double *theta,
 137                           size_t i)
 138 {
 139   double result = 0.0;
 140   size_t j;
 141   size_t k;
 142
 143
 144   result = est->cov[0];
 145   for (j = 0; j < i; j++)
 146     {
 147       k = i - j;
 148       result -= theta[k] * theta[k] * est->scale[j];
 149     }
 150   est->scale[i] = result;
 151 }
 152
 153 static void
 154 get_coef (const gsl_matrix *data,
 155           struct innovations_estimate **est, size_t max_lag)
 156 {
 157   size_t j;
 158   size_t i;
 159   size_t k;
 160   size_t n;
 161   double v;
 162   double **theta;
 163
 164   theta = xnmalloc (max_lag, sizeof (*theta));
 165   for (i = 0; i < max_lag; i++)
 166     {
 167       theta[i] = xnmalloc (i+1, sizeof (theta[i]));
 168
 169     }
 170   for (n = 0; n < data->size2; n++)
 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       innovations_update_scale (est[n], theta[0], 0);
 180       for (i = 0; i < max_lag; i++)
 181         {
 182           v = est[n]->cov[i];
 183           for (j = 0; j < i; j++)
 184             {
 185               k = i - j;
 186               theta[i-1][k-1] = est[n]->cov[k] -
 187                 innovations_convolve (theta, est[n], i, j);
 188             }
 189           innovations_update_scale (est[n], theta[i], i);
 190         }
 191       /* Copy the final row of coefficients into EST->COEFF.*/
 192       for (i = 0; i < max_lag; i++)
 193         {
 194           /*
 195             The order of storage here means that the best predicted value
 196             for the time series is computed as follows:
 197
 198             Let X[m], X[m-1],... denote the original series.
 199             Let X_hat[0] denote the best predicted value of X[0],
 200             X_hat[1] denote the projection of X[1] onto the subspace
 201             spanned by {X[0] - X_hat[0]}. Let X_hat[m] denote the
 202             projection of X[m] onto the subspace spanned by {X[m-1] - X_hat[m-1],
 203             X[m-2] - X_hat[m-2],...,X[0] - X_hat[0]}.
 204
 205             Then X_hat[m] = est->coeff[m-1] * (X[m-1] - X_hat[m-1])
 206                           + est->coeff[m-1] * (X[m-2] - X_hat[m-2])
 207                           ...
 208                           + est->coeff[m-max_lag] * (X[m - max_lag] - X_hat[m - max_lag])
 209
 210             (That is what X_hat[m] SHOULD be, anyway. These routines need
 211             to be tested.)
 212            */
 213           pspp_coeff_set_estimate (est[n]->coeff[i], theta[max_lag - 1][i]);
 214         }
 215     }
 216   for (i = 0; i < max_lag; i++)
 217     {
 218       free (theta[i]);
 219     }
 220   free (theta);
 221 }
 222
 223 struct innovations_estimate **
 224 pspp_innovations (const gsl_matrix *data, size_t lag)
 225 {
 226   struct innovations_estimate **est;
 227   size_t i;
 228   size_t j;
 229
 230   est = xnmalloc (data->size2, sizeof *est);
 231   for (i = 0; i < data->size2; i++)
 232     {
 233       est[i] = xmalloc (sizeof **est);
 234 /*       est[i]->variable = vars[i]; */
 235       est[i]->mean = 0.0;
 236       est[i]->variance = 0.0;
 237       est[i]->cov = xnmalloc (lag, sizeof (*est[i]->cov));
 238       est[i]->scale = xnmalloc (lag, sizeof (*est[i]->scale));
 239       est[i]->coeff = xnmalloc (lag, sizeof (*est[i]->coeff));
 240       for (j = 0; j < lag; j++)
 241         {
 242           est[i]->coeff[j] = xmalloc (sizeof (*(est[i]->coeff + j)));
 243         }
 244     }
 245
 246   get_mean_variance (data, est);
 247   get_covariance (data, est, lag);
 248   get_coef (data, est, lag);
 249
 250   return est;
 251 }