[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Pspp-cvs] experimental/src/math/ts acf.c arma.c dacf.c in...
|
From: |
Jason H Stover |
|
Subject: |
[Pspp-cvs] experimental/src/math/ts acf.c arma.c dacf.c in... |
|
Date: |
Wed, 24 Jan 2007 17:28:35 +0000 |
CVSROOT: /sources/pspp
Module name: experimental
Changes by: Jason H Stover <jstover> 07/01/24 17:28:35
Added files:
src/math/ts : acf.c arma.c dacf.c innovations.c likelihood.c
predict.c tmp.c
Log message:
initial version
CVSWeb URLs:
http://cvs.savannah.gnu.org/viewcvs/experimental/src/math/ts/acf.c?cvsroot=pspp&rev=1.1
http://cvs.savannah.gnu.org/viewcvs/experimental/src/math/ts/arma.c?cvsroot=pspp&rev=1.1
http://cvs.savannah.gnu.org/viewcvs/experimental/src/math/ts/dacf.c?cvsroot=pspp&rev=1.1
http://cvs.savannah.gnu.org/viewcvs/experimental/src/math/ts/innovations.c?cvsroot=pspp&rev=1.1
http://cvs.savannah.gnu.org/viewcvs/experimental/src/math/ts/likelihood.c?cvsroot=pspp&rev=1.1
http://cvs.savannah.gnu.org/viewcvs/experimental/src/math/ts/predict.c?cvsroot=pspp&rev=1.1
http://cvs.savannah.gnu.org/viewcvs/experimental/src/math/ts/tmp.c?cvsroot=pspp&rev=1.1
Patches:
Index: acf.c
===================================================================
RCS file: acf.c
diff -N acf.c
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ acf.c 24 Jan 2007 17:28:35 -0000 1.1
@@ -0,0 +1,262 @@
+/*
+ src/math/ts/acf.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.
+ */
+/*
+ Compute the autocovariance function of an ARMA (p,q) process. The
+ process must be causal and invertible.
+ */
+#include <assert.h>
+#include "acf.h"
+#include "arma-optimizer.h"
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_blas.h>
+#include <stdlib.h>
+#include <libpspp/alloc.h>
+#include <libpspp/compiler.h>
+
+/*
+ Find the coefficients for the infinite-order MA process.
+ */
+void
+set_inf_ma_coef (struct pspp_arma_state *s)
+{
+ size_t i;
+ size_t j;
+ pspp_arma *x;
+ double tmp;
+
+ x = s->arma;
+ assert (s->inf_ma_coeff != NULL);
+ gsl_vector_set (s->inf_ma_coeff, 0, 1.0);
+
+ for (i = 1; i <= x->ar_order; i++)
+ {
+ tmp = arma_get_ma_est (x, i - 1);
+ for (j = 0; j < i; j++)
+ {
+ tmp += gsl_vector_get (s->inf_ma_coeff, j) *
+ arma_get_ar_est (x, i - j - 1);
+ }
+ gsl_vector_set (s->inf_ma_coeff, i, tmp);
+ }
+}
+static double
+get_phi (int i, const pspp_arma *x)
+{
+ /*
+ Be careful here: we are using the sign of i,
+ so don't use a size_t.
+ */
+ if (i == 0)
+ {
+ return -1.0;
+ }
+ else if (i < 0)
+ {
+ return 0.0;
+ }
+ return arma_get_ar_est (x, (size_t) i - 1);
+}
+
+void
+get_matrix_elements (gsl_matrix *A, const pspp_arma *x)
+{
+ size_t i;
+ size_t j;
+ double phi1;
+ double phi2;
+
+ gsl_matrix_set (A, 0, 0, 1.0);
+ for (j = 1; j < A->size2; j++)
+ {
+ phi1 = -1 * arma_get_ar_est (x, j - 1);
+ gsl_matrix_set (A, 0, j, phi1);
+ gsl_matrix_set (A, j, 0, phi1);
+ }
+ for (i = 1; i < A->size1; i++)
+ {
+ for (j = 1; j <= i; j++)
+ {
+ phi1 = -1.0 * get_phi (i-j, x);
+ phi2 = -1.0 * get_phi (i+j, x);
+ gsl_matrix_set (A, i, j, phi1 + phi2);
+ }
+ for (j = i + 1; j < x->ar_order; j++)
+ {
+ phi1 = arma_get_ar_est (x, 2 * i);
+ gsl_matrix_set (A, i, j, phi2);
+ }
+ }
+}
+static double
+get_vector_element (const gsl_vector *psi, const pspp_arma *x, size_t i)
+{
+ size_t j;
+ double result;
+ double tmp;
+
+ result = (i == 0) ? 1.0 : 0.0;
+ for (j = 0; j < x->ma_order; j++)
+ {
+ tmp = ((j + 1) < i) ? 0.0 : gsl_vector_get (psi, j + 1 - i);
+ result += arma_get_ma_est (x, j) * tmp;
+ }
+
+ return result;
+}
+static void
+get_vector_elements (gsl_vector *b, const struct pspp_arma_state *s)
+{
+ size_t i;
+ double element;
+ pspp_arma *x;
+
+ x = s->arma;
+ for (i = 0; i < b->size; i++)
+ {
+ element = get_vector_element (s->inf_ma_coeff, x, i);
+ gsl_vector_set (b, i, element * x->mse);
+ }
+}
+#if 0
+/*
+ The acf is found by first solving A \gamma = b, where
+ A is an ar->order-by-ar->order matrix, \gamma is a column
+ vector with entry i equal to the acf at lag i, and b is
+ column vector containing a convolution of ar->ma_coef and
+ the infinite-order moving average coefficients.
+ */
+void
+pspp_acf_vector (struct pspp_arma_state *s)
+{
+ gsl_vector *b;
+ gsl_matrix *A;
+ gsl_vector *tau;
+ pspp_arma *x = NULL;
+
+ x = s->arma;
+ assert (x != NULL);
+ set_inf_ma_coef (s);
+ A = gsl_matrix_calloc (x->ar_order + 1, x->ar_order + 1);
+ get_matrix_elements (A, x);
+ b = gsl_vector_calloc (x->ar_order + 1);
+ get_vector_elements (b, (const struct pspp_arma_state *) s);
+ tau = gsl_vector_alloc (A->size1);
+ gsl_linalg_QR_decomp (A, tau);
+ assert (x->acf != NULL);
+ assert (x->acf->size == x->ar_order + 1);
+ gsl_linalg_QR_solve ((const gsl_matrix *) A, (const gsl_vector *) tau,
+ (const gsl_vector *) b, x->acf);
+ gsl_matrix_free (A);
+ gsl_vector_free (tau);
+ gsl_vector_free (b);
+}
+#endif
+void
+pspp_acf_vector (struct pspp_arma_state *s)
+{
+ size_t i;
+ size_t j;
+ pspp_arma *x = NULL;
+ gsl_vector *shift_coeff;
+ double tmp;
+
+ x = s->arma;
+ assert (x != NULL);
+ set_inf_ma_coef (s);
+ assert (x->acf != NULL);
+ assert (x->acf->size == x->ar_order + 1);
+ j = 0;
+ shift_coeff = gsl_vector_alloc (s->inf_ma_coeff->size);
+ for (i = 0; i < x->acf->size; i++)
+ {
+ gsl_vector_set_zero (shift_coeff);
+ for (j = i; j < s->inf_ma_coeff->size; j++)
+ {
+ gsl_vector_set (shift_coeff, j - i, gsl_vector_get (s->inf_ma_coeff,
i));
+ }
+ gsl_blas_ddot (s->inf_ma_coeff, shift_coeff, &tmp);
+ gsl_vector_set (x->acf, i, tmp);
+ }
+ gsl_vector_scale (x->acf, x->mse);
+}
+double
+get_next_acf (pspp_arma *x, gsl_vector *past_acf)
+{
+ double result = 0.0;
+ size_t i;
+
+ assert (x != NULL);
+ assert (past_acf != NULL);
+
+ if (x->ar_order < past_acf->size)
+ {
+ return gsl_vector_get (x->acf, past_acf->size);
+ }
+ for (i = 0; i < x->ar_order; i++)
+ {
+ result += arma_get_ar_est (x, i) *
+ gsl_vector_get (past_acf, past_acf->size - i);
+ }
+ return result;
+}
+void
+copy_acf (pspp_arma *x, gsl_vector *acf)
+{
+ gsl_vector_memcpy (x->acf, acf);
+}
+double
+pspp_arma_acf (const pspp_arma *x, int lag)
+{
+ size_t i;
+ size_t j;
+ double result = 0.0;
+ gsl_vector *acf;
+
+ lag = abs (lag);
+ if (lag <= x->ar_order)
+ {
+ assert (x->acf != NULL);
+ result = gsl_vector_get (x->acf, (size_t) lag);
+ }
+ else
+ {
+ /*
+ Store the next ar_order of the autocovariance
+ function.
+ */
+ acf = gsl_vector_calloc (x->ar_order);
+ copy_acf ((pspp_arma *) x, acf);
+ i = x->ar_order;
+ while (i < (size_t) lag)
+ {
+ result = get_next_acf ((pspp_arma *) x, acf);
+ for (j = 1; j < acf->size; j++)
+ {
+ gsl_vector_set (acf, (j-1), gsl_vector_get (acf, j));
+ }
+ gsl_vector_set (acf, acf->size - 1, result);
+ i++;
+ }
+ result *= x->std * gsl_vector_get (acf, (size_t) lag);
+ }
+ return result;
+}
+
Index: arma.c
===================================================================
RCS file: arma.c
diff -N arma.c
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ arma.c 24 Jan 2007 17:28:35 -0000 1.1
@@ -0,0 +1,507 @@
+/*
+ src/math/ts/arma.c
+
+ Copyright (C) 2006 Free Software Foundation, Inc.
+
+ 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 ARMA coefficients via re-paramaterization and gradient-based
+ EM algorithm.
+
+ 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 <gsl/gsl_math.h>
+#include <gsl/gsl_multimin.h>
+#include <gsl/gsl_multifit.h>
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_errno.h>
+#include <stdlib.h>
+#include <libpspp/alloc.h>
+#include <libpspp/compiler.h>
+#include <math/ts/innovations.h>
+#include "arma.h"
+#include "acf.h"
+#include <assert.h>
+#include <libpspp/alloc.h>
+#include "arma-optimizer.h"
+#include <libpspp/alloc.h>
+static
+void initialize_minimizer (gsl_vector *, gsl_multimin_function_fdf *,
+ size_t, size_t);
+static void get_initial_estimates (struct pspp_arma_state *, const gsl_vector
*);
+static void coef_to_arma (const gsl_vector *, pspp_arma *);
+static void arma_to_coef (gsl_vector *, pspp_arma *);
+static void get_residuals (struct pspp_arma_state *, const gsl_vector *);
+
+/*
+ Allocate space for an ARMA model.
+ */
+pspp_arma *
+pspp_arma_alloc (size_t n_obs, size_t p, size_t q)
+{
+ pspp_arma *arma;
+
+ arma = xmalloc (sizeof (*arma));
+ assert (arma != NULL);
+ arma->n_obs = n_obs;
+ arma->ar_order = p;
+ arma->ma_order = q;
+ arma->ar_coeff = gsl_vector_calloc (arma->ar_order);
+ arma->ma_coeff = gsl_vector_calloc (arma->ma_order);
+ arma->cov = gsl_matrix_calloc (p + q, p + q);
+ arma->acf = gsl_vector_calloc (p + 1);
+ return arma;
+}
+
+double
+arma_get_ar_est (const pspp_arma *x, size_t i)
+{
+ double result = 0.0;
+
+ assert (x != NULL);
+ if (i <= x->ar_order)
+ {
+ result = gsl_vector_get (x->ar_coeff, i);
+ }
+ return result;
+}
+double
+arma_get_ma_est (const pspp_arma *x, size_t i)
+{
+ double result = 0.0;
+
+ assert (x != NULL);
+ if (i <= x->ma_order)
+ {
+ result = gsl_vector_get (x->ma_coeff, i);
+ }
+ return result;
+}
+static void
+get_scale (struct pspp_arma_state *as)
+{
+ size_t i;
+ size_t j;
+ size_t m;
+ double tmp;
+ double ma;
+
+ m = (as->arma->ar_order > as->arma->ma_order) ? as->arma->ar_order :
as->arma->ma_order;
+ for (i = 0; i < m; i++)
+ {
+ gsl_vector_set (as->scale, i, gsl_vector_get (as->arma->acf, 0));
+ }
+ for (i = m; i < as->scale->size; i++)
+ {
+ tmp = gsl_vector_get (as->arma->acf, 0);
+ for (j = 0; j < as->arma->ma_order; j++)
+ {
+ ma = arma_get_ma_est (as->arma, j);
+ tmp -= gsl_vector_get (as->scale, i - j - 1) * ma * ma;
+ assert (tmp > 0.0);
+ assert (!gsl_isnan (tmp));
+ }
+ gsl_vector_set (as->scale, i, tmp);
+ }
+}
+static void
+get_mse (struct pspp_arma_state *as)
+{
+ size_t i;
+ double tmp;
+
+ as->arma->mse = 0.0;
+ for (i = 0; i < as->residuals->size; i++)
+ {
+ tmp = gsl_vector_get (as->residuals, i);
+ as->arma->mse += tmp * tmp / gsl_vector_get (as->scale, i);
+ }
+ as->arma->mse /= as->residuals->size;
+ as->arma->std = sqrt (as->arma->mse);
+}
+/*
+ Initial values are derived from the innovations estimates
+ via least-squares regression.
+ */
+static void
+get_initial_estimates (struct pspp_arma_state *as, const gsl_vector *data)
+{
+ gsl_multifit_linear_workspace *ws;
+ pspp_arma *a;
+ gsl_vector *pred;
+ gsl_vector *resid;
+ gsl_vector_view y;
+ gsl_vector *coef;
+ gsl_matrix *obs;
+ gsl_matrix *cov;
+ double mse;
+ size_t i;
+ size_t j;
+ size_t order;
+ int rc;
+
+ a = as->arma;
+ as->in = pspp_innovations (data, NULL, NULL);
+ pred = pspp_innovations_predicted_values (as->in, data);
+ resid = pspp_innovations_residuals (as->in, data);
+ order = a->ar_order + a->ma_order;
+ assert (order < as->in->n_obs);
+ obs = gsl_matrix_alloc (as->in->n_obs - order, order);
+ for (j = 0; j < obs->size1; j++)
+ {
+ for (i = 0; i < a->ar_order; i++)
+ {
+ gsl_matrix_set (obs, j, i, gsl_vector_get (data, order - 1 - i + j));
+ }
+ for (i = a->ar_order; i < obs->size2; i++)
+ {
+ gsl_matrix_set (obs, j, i, gsl_vector_get (resid, order - 1 - i +
a->ar_order + j));
+ }
+ }
+ coef = gsl_vector_alloc (obs->size2);
+ cov = gsl_matrix_alloc (coef->size, coef->size);
+ y = gsl_vector_subvector (pred, order, obs->size1);
+ ws = gsl_multifit_linear_alloc (obs->size1, order);
+ rc = gsl_multifit_linear (obs, &y.vector, coef, cov, &mse, ws);
+ assert (rc == GSL_SUCCESS);
+ coef_to_arma (coef, a);
+ gsl_multifit_linear_free (ws);
+ a->mse = 0.0;
+ for (i = 0; i < resid->size; i++)
+ {
+ mse = gsl_vector_get (resid, i);
+ a->mse += mse * mse / gsl_vector_get (as->in->scale, i);
+ }
+ a->mse /= resid->size;
+ a->std = sqrt (a->mse);
+ gsl_vector_free (resid);
+ gsl_vector_free (pred);
+ gsl_vector_free (coef);
+ gsl_matrix_free (cov);
+}
+/*
+ Get the smallest subvector of size at least ORDER
+ with non-constant entries.
+*/
+static gsl_vector_view
+get_subvector (gsl_vector *y, size_t offset, size_t size)
+{
+ gsl_vector_view result;
+
+ assert ((size + offset) < y->size);
+ result = gsl_vector_subvector (y, offset, size);
+ while ((gsl_vector_max (&result.vector) - gsl_vector_min (&result.vector)) <
GSL_DBL_EPSILON)
+ {
+ size++;
+ result = gsl_vector_subvector (y, offset, size);
+ }
+ return result;
+}
+/*
+ COEF should store the AR and MA coefficients, in that
+ order.
+ */
+double
+pspp_arma_loglikelihood (const gsl_vector *coef, void *parms)
+{
+ size_t i;
+ struct pspp_arma_state *s;
+ double r = 0.0;
+ double tmp;
+ gsl_vector_view x;
+
+ s = (struct pspp_arma_state *) parms;
+ coef_to_arma (coef, s->arma);
+ i = (s->arma->ar_order > s->arma->ma_order) ? s->arma->ar_order :
s->arma->ma_order;
+ x = get_subvector (s->data, 0, i);
+ s->in = pspp_innovations (&x.vector, NULL, NULL);
+ get_residuals (s, &x.vector);
+ pspp_acf_vector (s);
+ for (i = 0; i < s->residuals->size; i++)
+ {
+ tmp = gsl_vector_get (s->scale, i);
+ assert (tmp > 0.0);
+ r += log (tmp);
+ }
+ r /= s->residuals->size;
+ r += log (s->arma->mse);
+ return r;
+}
+
+void
+pspp_arma_d_loglikelihood (const gsl_vector *coef, void *parms, gsl_vector *g)
+{
+ size_t m;
+ size_t i;
+ size_t j;
+ double tmp;
+ gsl_vector_view x;
+ gsl_vector_view resid;
+ gsl_vector_view data;
+ gsl_vector_view scale;
+ struct pspp_arma_state *s;
+
+ s = (struct pspp_arma_state *) parms;
+ coef_to_arma (coef, s->arma);
+ m = (s->arma->ar_order > s->arma->ma_order) ? s->arma->ar_order :
s->arma->ma_order;
+ x = get_subvector (s->data, 0, m);
+ s->in = pspp_innovations (&x.vector, NULL, NULL);
+ get_residuals (s, &x.vector);
+ pspp_acf_vector (s);
+ for (i = 0; i < s->arma->ar_order; i++)
+ {
+ resid = gsl_vector_subvector (s->residuals, m, s->residuals->size - m);
+ scale = gsl_vector_subvector (s->scale, m, s->residuals->size - m);
+ data = gsl_vector_subvector (s->data, m - i - 1, s->residuals->size - m);
+ tmp = 0.0;
+ for (j = 0; j < resid.vector.size; j++)
+ {
+ tmp += gsl_vector_get (&resid.vector, j) * gsl_vector_get
(&data.vector, j) /
+ gsl_vector_get (&scale.vector, j);
+ }
+ gsl_vector_set (g, i, -2.0 * tmp);
+ }
+ for (i = 0; i < s->arma->ma_order; i++)
+ {
+ resid = gsl_vector_subvector (s->residuals, m, s->residuals->size - m);
+ scale = gsl_vector_subvector (s->scale, m, s->residuals->size - m);
+ data = gsl_vector_subvector (s->residuals, m - i - 1, s->residuals->size
- m);
+ tmp = 0.0;
+ for (j = 0; j < resid.vector.size; j++)
+ {
+ tmp += gsl_vector_get (&resid.vector, j) * gsl_vector_get
(&data.vector, j) /
+ gsl_vector_get (&scale.vector, j);
+ }
+ gsl_vector_set (g, i + s->arma->ar_order, -2.0 * tmp);
+ }
+}
+/*
+ Likelihood and its gradient together.
+ */
+void
+pspp_arma_l_dl (const gsl_vector *coef, void *parms, double *like, gsl_vector
*grad)
+{
+ *like = pspp_arma_loglikelihood (coef, parms);
+ pspp_arma_d_loglikelihood (coef, parms, grad);
+}
+/*
+ COEF stores the autoregressive coefficients, then the
+ moving average coefficients.
+ */
+static void
+coef_to_arma (const gsl_vector *coef, pspp_arma *a)
+{
+ size_t i;
+
+ for (i = 0; i < a->ar_order; i++)
+ {
+ gsl_vector_set (a->ar_coeff, i, gsl_vector_get (coef, i));
+ }
+ for (i = a->ar_order; i < coef->size; i++)
+ {
+ gsl_vector_set (a->ma_coeff, i - a->ar_order, gsl_vector_get (coef, i));
+ }
+}
+static void
+arma_to_coef (gsl_vector *coef, pspp_arma *a)
+{
+ size_t i;
+ for (i = 0; i < a->ar_order; i++)
+ {
+ gsl_vector_set (coef, i, gsl_vector_get (a->ar_coeff, i));
+ }
+ for (i = a->ar_order; i < coef->size; i++)
+ {
+ gsl_vector_set (coef, i, gsl_vector_get (a->ma_coeff, i - a->ar_order));
+ }
+}
+
+/*
+ This function assumes DATA is the reparameterized process. This
+ function will transform to give the residuals for the original
+ process. This function also computes and stores the mean squared
+ error.
+ */
+static
+void get_residuals (struct pspp_arma_state *as, const gsl_vector *data)
+{
+ gsl_vector *tmp;
+ pspp_arma *a;
+ double resid;
+ double pred;
+ size_t i;
+ size_t j;
+
+ assert (as != NULL);
+ a = as->arma;
+ assert (as->residuals != NULL);
+ tmp = pspp_innovations_residuals (as->in, data);
+ for (i = 0; i < tmp->size; i++)
+ {
+ resid = gsl_vector_get (tmp, i);
+ gsl_vector_set (as->residuals, i, resid);
+ }
+ for (i = tmp->size; i < as->residuals->size; i++)
+ {
+ pred = 0.0;
+ for (j = 0; j < a->ar_order; j++)
+ {
+ pred += arma_get_ar_est (a, j) * gsl_vector_get (as->data, i - j - 1);
+ }
+ for (j = 0; j < a->ma_order; j++)
+ {
+ pred += arma_get_ma_est (a, j) * gsl_vector_get (as->residuals, i - j
- 1);
+ }
+ resid = gsl_vector_get (as->data, i) - pred;
+ gsl_vector_set (as->residuals, i, resid);
+ }
+ get_scale (as);
+ get_mse (as);
+ gsl_vector_free (tmp);
+}
+static
+void initialize_minimizer (gsl_vector *data,
+ gsl_multimin_function_fdf *m, size_t ar_order,
+ size_t ma_order)
+{
+ size_t order;
+ gsl_vector_view x;
+ struct pspp_arma_state *as = NULL;
+ assert (m != NULL);
+
+ as = xmalloc (sizeof (*as));
+ assert (as != NULL);
+ as->arma = pspp_arma_alloc (data->size, ar_order, ma_order);
+ assert (as->arma != NULL);
+ order = as->arma->ar_order + as->arma->ma_order;
+ m->f = &pspp_arma_loglikelihood;
+ m->df = &pspp_arma_d_loglikelihood;
+ m->fdf = &pspp_arma_l_dl;
+ m->n = order;
+ as->data = data;
+ as->inf_ma_coeff = gsl_vector_calloc (data->size);
+ as->scale = gsl_vector_calloc (data->size);
+ get_initial_estimates (as, data);
+ pspp_acf_vector (as);
+ as->residuals = gsl_vector_calloc (data->size);
+ x = get_subvector (as->data, 0, order);
+ get_residuals (as, &x.vector);
+ m->params = as;
+}
+
+/*
+ Find the maximum likelihood estimates of an ARMA process. The
+ function uses the EM algorithm with a gradient-descent. The
+ gradient is computed via an ad nauseum application of the chain rule
+ to the reparameterization outlined in Brockwell and Davis.
+ */
+pspp_arma *
+pspp_arma_fit (gsl_vector *data, size_t ar_order, size_t ma_order)
+{
+ pspp_arma *result = NULL;
+ struct pspp_arma_state *as = NULL;
+ const gsl_multimin_fdfminimizer_type *t;
+ gsl_multimin_fdfminimizer *s = NULL;
+ gsl_multimin_function_fdf *m = NULL;
+ gsl_vector *coef = NULL;
+ gsl_vector *gradient = NULL;
+ size_t order;
+ double tol = 0.01;
+ double step_size;
+ double epsabs = 1e-2;
+ double m_stage_epsabs = 0.1;
+ int rc;
+
+ t = gsl_multimin_fdfminimizer_vector_bfgs;
+ m = xmalloc (sizeof (*m));
+ initialize_minimizer (data, m, ar_order, ma_order);
+ as = (struct pspp_arma_state *) m->params;
+ order = ar_order + ma_order;
+ coef = gsl_vector_calloc (order);
+ s = gsl_multimin_fdfminimizer_alloc (t, order);
+ arma_to_coef (coef, as->arma);
+
+ step_size = 0.1 * gsl_blas_dnrm2 ((const gsl_vector *) as->residuals);
+ rc = gsl_multimin_fdfminimizer_set (s, m, coef, step_size, tol);
+ assert (rc != GSL_EBADLEN);
+
+ gsl_multimin_fdfminimizer_iterate (s);
+ gradient = gsl_multimin_fdfminimizer_gradient (s);
+ /*
+ EM algorithm.
+ */
+ while (gsl_multimin_test_gradient (gradient, epsabs) ==
+ GSL_CONTINUE)
+ {
+ while (gsl_multimin_test_gradient (gradient, m_stage_epsabs)
+ == GSL_CONTINUE)
+ {
+ /* Minimization */
+ gsl_multimin_fdfminimizer_iterate (s);
+ gradient = gsl_multimin_fdfminimizer_gradient (s);
+ }
+#if 0
+ coef_to_arma ((const gsl_vector *) coef, as->arma);
+ /* Expectation */
+ get_residuals (as, data);
+ pspp_acf_vector (as);
+#endif
+ gsl_multimin_fdfminimizer_restart (s);
+ }
+ gsl_multimin_fdfminimizer_free (s);
+ gsl_vector_free (as->residuals);
+ gsl_vector_free (as->inf_ma_coeff);
+ gsl_vector_free (as->scale);
+ gsl_vector_free (coef);
+ result = as->arma;
+ free (as);
+ free (m);
+ return result;
+}
+bool
+pspp_arma_free (pspp_arma *a)
+{
+ if (a != NULL)
+ {
+ if (a->ar_coeff != NULL)
+ {
+ gsl_vector_free (a->ar_coeff);
+ }
+ if (a->ma_coeff != NULL)
+ {
+ gsl_vector_free (a->ma_coeff);
+ }
+ if (a->acf != NULL)
+ {
+ gsl_vector_free (a->acf);
+ }
+ if (a->cov != NULL)
+ {
+ gsl_matrix_free (a->cov);
+ }
+ free (a);
+ return true;
+ }
+ return false;
+}
+
Index: dacf.c
===================================================================
RCS file: dacf.c
diff -N dacf.c
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ dacf.c 24 Jan 2007 17:28:35 -0000 1.1
@@ -0,0 +1,603 @@
+/*
+ src/math/ts/dacf.c
+
+ Copyright (C) 2007 Free Software Foundation, Inc.
+
+ 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.
+ */
+/*
+ Compute the derivative of autocovariance function of an ARMA (p,q)
+ process with respect to the autoregressive and moving average
+ parameter estimates. The process must be causal and invertible.
+ */
+#include <assert.h>
+#include <src/math/ts/acf.h>
+
+static double
+get_inf_ma_coef (const struct pspp_arma_state *x, int i)
+{
+ if (i < 0)
+ {
+ return 0.0;
+ }
+ if (i == 0)
+ {
+ return 1.0;
+ }
+ if (i < x->arma->ar_order)
+ {
+ return gsl_vector_get (x->inf_ma_coeff, (size_t) i);
+ }
+ return 0.0;
+}
+/*
+ Get the derivative of the infinite order moving average
+ coefficient with respect to the given autoregressive
+ coefficient.
+ */
+void
+set_d_psi_d_ar (struct pspp_arma_state *s)
+{
+ size_t i;
+ size_t phi;
+ size_t psi;
+ double deriv;
+
+ assert (s->d_inf_ma_d_ar != NULL);
+ assert (s->d_inf_ma_d_ar->size2 >= s->arma->ar_order);
+ gsl_matrix_set_zero (s->d_inf_ma_d_ar);
+ for (psi = 1; psi < s->d_inf_ma_d_ar->size1; psi++)
+ {
+ for (phi = 0; phi < s->d_inf_ma_d_ar->size2; phi++)
+ {
+ deriv = (psi >= phi + 1) ? get_inf_ma_coef (s, psi - phi - 1) : 0;
+ for (i = 0; i < psi - 1; i++)
+ {
+ deriv += arma_get_ar_est (s->arma, i) *
+ gsl_matrix_get (s->d_inf_ma_d_ar, psi - i - 1, i);
+ }
+ gsl_matrix_set (s->d_inf_ma_d_ar, psi, phi, deriv);
+ }
+ }
+}
+static double
+get_d_psi_d_ar (const struct pspp_arma_state *s, size_t psi, size_t phi)
+{
+ if ((psi >= s->d_inf_ma_d_ar->size1) || (phi >= s->d_inf_ma_d_ar->size2))
+ {
+ return 0.0;
+ }
+ return gsl_matrix_get (s->d_inf_ma_d_ar, psi, phi);
+}
+
+/*
+ Get the partial derivatives of the infinite-order MA process
+ with respect to the moving average parameters.
+ */
+void
+set_d_psi_d_ma (struct pspp_arma_state *s)
+{
+ size_t i;
+ size_t theta;
+ size_t psi;
+ double deriv;
+
+ assert (s->d_inf_ma_d_ar != NULL);
+ assert (s->d_inf_ma_d_ar->size2 >= s->arma->ar_order);
+ gsl_matrix_set_zero (s->d_inf_ma_d_ar);
+ for (psi = 1; psi < s->d_inf_ma_d_ar->size1; psi++)
+ {
+ for (theta = 0; theta < s->d_inf_ma_d_ar->size2; theta++)
+ {
+ deriv = (psi == theta + 1) ? 1.0 : 0.0;
+ for (i = 0; i < psi - 1; i++)
+ {
+ deriv += arma_get_ar_est (s->arma, i) *
+ gsl_matrix_get (s->d_inf_ma_d_ar, psi - i - 1, i);
+ }
+ gsl_matrix_set (s->d_inf_ma_d_ar, psi, theta, deriv);
+ }
+ }
+
+}
+/*
+ Retrieve the partial derivative of the infinite-order MA
+ coefficient with respect to the MA parameter THETA.
+ */
+static double
+get_d_psi_d_ma (const struct pspp_arma_state *s, int psi_coef, int theta)
+{
+ if ((theta < s->arma->ma_order) && (psi_coef >= 0) &&
+ (theta < s->arma->ma_order) && (theta >= 0))
+ {
+ return gsl_matrix_get (s->d_inf_ma_d_ma, psi_coef, theta);
+ }
+ return 0.0;
+}
+static double
+get_vector_element (const struct pspp_arma_state *s, int i, int phi)
+{
+ int j;
+ double result = 0.0;
+
+ for (j = i; j < s->arma->ma_order; j++)
+ {
+ result += arma_get_ma_est (s->arma, j) * get_d_psi_d_ma (s, j-i, phi);
+ }
+
+ return result;
+}
+static void
+get_vector_elements_ar (gsl_vector *b, const struct pspp_arma_state *x, size_t
phi)
+{
+ size_t i;
+ size_t j;
+ int k;
+ double tmp;
+ double element;
+
+ k = (int) phi + 1;
+ for (i = 0; i < x->arma->ar_order + 1; i++)
+ {
+ tmp = 0.0;
+ for (j = 0; j < x->arma->ma_order; j++)
+ {
+ tmp += arma_get_ma_est (x->arma, j + i) * get_d_psi_d_ar (x, j + 1,
phi);
+ }
+ element = pspp_arma_acf (x->arma, (size_t) abs (k)) + x->arma->mse * tmp;
+ gsl_vector_set (b, i, element);
+ k--;
+ }
+}
+static void
+get_vector_elements_ma (gsl_vector *b, const struct pspp_arma_state *x, size_t
ma)
+{
+ size_t i;
+ size_t j;
+ size_t k;
+ double tmp;
+ double element;
+
+ k = (int) ma + 1;
+ for (i = 0; i < x->arma->ar_order + 1; i++)
+ {
+ tmp = get_inf_ma_coef (x, k);
+ for (j = 0; j < x->arma->ma_order; j++)
+ {
+ tmp += arma_get_ma_est (x->arma, j + i) * get_d_psi_d_ma (x, j, ma);
+ }
+ element = x->arma->mse * tmp;
+ gsl_vector_set (b, i, element);
+ k--;
+ }
+}
+/*
+ Generic partial derivative for a parameter.
+ */
+static double
+d_acf_d_param (pspp_arma *x, gsl_vector *acf, size_t i, int lag)
+{
+ double result = 0.0;
+ size_t j;
+
+ while (i < (size_t) lag)
+ {
+ result = get_next_acf (x, acf);/*PROBLEM HERE, FIX THIS CALL, MAYBE USE
A DIFFERENT FUNCTION*/
+ for (j = 1; j < acf->size; j++)
+ {
+ gsl_vector_set (acf, (j-1), gsl_vector_get (acf, j));
+ }
+ gsl_vector_set (acf, acf->size - 1, result);
+ i++;
+ }
+ return result;
+}
+
+/*
+ Find the gradient of the autocovariance function with respect
+ to the autoregressive coefficients.
+
+ The acf is found by first solving A \gamma = b, where
+ A is an ar->order-by-ar->order matrix, \gamma is a column
+ vector with entry i equal to the acf at lag i, and b is
+ column vector containing a convolution of ar->ma_coef and
+ the infinite-order moving average coefficients.
+ */
+void
+pspp_d_acf_d_ar (struct pspp_arma_state *s)
+{
+ gsl_vector *result = NULL;
+ gsl_vector *b = NULL;
+ gsl_matrix *A = NULL;
+ gsl_vector *tau = NULL;
+ pspp_arma *x = NULL;
+ size_t i;
+
+ assert (s != NULL);
+ x = s->arma;
+ assert (x != NULL);
+ A = gsl_matrix_calloc (x->ar_order + 1, x->ar_order + 1);
+ result = gsl_vector_calloc (A->size1);
+ tau = gsl_vector_alloc (A->size1);
+ get_matrix_elements (A, x);
+ gsl_linalg_QR_decomp (A, tau);
+ b = gsl_vector_alloc (A->size1);
+ for (i = 0; i < x->ar_order; i++)
+ {
+ get_vector_elements_ar (b, (const struct pspp_arma_state *) s, i);
+ gsl_linalg_QR_solve ((const gsl_matrix *) A, (const gsl_vector *) tau,
+ (const gsl_vector *) b, result);
+ gsl_matrix_set_row (s->d_acf_d_phi, i, result);
+ }
+ gsl_matrix_free (A);
+ gsl_vector_free (tau);
+ gsl_vector_free (b);
+ gsl_vector_free (result);
+}
+static
+double get_d_sum_d_phi (struct pspp_arma_state *s, int lag, int coef)
+{
+ size_t i;
+ double result = 0.0;
+ pspp_arma *x;
+
+ if (lag >= s->arma->ar_order && lag >= s->arma->ma_order)
+ {
+ return 0.0;
+ }
+ x = s->arma;
+ for (i = 0; i < x->ma_order; i++)
+ {
+ result += arma_get_ma_est (x, i) *
+ get_d_psi_d_ma ((const struct pspp_arma_state *) s, i - lag, coef);
+ }
+ result *= x->mse;
+ return result;
+}
+
+static
+double get_d_sum_d_ma (struct pspp_arma_state *s, int lag, int coef)
+{
+ size_t i;
+ double result = 0.0;
+ pspp_arma *x;
+
+ if (lag >= s->arma->ar_order && lag >= s->arma->ma_order)
+ {
+ return 0.0;
+ }
+ x = s->arma;
+ for (i = 0; i < x->ma_order; i++)
+ {
+ result += arma_get_ma_est (x, i) *
+ get_d_psi_d_ma ((const struct pspp_arma_state *) s, i - lag, coef);
+ }
+ result += get_inf_ma_coef (s, coef - lag);
+ result *= x->mse;
+ return result;
+}
+static
+double get_d_sum_d_ar (struct pspp_arma_state *s, int lag, int coef)
+{
+ size_t i;
+ double result = 0.0;
+ pspp_arma *x;
+
+ if (lag >= s->arma->ar_order && lag >= s->arma->ma_order)
+ {
+ return 0.0;
+ }
+ x = s->arma;
+ for (i = 0; i < x->ma_order; i++)
+ {
+ result += arma_get_ma_est (x, i) *
+ get_d_psi_d_ar ((const struct pspp_arma_state *) s, i - lag, coef);
+ }
+ result += get_inf_ma_coef (s, coef - lag);
+ result *= x->mse;
+ return result;
+}
+static
+double get_next_d_acf_d_ar (struct pspp_arma_state *s, gsl_vector *x, int lag,
int coef)
+{
+ size_t i;
+ double result = 0.0;
+
+ for (i = 0; i < x->size; i++)
+ {
+ result += gsl_vector_get (x, i);
+ }
+ result += get_d_sum_d_ar (s, lag, coef) +
+ pspp_arma_acf (s->arma, (size_t) abs (lag - coef));
+
+ return result;
+}
+/*
+ Retrieve a partial derivative of the autocovariance
+ function.
+ */
+static double
+get_d_acf_d_ar (struct pspp_arma_state *s, int lag, int coef)
+{
+ double result = 0.0;
+ gsl_vector *acf;
+ pspp_arma *x;
+ size_t i;
+ size_t j;
+
+ x = s->arma;
+ lag = abs (lag);
+ if (coef >= s->d_acf_d_phi->size1)
+ {
+ return 0.0;
+ }
+ if (lag <= s->d_acf_d_phi->size2)
+ {
+ result = gsl_matrix_get (s->d_acf_d_phi, (size_t) coef, (size_t) lag);
+ }
+ else
+ {
+ /*
+ Store the next ar_order of the autocovariance
+ function.
+ */
+ acf = gsl_vector_calloc (s->d_acf_d_phi->size2 - 1);
+ for (i = 0; i < acf->size; i++)
+ {
+ gsl_vector_set (acf, i, gsl_matrix_get (s->d_acf_d_phi, coef, i + 1));
+ }
+
+ for (i = s->d_acf_d_phi->size2 + 1; i <= lag; i++)
+ {
+ for (j = 0; j < acf->size; j++)
+ {
+ gsl_vector_set (acf, j, gsl_vector_get (acf, j + 1));
+ }
+ result = get_next_d_acf_d_ar (s, acf, i, coef);
+ gsl_vector_set (acf, acf->size - 1, result);
+ }
+ gsl_vector_free (acf);
+ }
+ return result;
+}
+
+/*
+ Find the gradient of the autocovariance function with respect
+ to the moving average coefficients.
+
+ The acf is found by first solving A \gamma = b, where
+ A is an ma->order-by-ma->order matrix, \gamma is a column
+ vector with entry i equal to the acf at lag i, and b is
+ column vector containing a convolution of ar->ma_coef and
+ the infinite-order moving average coefficients.
+ */
+void
+pspp_d_acf_d_ma (struct pspp_arma_state *s)
+{
+ size_t i;
+ gsl_vector *result = NULL;
+ gsl_vector *b = NULL;
+ gsl_matrix *A = NULL;
+ gsl_vector *tau = NULL;
+ pspp_arma *x;
+
+ assert (s != NULL);
+ x = s->arma;
+ assert (x != NULL);
+ A = gsl_matrix_calloc (x->ar_order + 1, x->ar_order + 1);
+ result = gsl_vector_calloc (A->size1);
+ tau = gsl_vector_alloc (A->size1);
+ get_matrix_elements (A, x);
+ gsl_linalg_QR_decomp (A, tau);
+ b = gsl_vector_alloc (A->size1);
+ for (i = 0; i < x->ma_order; i++)
+ {
+ get_vector_elements_ma (b, (const struct pspp_arma_state *) s, i);
+ gsl_linalg_QR_solve ((const gsl_matrix *) A, (const gsl_vector *) tau,
+ (const gsl_vector *) b, result);
+ gsl_matrix_set_row (s->d_acf_d_ma, i, result);
+ }
+ gsl_matrix_free (A);
+ gsl_vector_free (tau);
+ gsl_vector_free (b);
+ gsl_vector_free (result);
+}
+static
+double get_next_d_acf_d_ma (struct pspp_arma_state *s, gsl_vector *x, int lag,
int coef)
+{
+ size_t i;
+ double result = 0.0;
+
+ for (i = 0; i < x->size; i++)
+ {
+ result += gsl_vector_get (x, i);
+ }
+ result += get_d_sum_d_ma (s, lag, coef);
+
+ return result;
+}
+/*
+ Retrieve a partial derivative of the autocovariance
+ function with respect to the moving average coefficients.
+ */
+static double
+get_d_acf_d_ma (struct pspp_arma_state *s, int lag, int coef)
+{
+ double result = 0.0;
+ gsl_vector *acf;
+ pspp_arma *x;
+ size_t i;
+ size_t j;
+
+ x = s->arma;
+ lag = abs (lag);
+ if (lag < 0)
+ {
+ return 0.0;
+ }
+ if (lag <= x->ma_order)
+ {
+ result = gsl_matrix_get (s->d_acf_d_ma, (size_t) coef, (size_t) lag);
+ }
+ else
+ {
+ /*
+ Store the next ma_order of the autocovariance
+ function.
+ */
+ acf = gsl_vector_calloc (s->d_acf_d_ma->size2 - 1);
+ for (i = 0; i < acf->size; i++)
+ {
+ gsl_vector_set (acf, i, gsl_matrix_get (s->d_acf_d_ma, coef, i + 1));
+ }
+
+ for (i = s->d_acf_d_ma->size2 + 1; i <= lag; i++)
+ {
+ for (j = 0; j < acf->size; j++)
+ {
+ gsl_vector_set (acf, j, gsl_vector_get (acf, j + 1));
+ }
+ result = get_next_d_acf_d_ma (s, acf, i, coef);
+ gsl_vector_set (acf, acf->size - 1, result);
+ }
+ gsl_vector_free (acf);
+ }
+ return result;
+}
+/*
+ Partial derivative of the autocovariance function of
+ the reparameterized process with respect to the autoregressive
+ parameters.
+ */
+static double
+d_reparam_acf_d_ar_form1 (struct pspp_arma_state *s, int i, int j,
+ int coeff)
+{
+ double result;
+
+ result = get_d_acf_d_ar (s, i - j, coeff) / s->arma->mse;
+
+ return result;
+}
+static double
+d_reparam_acf_d_ar_form2 (struct pspp_arma_state *s, size_t i, size_t j,
+ int coeff)
+{
+ size_t k;
+ size_t n;
+ double result;
+
+ n = (i > j) ? (i - j) : (j - i);
+ k = (coeff > n) ? (coeff + 1 - n) : (n - coeff - 1);
+ result = get_d_acf_d_ar (s, n, coeff) - pspp_arma_acf (s->arma, k);
+ for (k = 0; k < s->arma->ar_order; k++)
+ {
+ n = (i > j) ? (i - j) : (j - i);
+ n = ((k + 1) > n) ? (k + 1 - n) : (n - k - 1);
+ result -= get_d_acf_d_ar (s, (int) n, coeff) *
+ arma_get_ar_est (s->arma, k);
+ }
+ result /= s->arma->mse;
+ return result;
+}
+
+double
+d_reparam_acf_d_ar (struct pspp_arma_state *s, int i,
+ int j, int coef)
+{
+ pspp_arma *a;
+ double result = 0.0;
+ size_t m;
+ size_t minij;
+ size_t maxij;
+
+ a = s->arma;
+ m = (a->ar_order > a->ma_order) ? a->ar_order : a->ma_order;
+ minij = (i < j) ? i : j;
+ maxij = (i < j) ? j : i;
+ if ((i <= m) && (j <= m))
+ {
+ result = d_reparam_acf_d_ar_form1 (s, i, j, coef);
+ }
+ else if ((minij <= m) && (m < maxij) && (maxij <= 2 * m))
+ {
+ result = d_reparam_acf_d_ar_form2 (s, i, j, coef);
+ }
+
+ return result;
+}
+/*
+ Partial derivative of the autocovariance function of
+ the reparameterized process with respect to the moving average
+ parameters.
+ */
+static double
+d_reparam_acf_d_ma_form1 (struct pspp_arma_state *s, int i, int j,
+ int coeff)
+{
+ double result;
+
+ result = get_d_acf_d_ma (s, i - j, coeff) / s->arma->mse;
+
+ return result;
+}
+static double
+d_reparam_acf_d_ma_form2 (struct pspp_arma_state *s, size_t i, size_t j,
+ int coeff)
+{
+ int k;
+ double result;
+
+ result = get_d_acf_d_ma (s, i - j, coeff);
+ for (k = 0; k < s->arma->ma_order; k++)
+ {
+ result -= get_d_acf_d_ma (s, k - abs (i - j), coeff) *
+ arma_get_ma_est (s->arma, k);
+ }
+ result /= s->arma->mse;
+ return result;
+}
+
+double
+d_reparam_acf_d_ma (struct pspp_arma_state *s, int i,
+ int j, int coef)
+{
+ pspp_arma *a;
+ double result = 0.0;
+ size_t m;
+ size_t minij;
+ size_t maxij;
+
+ a = s->arma;
+ m = (a->ma_order > a->ma_order) ? a->ar_order : a->ma_order;
+ minij = (i < j) ? i : j;
+ maxij = (i < j) ? j : i;
+ if ((i <= m) && (j <= m))
+ {
+ result = d_reparam_acf_d_ma_form1 (s, i, j, coef);
+ }
+ else if ((minij <= m) && (m < maxij) && (maxij <= 2 * m))
+ {
+ result = d_reparam_acf_d_ma_form2 (s, i, j, coef);
+ }
+ else if (minij > m)
+ {
+ result = arma_get_ma_est (s->arma, coef) +
+ arma_get_ma_est (s->arma, coef - abs (i - j));
+ }
+
+ return result;
+}
Index: innovations.c
===================================================================
RCS file: innovations.c
diff -N innovations.c
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ innovations.c 24 Jan 2007 17:28:35 -0000 1.1
@@ -0,0 +1,322 @@
+/*
+ src/math/ts/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 <gsl/gsl_matrix.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_linalg.h>
+#include <stdlib.h>
+#include <libpspp/compiler.h>
+#include <libpspp/alloc.h>
+#include <math/ts/innovations.h>
+#include <assert.h>
+
+static void
+get_mean (gsl_vector *data,
+ struct innovations_estimate *est)
+
+{
+ size_t i;
+ double d;
+ double tmp;
+
+ est->n_obs = 0.0;
+ est->mean = 0.0;
+ for (i = 0; i < data->size; i++)
+ {
+ tmp = gsl_vector_get (data, i);
+ if (!gsl_isnan (tmp))
+ {
+ est->n_obs += 1.0;
+ d = (tmp - est->mean) / est->n_obs;
+ est->mean += d;
+ }
+ }
+}
+/*
+ Assumes the mean of X is 0.
+ */
+static double
+sample_acf (const gsl_vector *x, size_t lag)
+{
+ size_t i;
+ size_t j;
+ double xi;
+ double xj;
+ double result = 0.0;
+
+ assert (lag <= x->size);
+ for (i = 0; i < x->size - lag; i++)
+ {
+ j = i + lag;
+ xi = gsl_vector_get ((gsl_vector *) x, i);
+ xj = gsl_vector_get ((gsl_vector *) x, j);
+ result += xi * xj;
+ }
+ result /= (double) x->size;
+
+ return result;
+}
+/*
+ Given two size_t's, get the lag.
+ */
+static size_t
+get_lag (size_t x, size_t y)
+{
+ return ((x >= y) ? (x - y):(y - x));
+}
+
+/*
+ If ACF is not null, it will be used to compute the autocovariance at
+ lags 0,..., max_lag. The function ACF should accept ACF_ARG
+ and two subscripts, i and j, and return a double equal to the
+ autocovariance at (i, j).
+ */
+static int
+get_covariance (gsl_vector *data,
+ struct innovations_estimate *est,
+ double (*acf) (size_t, size_t, void *),
+ void *acf_arg)
+{
+ size_t i;
+ size_t j;
+ size_t k;
+ int rc = 1;
+ double tmp;
+
+ assert (data != NULL);
+ assert (est != NULL);
+
+ for (i = 0; i < est->cov->size1; i++)
+ {
+ for (j = 0; j <= i; j++)
+ {
+ /*
+ The entire covariance matrix is filled in here (as
+ opposed to only the upper triangle) so we can be stupid
+ later, and forget which part of the matrix has our data.
+ This isn't efficient, but it is more robust than
+ using only the upper triangle, since anyone who uses
+ the covariance matrix doesn't need to know how the data
+ are stored.
+ */
+ if (acf == NULL)
+ {
+ k = get_lag (i, j);
+ gsl_matrix_set (est->cov, i, j, sample_acf (data, k));
+ }
+ else
+ {
+ tmp = (*acf) (i, j, acf_arg);/*PROBLEM HERE: GIVES A NAN*/
+ assert (!gsl_isnan ((const double) tmp));
+ gsl_matrix_set (est->cov, i, j, tmp);
+ }
+ gsl_matrix_set (est->cov, j, i, gsl_matrix_get (est->cov, i, j));
+ }
+ assert (gsl_matrix_get (est->cov, i, i) > 0.0);
+ }
+ return rc;
+}
+
+static void
+get_coef (struct innovations_estimate *est)
+{
+ int rc;
+ size_t i;
+ size_t j;
+ double tmp;
+ double y;
+ /*
+ The coefficient matrix is stored with a '1'
+ on the main diagonal, so that the vector of
+ predicted values X_hat has this form:
+
+ X_hat = (I - inverse (EST->COEFF)) X
+
+ where X is the vector of input data and I is the
+ EST->N_OBS-dimensional square identity matrix.
+ */
+ gsl_matrix_memcpy (est->coeff, est->cov);
+ rc = gsl_linalg_cholesky_decomp (est->coeff);
+ assert (rc != GSL_EDOM);
+ for (i = 0; i < est->scale->size; i++)
+ {
+ tmp = gsl_matrix_get (est->coeff, i, i);
+ gsl_vector_set (est->scale, i, tmp * tmp);
+ /*
+ Scale the lower triangle and set the upper
+ triangle to 0.
+ */
+ assert (i < est->coeff->size1);
+ for (j = 0; j <= i; j++)
+ {
+ assert (j < est->coeff->size2);
+ y = gsl_matrix_get (est->coeff, i, j) / fabs (tmp);
+ gsl_matrix_set (est->coeff, i, j, y);
+ }
+ for (j = i + 1; j < est->coeff->size2; j++)
+ {
+ gsl_matrix_set (est->coeff, i, j, 0.0);
+ }
+ }
+}
+
+static void
+innovations_struct_init (struct innovations_estimate *est, size_t n_rows)
+{
+ est->mean = 0.0;
+ /*
+ COV[i][j] holds the covariance between X_i and X_j.
+ */
+ est->cov = gsl_matrix_calloc (n_rows, n_rows);
+ est->scale = gsl_vector_alloc (n_rows);
+ est->coeff = gsl_matrix_calloc (n_rows, n_rows); /* No intercept. */
+}
+/*
+ The mean is subtracted from the original data before computing the
+ coefficients. The mean is NOT added back, so if you want to predict
+ a new value, you must add the mean to X_hat[m] to get the correct
+ value.
+ */
+static void
+subtract_mean (gsl_vector *m, struct innovations_estimate *est)
+{
+ size_t i;
+ double tmp;
+
+ for (i = 0; i < m->size; i++)
+ {
+ tmp = gsl_vector_get ((gsl_vector *) m, i) - est->mean;
+ gsl_vector_set (m, i, tmp);
+ }
+}
+/*
+ Use the innovations algorithm to compute projections of
+ observations onto the span of the past observations.
+ If ACF is not null, it will be used to compute covariances
+ with the data stored in ACF_ARG.
+ If ACF is null, the sample autocovariance will be used.
+ */
+struct innovations_estimate *
+pspp_innovations (const gsl_vector *x, double (*acf) (size_t, size_t, void *),
+ void *acf_arg)
+{
+ struct innovations_estimate *est;
+ gsl_vector *data;
+
+ data = gsl_vector_alloc (x->size);
+ gsl_vector_memcpy (data, x);
+ est = xmalloc (sizeof (*est));
+ innovations_struct_init (est, data->size);
+ get_mean (data, est);
+ subtract_mean (data, est);
+ get_covariance (data, est, acf, acf_arg);
+ get_coef (est);
+ gsl_vector_free (data);
+ return est;
+}
+
+void
+pspp_innovations_free (struct innovations_estimate *est)
+{
+ assert (est != NULL);
+ gsl_matrix_free (est->coeff);
+ gsl_matrix_free (est->cov);
+ gsl_vector_free (est->scale);
+ free (est);
+}
+
+/*
+ One-step prediction. Element i of RESULT is the predicted
+ value for element i of X.
+ */
+gsl_vector *
+pspp_innovations_predicted_values (struct innovations_estimate *est,
+ gsl_vector *x)
+{
+ size_t i;
+ size_t j;
+ double tmp;
+ gsl_matrix *coef_m_i;
+ gsl_vector *result;
+
+ result = gsl_vector_alloc (x->size);
+ gsl_vector_memcpy (result, x);
+ coef_m_i = gsl_matrix_alloc (est->coeff->size1, est->coeff->size2);
+ for (i = 1; i < coef_m_i->size1; i++)
+ {
+ for (j = 0; j < i; j++)
+ {
+ gsl_matrix_set (coef_m_i, i, j, gsl_matrix_get (est->coeff, i, j));
+ }
+ }
+ gsl_blas_dtrmv (CblasLower, CblasNoTrans, CblasNonUnit, coef_m_i, result);
+ gsl_vector_set (result, 0, 0.0);
+ for (i = 1; i < result->size; i++)
+ {
+ tmp = gsl_vector_get (result, i);
+ for (j = 0; j < i; j++)
+ {
+ tmp -= gsl_matrix_get (est->coeff, i, j) * gsl_vector_get (result, j);
+ }
+ gsl_vector_set (result, i, tmp);
+ }
+ gsl_matrix_free (coef_m_i);
+ return result;
+}
+/*
+ Residuals from one-step prediction. Element i of RESULT is the
+ residual value for element i of X.
+ */
+gsl_vector *
+pspp_innovations_residuals (struct innovations_estimate *est,
+ const gsl_vector *x)
+{
+ gsl_vector *resid;
+ gsl_vector_view x1;
+
+ /*
+ If X is too long, just use the first EST->N_OBS entries.
+ */
+ if (x->size >= est->n_obs)
+ {
+ x1 = gsl_vector_subvector (x, 0, est->n_obs);
+ resid = pspp_innovations_predicted_values (est, &x1.vector);
+ }
+ else
+ {
+ resid = pspp_innovations_predicted_values (est, x);
+ }
+
+ return resid;
+}
+
Index: likelihood.c
===================================================================
RCS file: likelihood.c
diff -N likelihood.c
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ likelihood.c 24 Jan 2007 17:28:35 -0000 1.1
@@ -0,0 +1,428 @@
+/*
+ src/math/ts/likelihood.c
+
+ Copyright (C) 2006 Free Software Foundation, Inc.
+
+ 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.
+ */
+
+/*
+ Compute the likelihood, its gradient, and other related values for
+ an ARMA process.
+ */
+#include <assert.h>
+#include <gsl/gsl_math.h>
+#include <stdlib.h>
+#include <libpspp/alloc.h>
+#include <libpspp/compiler.h>
+#include <src/math/ts/arma-optimizer.h>
+#include <src/math/ts/innovations.h>
+#include <src/math/ts/acf.h>
+/*
+ Each gsl_matrix in a d_innovations contains the
+ partial derivatives of the innovations coefficients
+ with respect to a single parameter.
+ */
+struct d_innovations
+{
+ gsl_matrix **d_ar;
+ gsl_matrix **d_ma;
+};
+
+static double
+get_scale (struct innovations_estimate *in, size_t i)
+{
+ return (i < in->n_obs) ? gsl_vector_get (in->scale, i) : gsl_matrix_get
(in->cov, 0, 0);
+}
+
+/*
+ RESID should be the residuals for the reparameterized
+ process.
+ */
+static double
+log_unnormalized_likelihood (const gsl_vector *resid,
+ struct innovations_estimate *in)
+{
+ double result = 0.0;
+ double r1 = 0.0;
+ double r2 = 0.0;
+ double tmp;
+ double tmp2;
+ size_t i;
+
+ for (i = 0; i < resid->size - 1; i++)
+ {
+ tmp = gsl_vector_get ((gsl_vector *) resid, i);
+ tmp2 = get_scale (in, i);
+ r1 += tmp * tmp / tmp2;
+ r2 += log (tmp2);
+ }
+ result = log (r1 / resid->size) + r2 / resid->size;
+ assert (!isnan (result));
+
+ return result;
+}
+/*
+ Reduced likelihood (unnormalized and computed with the maximum likelihood
+ estimate of the variance).
+ */
+double
+pspp_arma_log_likelihood (const gsl_vector *resid,
+ struct innovations_estimate *in)
+{
+ double result = 0.0;
+
+ result = log_unnormalized_likelihood (resid, in);
+
+ return result;
+}
+static double
+theta_summand_1 (gsl_matrix *deriv, gsl_matrix *coef, double r,
+ gsl_matrix *d_r_d_ar,
+ size_t n, size_t k, size_t j, size_t phi)
+{
+ double result;
+
+ result = r * (gsl_matrix_get (deriv, k, k - j) *
+ gsl_matrix_get (coef, n, n - j) +
+ gsl_matrix_get (deriv, n, n - j) *
+ gsl_matrix_get (coef, k, k - j)) +
+ gsl_matrix_get (d_r_d_ar, j, phi) * gsl_matrix_get (coef, k, k - j) *
+ gsl_matrix_get (coef, n, n - j);
+ return result;
+}
+static double
+theta_summand_2 (gsl_matrix *coef, double r, size_t n, size_t k, size_t j)
+{
+ double result;
+
+ result = gsl_matrix_get (coef, k, k - j) * gsl_matrix_get (coef, n, n - j)
+ * r;
+
+ return result;
+}
+/*
+ Partial derivative of the innovations estimated coefficients with
+ respect to the autoregressive parameter. Make sure IN points to the
+ innovations estimate of the reparameterized process. D_REPARAM_ACF
+ should be one of D_REPARAM_ACF_D_AR or D_REPARAM_ACF_D_MA.
+*/
+static gsl_vector *
+set_d_innov (struct pspp_arma_state *s,
+ struct innovations_estimate *in,
+ gsl_matrix *d_theta,
+ gsl_matrix *d_r_d_ar, size_t n, size_t phi,
+ double (*d_reparam_acf) (struct pspp_arma_state*, int, int, int))
+{
+ size_t k;
+ size_t j;
+ double sum1;
+ double sum2;
+ double tmp;
+ double tmp2;
+ gsl_vector *result;
+
+ result = gsl_vector_calloc (in->coeff->size2);
+
+#if 0
+ sum1 = (*d_reparam_acf) (s, n + 1, 0, (int) phi) *
+ reparam_acf (1, 1, s->arma) -
+ reparam_acf (n + 1, 1, s->arma) *
+ (*d_reparam_acf) (s, 1, 1, (int) phi); /*CHECK THIS*/
+ tmp = get_scale (in, 0);
+ tmp = sum1 / (tmp * tmp);
+ gsl_vector_set (result, 0, tmp);
+#endif
+ for (k = 0; (k <= n) && (n < k + result->size); k++)
+ {
+ sum1 = 0.0;
+ sum2 = 0.0;
+ for (j = 0; j < k; j++)
+ {
+ sum1 += theta_summand_1 (d_theta, in->coeff, get_scale (in, j),
+ d_r_d_ar, n, k, j, phi);
+ sum2 += theta_summand_2 (in->coeff, get_scale (in, j), n, k, j);
+ }
+ tmp = (*d_reparam_acf) (s, n + 1, k, phi) - sum1;
+ tmp2 = get_scale (in, k);
+ tmp -= (reparam_acf (n + 1, k, s->arma) - sum2) *
+ gsl_matrix_get (d_r_d_ar, k, phi) / tmp2;
+ tmp /= tmp2;
+ gsl_vector_set (result, n - k, tmp);
+ }
+
+ return result;
+}
+/*
+ Get the derivative of the scale parameter for the innovations
+ estimates with respect to either the MA or the AR parameter.
+ D_REPARAM_ACF should be D_REPARAM_ACF_D_MA or D_REPARAM_ACF_D_AR.
+
+ Element J, COEF of D_R is the partial derivative of IN->SCALE[J] with
+ parameter COEF.
+ */
+static void
+set_d_r (struct pspp_arma_state *s, size_t coef, size_t j,
+ struct innovations_estimate *in, gsl_matrix *d_theta, gsl_matrix *d_r,
+ double (*d_reparam_acf) (struct pspp_arma_state *, int, int, int))
+{
+ double result = 0.0;
+ double tmp;
+ size_t k;
+
+ assert (j < d_r->size1);
+ result = (*d_reparam_acf) (s, (int) j + 1, (int) j + 1, coef);
+
+ for (k = 0; (k + 1) < j; k++)
+ {
+ tmp = gsl_matrix_get (in->coeff, j, j - k - 1);
+ result -= tmp *
+ (2.0 * gsl_matrix_get (d_theta, j, j - k - 1) *
+ get_scale (in, k) + tmp * gsl_matrix_get (d_r, k, coef));
+ }
+ gsl_matrix_set (d_r, j, coef, result);
+}
+
+static double
+d_S_d_ar (struct pspp_arma_state *s, struct innovations_estimate *in,
gsl_vector *d_r, size_t coef)
+{
+ double result = 0.0;
+ double resid;
+ double tmp;
+ size_t m;
+ size_t i;
+
+ m = (s->arma->ar_order > s->arma->ma_order) ? s->arma->ar_order :
s->arma->ma_order;
+ for (i = 0; i < m; i++)
+ {
+ resid = gsl_vector_get (s->residuals, i);
+ tmp = get_scale (in, i);
+ result -= resid * resid * gsl_vector_get (d_r, i) /
+ (tmp * tmp);
+ }
+ for (i = m; i < d_r->size; i++)
+ {
+ resid = gsl_vector_get (s->residuals, i);
+ tmp = get_scale (in, i);
+ result -= resid * resid * gsl_vector_get (d_r, i) /
+ (tmp * tmp);
+ result -= 2.0 * resid * gsl_vector_get (s->data, i - coef - 1) /
+ tmp;
+ }
+ return result;
+}
+static double
+d_S_d_ma (struct pspp_arma_state *s, struct innovations_estimate *in,
gsl_vector *d_r, size_t coef)
+{
+ double result = 0.0;
+ double resid;
+ double tmp;
+ size_t m;
+ size_t i;
+
+ m = (s->arma->ar_order > s->arma->ma_order) ? s->arma->ar_order :
s->arma->ma_order;
+ for (i = 0; i < m; i++)
+ {
+ resid = gsl_vector_get (s->residuals, i);
+ tmp = get_scale (in, i);
+ result -= resid * resid * gsl_vector_get (d_r, i) /
+ (tmp * tmp);
+ }
+ for (i = m; i < d_r->size; i++)
+ {
+ resid = gsl_vector_get (s->residuals, i);
+ tmp = get_scale (in, i);
+ result -= resid * resid * gsl_vector_get (d_r, i) /
+ (tmp * tmp);
+ result -= 2.0 * resid * gsl_vector_get (s->residuals, i - coef) /
+ tmp;
+ }
+ return result;
+}
+static double
+get_d_log_normalizer (gsl_vector *dr, struct innovations_estimate *in)
+{
+ double result = 0.0;
+ size_t i;
+
+ for (i = 0; i < dr->size; i++)
+ {
+ result += gsl_vector_get (dr, i) / get_scale (in, i);
+ }
+ result /= dr->size;
+ return result;
+}
+/*
+ Actual computation of the gradient. This function is
+ generic, meaning it will compute the gradient of the log
+ likelihood for either the AR parameters or the MA parameters,
+ depending on the function D_REPARAM_ACF.
+ */
+static void
+grad_ar (struct pspp_arma_state *s,
+ gsl_vector *result,
+ gsl_matrix **d_theta, size_t order,
+ double (*d_reparam_acf) (struct pspp_arma_state *, int, int, int))
+{
+ struct innovations_estimate *in;
+ gsl_matrix *d_r = NULL;
+ gsl_vector *theta_row = NULL;
+ gsl_vector_view d_ri;
+ double d_log_normalizer;
+ double partial_deriv;
+ double dS;
+ size_t coef;
+ size_t n;
+
+ in = s->in;
+ d_r = gsl_matrix_calloc (in->coeff->size1, order);
+ for (coef = 0; coef < order; coef++)
+ {
+ set_d_r (s, coef, 0, in, d_theta[coef], d_r, d_reparam_acf);
+ for (n = 0; n < d_theta[coef]->size1; n++)
+ {
+ theta_row = set_d_innov (s, in, d_theta[coef],
+ d_r, n, coef, d_reparam_acf);
+ gsl_matrix_set_row (d_theta[coef], n,
+ (const gsl_vector *) theta_row);
+ set_d_r (s, coef, n, in, d_theta[coef], d_r, d_reparam_acf);
+ }
+ }
+ for (coef = 0; coef < order; coef++)
+ {
+ d_ri = gsl_matrix_column (d_r, coef);/*CHECK THIS*/
+/* tmp = S (resid_sq, in); */
+/* assert (fabs (tmp) > GSL_DBL_EPSILON); */
+ dS = d_S_d_ar (s, in, &d_ri.vector, coef);
+ d_log_normalizer = get_d_log_normalizer (&d_ri.vector, in);
+ partial_deriv = dS / s->arma->mse + d_log_normalizer;
+ gsl_vector_set (result, coef, partial_deriv);
+ }
+ gsl_matrix_free (d_r);
+}
+static void
+grad_ma (struct pspp_arma_state *s,
+ gsl_vector *result,
+ gsl_matrix **d_theta, size_t order,
+ double (*d_reparam_acf) (struct pspp_arma_state *, int, int, int))
+{
+ struct innovations_estimate *in;
+ gsl_matrix *d_r = NULL;
+ gsl_vector *theta_row = NULL;
+ gsl_vector_view d_ri;
+ double d_log_normalizer;
+ double partial_deriv;
+ double dS;
+ size_t coef;
+ size_t n;
+
+ in = s->in;
+ d_r = gsl_matrix_calloc (in->coeff->size1, order);
+ for (coef = 0; coef < order; coef++)
+ {
+ set_d_r (s, coef, 0, in, d_theta[coef], d_r, d_reparam_acf);
+ for (n = 0; n < d_theta[coef]->size1; n++)
+ {
+ theta_row = set_d_innov (s, in, d_theta[coef],
+ d_r, n, coef, d_reparam_acf);
+ gsl_matrix_set_row (d_theta[coef], n,
+ (const gsl_vector *) theta_row);
+ set_d_r (s, coef, n, in, d_theta[coef], d_r, d_reparam_acf);
+ }
+ }
+ for (coef = 0; coef < order; coef++)
+ {
+ d_ri = gsl_matrix_column (d_r, coef);
+ dS = d_S_d_ma (s, in, &d_ri.vector, coef);
+ d_log_normalizer = get_d_log_normalizer (&d_ri.vector, in);
+ partial_deriv = dS / s->arma->mse + d_log_normalizer;
+ gsl_vector_set (result, coef, partial_deriv);
+ }
+ gsl_matrix_free (d_r);
+}
+/*
+ Make sure RESID contains the residuals of the reparameterized
+ process, not the original data. This function computes the exact
+ gradient of the log-likelihood by applying the chain many times to a
+ tangled mess of intermediate variables and functions. It overwrites
+ RESULT with the gradient.
+ */
+void
+pspp_arma_likelihood_gradient (struct pspp_arma_state *s,
+ gsl_vector *result)
+{
+ struct d_innovations d_theta;
+ struct innovations_estimate *in;
+ gsl_matrix *d_r;
+ size_t n;
+ gsl_vector_view d_ar;
+ gsl_vector_view d_ma;
+
+ assert (result != NULL);
+ /*
+ Update the matrices storing the gradient of the acf with respect to
+ the parameters.
+ */
+ set_d_psi_d_ar (s);
+ set_d_psi_d_ma (s);
+ pspp_d_acf_d_ma (s);
+ pspp_d_acf_d_ar (s);
+ in = s->in;
+ /*
+ Row i contains the partial derivative of IN->SCALE[J]
+ with respect to autoregressive coefficient i.
+ */
+ d_r = gsl_matrix_calloc (in->coeff->size2, in->coeff->size2);
+ d_theta.d_ar = xmalloc (sizeof (*d_theta.d_ar));
+ for (n = 0; n < s->arma->ar_order; n++)
+ {
+ /* Each matrix contains the partial derivative of the
+ innovation coefficients with respect to autoregressive
+ parameter n.
+ */
+ d_theta.d_ar[n] = gsl_matrix_calloc (in->coeff->size1, in->coeff->size2);
+ }
+ d_theta.d_ma = xmalloc (sizeof (*d_theta.d_ma));
+ for (n = 0; n < s->arma->ma_order; n++)
+ {
+ /* Each matrix contains the partial derivative of the
+ innovation coefficients with respect to moving average
+ parameter n.
+ */
+ d_theta.d_ma[n] = gsl_matrix_calloc (in->coeff->size1, in->coeff->size2);
+ }
+ /*
+ The gradient vector stores the partial derivatives with respect to
+ the autoregressive parameters first, then the partial derivative
+ with respect to the moving average parameters.
+ */
+ d_ar = gsl_vector_subvector (result, 0, s->arma->ar_order);
+ grad_ar (s, &d_ar.vector, d_theta.d_ar, s->arma->ar_order,
+ &d_reparam_acf_d_ar);
+ d_ma = gsl_vector_subvector (result, s->arma->ar_order, s->arma->ma_order);
+ grad_ma (s, &d_ma.vector, d_theta.d_ma, s->arma->ma_order,
+ &d_reparam_acf_d_ma);
+
+ gsl_matrix_free (d_r);
+ for (n = 0; n < s->arma->ar_order; n++)
+ {
+ gsl_matrix_free (d_theta.d_ar[n]);
+ }
+ for (n = 0; n < s->arma->ma_order; n++)
+ {
+ gsl_matrix_free (d_theta.d_ma[n]);
+ }
+}
+
Index: predict.c
===================================================================
RCS file: predict.c
diff -N predict.c
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ predict.c 24 Jan 2007 17:28:35 -0000 1.1
@@ -0,0 +1,21 @@
+/*
+ For an ARMA model, the prediction function can ignore the
+ VARIABLE argument. The argument is still in the function's
+ definition to keep the argument list of this prediction
+ function the same as that for other models.
+ */
+double
+pspp_arma_predict (const struct variable **var, const union value **val,
+ const void *arma_, int n)
+{
+ double result = 0.0;
+ pspp_arma *arma;
+ size_t i;
+ size_t j;
+
+ arma = (pspp_arma *) arma_;
+ result = val[0].f;
+ for (i = 0; i < arma->ar_order; i++)
+ {
+ result += val[i] * arma_get_ar_est (arma, i);
+ }
Index: tmp.c
===================================================================
RCS file: tmp.c
diff -N tmp.c
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ tmp.c 24 Jan 2007 17:28:35 -0000 1.1
@@ -0,0 +1,7 @@
+#include <stdio.h>
+int main()
+{
+ printf ("%lf\n", 10e+9);
+ return;
+}
+
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Pspp-cvs] experimental/src/math/ts acf.c arma.c dacf.c in...,
Jason H Stover <=