From 9368d05aae3a64fd81b237dcf77fb1376af2fefd Mon Sep 17 00:00:00 2001 From: Jason Stover Date: Thu, 25 May 2006 22:15:03 +0000 Subject: [PATCH] new file --- src/math/time-series/ChangeLog | 4 + src/math/time-series/innovations.c | 123 +++++++++++++++++++++++++++++ 2 files changed, 127 insertions(+) create mode 100644 src/math/time-series/ChangeLog create mode 100644 src/math/time-series/innovations.c diff --git a/src/math/time-series/ChangeLog b/src/math/time-series/ChangeLog new file mode 100644 index 00000000..bd14725a --- /dev/null +++ b/src/math/time-series/ChangeLog @@ -0,0 +1,4 @@ +2006-05-25 Jason Stover + + * innovations.c: New file + diff --git a/src/math/time-series/innovations.c b/src/math/time-series/innovations.c new file mode 100644 index 00000000..066530d2 --- /dev/null +++ b/src/math/time-series/innovations.c @@ -0,0 +1,123 @@ +/* + src/math/time-series/arma/innovations.c + + Copyright (C) 2006 Free Software Foundation, Inc. Written by Jason H. Stover. + + 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, write to the Free Software Foundation, Inc., 51 + Franklin Street, Fifth Floor, Boston, MA 02111-1307, USA. + */ +/* + Find preliminary ARMA coefficients via the innovations algorithm. + Also compute the sample mean and covariance matrix for each series. + + Reference: + + P. J. Brockwell and R. A. Davis. Time Series: Theory and + Methods. Second edition. Springer. New York. 1991. ISBN + 0-387-97429-6. Sections 5.2, 8.3 and 8.4. + */ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#define _(msgid) gettext (msgid) + +static void +get_mean_variance (size_t n_vars, const struct casefile *cf, + struct innovations_estimate **est) + +{ + struct casereader *r; + struct ccase *c; + size_t n; + double *x; + double d; + double tmp; + double variance; + + x = xnmalloc (n_vars, sizeof *j); + + for (n = 0; n < n_vars; n++) + { + x[n] = 2.0; + est[n]->mean = 0.0; + est[n]->variance = 0.0; + } + for (r = casefile_get_reader (cf); casereader_read (r, &c); + case_destroy (&c)) + { + for (n = 0; n < n_vars; n++) + { + if (!mv_is_value_missing (&v->miss, val)) + { + tmp = case_data (&c, est[n]->variable->fv); + d = (tmp - est[n]->mean) / x[n]; + est[n]->mean += d; + est[n]->variance += x[n] * x[n] * d * d; + x[n] += 1.0; + } + } + } + for (n = 0; n < n_vars; n++) + { + est[n]->variance /= x[n]; + } + free (x); +} + +struct innovations_estimate ** pspp_innovations (const struct variable **vars, size_t *n_vars, + size_t max_lag, const struct casefile *cf) +{ + struct innovations_estimate **est; + struct casereader *r; + struct ccase *c; + size_t i; + + est = xnmalloc (*n_vars, sizeof *est); + for (i = 0; i < *n_vars; i++) + { + if (vars[i]->type == NUMERIC) + { + est[i] = xmalloc (sizeof **est); + est[i]->variable = vars[i]; + est[i]->mean = 0.0; + est[i]->variance = 0.0; + est[i]->cov = gsl_matrix_calloc (max_lag, max_lag); + est[i]->coeff = xnmalloc (max_lag, sizeof (*est[i]->coeff)); + for (j = 0; j < max_lag; j++) + { + est[i]->coeff + j = xmalloc (sizeof (*(est[i]->coeff + j))); + } + } + else + { + *n_vars--; + msg (MW, _("Cannot compute autocovariance for a non-numeric variable %s"), + var_to_string (vars[i])); + } + } + + /* + First data pass to get the mean and variance. + */ + get_mean_variance (*n_vars, cf, est); +} -- 2.30.2