[Pspp-cvs] experimental/src/math/ts arma.c

[Jason H Stover]

CVSROOT:        /sources/pspp
Module name:    experimental
Changes by:     Jason H Stover <jstover>        07/04/06 21:03:34

Modified files:
        src/math/ts    : arma.c 

Log message:
        reimplementing reparam_acf

CVSWeb URLs:
http://cvs.savannah.gnu.org/viewcvs/experimental/src/math/ts/arma.c?cvsroot=pspp&r1=1.1&r2=1.2

Patches:
Index: arma.c
===================================================================
RCS file: /sources/pspp/experimental/src/math/ts/arma.c,v
retrieving revision 1.1
retrieving revision 1.2
diff -u -b -r1.1 -r1.2
--- arma.c      24 Jan 2007 17:28:35 -0000      1.1
+++ arma.c      6 Apr 2007 21:03:33 -0000       1.2
@@ -69,11 +69,172 @@
   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);
+  arma->acf = gsl_vector_calloc (n_obs);
   return arma;
 }
+/*
+  Reparameterized acf, as defined in section 5.3 of Brockwell and Davis.
+*/
+static double
+reparam_acf_form1 (size_t i, size_t j, pspp_arma *a, gsl_vector *acf)
+{
+  double result;
+  int lag;
+  lag = (i > j) ? (i - j) : (j - i);
+  result = gsl_vector_get (acf, lag) / a->mse;
+  return result;
+}
+static double
+reparam_acf_form2 (size_t i, size_t j, pspp_arma *a, gsl_vector *acf)
+{
+  double result;
+  double phi;
+  int lag;
+  int n;
+  size_t k;
+  
+  n =  (i > j) ? (i - j) : (j - i);
+  result = pspp_arma_acf (a, n);
+  for (k = 0; k < a->ar_order; k++)
+    {
+      lag = k - n;
+      phi = arma_get_ar_est (a, k);
+      result -= phi * gsl_vector_get (acf, lag);
+    }
+  result /= a->mse;
+  return result;
+}
+static double
+reparam_acf_form3 (size_t i, size_t j, pspp_arma *a)
+{
+  double result = 0.0;
+  double *v1;
+  double *v2;
+  int k;
+  int imj;
+
+
+  v1 = nmalloc (a->ma_order + 1, sizeof (*v1));
+  v2 = nmalloc (a->ma_order + 1, sizeof (*v2));
+  v1[0] = 1.0;
+  for (k = 1; k < a->ma_order + 1; k++)
+    {
+      v1[k] = arma_get_ma_est (a, (size_t) k);
+    }
+  if (i == j)
+    {
+      imj = 1;
+      v2[0] = 1.0;
+    }
+  else
+    {
+      imj = (i > j) ? (i - j) : (j - i);
+      imj--;
+    }
+  for (k = imj; k < a->ma_order + 1; k++)
+    {
+      v2[k] = arma_get_ma_est (a, (size_t) k);
+    }
+  for (k = 0; k < a->ma_order + 1; k++)
+    {
+      result += v1[k] * v2[k];
+    }
+  free (v1);
+  free (v2);
+  return result;
+}
 
 double
+reparam_acf (size_t i, size_t j, pspp_arma *a)
+{
+  pspp_arma *a;
+  double result = 0.0;
+  size_t m;
+  size_t minij;
+  size_t maxij;
+  
+  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 = reparam_acf_form1 (i, j, a);
+    }
+  else if ((minij <= m) && (m < maxij) && (maxij <= 2 * m))
+    {
+      result = reparam_acf_form2 (i, j, a);
+    }
+  else if (minij > m)
+    {
+      result = reparam_acf_form3 (i, j, a);
+    }
+  return result;
+}
+/*
+  Lose this function once you figure out how to make this work right.
+*/
+static void 
+get_inf_ma (struct pspp_arma_state *s)
+{
+  size_t i;
+  double tmp;
+  double tmp2;
+  gsl_vector_set (s->inf_ma_coeff, 0, 1.0);
+  tmp = gsl_vector_get (s->arma->ar_coeff, 0) + gsl_vector_get 
(s->arma->ma_coeff, 0);
+  gsl_vector_set (s->inf_ma_coeff, 1, tmp);
+  for (i = 2; i < s->data->size; i++)
+    {
+      tmp = gsl_vector_get (s->arma->ar_coeff, 0) * gsl_vector_get 
(s->inf_ma_coeff, i - 1);
+      tmp2 = gsl_vector_get (s->arma->ar_coeff, 1) * gsl_vector_get 
(s->inf_ma_coeff, i - 2);
+      gsl_vector_set (s->inf_ma_coeff, i, tmp+tmp2);
+    }
+}
+/*
+  Remove this function when you figure out how to make this all work.
+ */
+static void
+acf_vector (struct pspp_arma_state *s, gsl_vector *coef)
+{
+  gsl_matrix *m;
+  gsl_vector *v;
+  gsl_vector *tmp;
+  double x;
+  double phi1;
+  double phi2;
+  double theta;
+  int i;
+  size_t j;
+  /*BUSTED IN HERE. DOES NOT CORRECTLY SOLVE SYSTEM.*/
+  m = gsl_matrix_calloc (3,3);
+  v = gsl_vector_calloc (3);
+  tmp = gsl_vector_calloc (3);
+  gsl_matrix_set (m, 0,0, 1.0);
+  gsl_matrix_set (m, 0, 1, -1.0 * gsl_vector_get (coef, 0));
+  gsl_matrix_set (m, 1, 0, -1.0 * gsl_vector_get (coef, 0));
+  gsl_matrix_set (m, 2, 1, -1.0 * gsl_vector_get (coef, 0));
+  gsl_matrix_set (m, 0, 2, -1.0 * gsl_vector_get (coef, 1));
+  gsl_matrix_set (m, 2, 0, -1.0 * gsl_vector_get (coef, 1));
+  gsl_matrix_set (m, 1, 1, 1.0 - gsl_vector_get (coef, 0));
+  gsl_matrix_set (m, 2, 2, 1.0);
+  gsl_vector_set (v, 0, s->arma->mse * (1.0 + gsl_vector_get (coef, 2) * 
(gsl_vector_get (coef, 2) + gsl_vector_get (coef, 0))));
+  gsl_vector_set (v, 1, s->arma->mse * gsl_vector_get (coef, 2));
+  gsl_permutation *p = gsl_permutation_alloc (3);
+  gsl_linalg_LU_decomp (m, p, &i);
+  gsl_linalg_LU_solve (m, p, v, tmp);
+
+  for (j = 0; j < 3; j++)
+    {
+      gsl_vector_set (s->arma->acf, j, gsl_vector_get (tmp, j));
+    }
+
+  for (j = 3; j < s->arma->acf->size; j++)
+    {
+      x = gsl_vector_get (coef,0) * gsl_vector_get (s->arma->acf, j - 1) +
+       gsl_vector_get (coef, 1) * gsl_vector_get (s->arma->acf, j - 2);
+      gsl_vector_set (s->arma->acf, j, x);
+    }
+}
+double
 arma_get_ar_est (const pspp_arma *x, size_t i)
 {
   double result = 0.0;
@@ -124,6 +285,7 @@
       gsl_vector_set (as->scale, i, tmp);
     }
 }
+
 static void
 get_mse (struct pspp_arma_state *as)
 {
@@ -131,25 +293,28 @@
   double tmp;
   
   as->arma->mse = 0.0;
-  for (i = 0; i < as->residuals->size; i++)
+  for (i = 0; i < as->in->residuals->size; i++)
     {
-      tmp = gsl_vector_get (as->residuals, i);
+      tmp = gsl_vector_get (as->in->residuals, i);
       as->arma->mse += tmp * tmp / gsl_vector_get (as->scale, i);
     }
-  as->arma->mse /= as->residuals->size;
+  as->arma->mse /= as->in->residuals->size;
   as->arma->std = sqrt (as->arma->mse);
 }
+
 /*
-  Initial values are derived from the innovations estimates
-  via least-squares regression.
+  Initial values are derived 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;
+  double *pred;
+  double *resid;
+  double *d_data;
+  double c;
+  double ssq;
   gsl_vector_view y;
   gsl_vector *coef;
   gsl_matrix *obs;
@@ -161,43 +326,44 @@
   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);
+  order = a->ar_order;
+  obs = gsl_matrix_alloc (data->size - order, order);
   for (j = 0; j < obs->size1; j++)
     {
-      for (i = 0; i < a->ar_order; i++)
+      for (i = 0; i < 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);
+  y = gsl_vector_subvector (data, order, obs->size1);
   ws = gsl_multifit_linear_alloc (obs->size1, order);
   rc = gsl_multifit_linear (obs, &y.vector, coef, cov, &mse, ws);
+  a->mse = mse / (obs->size1 - order);
   assert (rc == GSL_SUCCESS);
-  coef_to_arma (coef, a);
+  gsl_vector_set (a->ar_coeff, 0, gsl_vector_get (coef, 0));
+  gsl_vector_set (a->ar_coeff, 1, gsl_vector_get (coef, 1)); 
   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);
+  pred = nmalloc (obs->size1 - 1, sizeof (*pred));
+  resid = nmalloc (obs->size1 - 1, sizeof (*resid));
+  d_data = nmalloc (obs->size1 - 1, sizeof (*d_data));
+  for (j = 0; j < (obs->size1 - 1); j++)
+    {
+      pred[j] = gsl_vector_get (coef, 0) * gsl_matrix_get (obs, j, 0) + 
+       gsl_vector_get (coef, 1) * gsl_matrix_get (obs, j, 1);
+      resid[j] = gsl_vector_get (&y.vector, j) - pred[j];
+      d_data[j] = gsl_vector_get (&y.vector, j + 1);
     }
-  a->mse /= resid->size;
+  gsl_fit_mul (resid, 1, d_data , 1, obs->size1, &c, &mse, &ssq);
   a->std = sqrt (a->mse);
-  gsl_vector_free (resid);
-  gsl_vector_free (pred);
+  gsl_vector_set (a->ma_coeff, 0, 0.0);
   gsl_vector_free (coef);
   gsl_matrix_free (cov);
+  gsl_matrix_free (obs);
+  free (resid);
+  free (pred);
+  free (d_data);
 }
 /*
   Get the smallest subvector of size at least ORDER
@@ -217,6 +383,7 @@
     }
   return result;
 }
+
 /*
   COEF should store the AR and MA coefficients, in that
   order. 
@@ -225,75 +392,155 @@
 pspp_arma_loglikelihood (const gsl_vector *coef, void *parms)
 {
   size_t i;
+  size_t j;
+  size_t lag;
+  int rc;
   struct pspp_arma_state *s;
-  double r = 0.0;
   double tmp;
-  gsl_vector_view x;
+  double result = 0.0;
+  double S = 0.0;
+  double log_normalizer = 0.0;
+  double mse;
+  double y;
+  gsl_matrix *gamma; /* covariance */
+  gsl_vector *scale; /* mean-squared errors */
+  gsl_vector *residuals;
 
   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++)
+
+  acf_vector (s, coef);
+  residuals = gsl_vector_calloc (s->data->size);
+  gamma = gsl_matrix_calloc (residuals->size, residuals->size);
+  scale = gsl_vector_calloc (residuals->size);
+  for (i = 0; i < residuals->size; i++)
     {
-      tmp = gsl_vector_get (s->scale, i);
-      assert (tmp > 0.0);
-      r += log (tmp);
+      for (j = 0; j < residuals->size; j++)
+       {
+         lag = (i > j) ? (i - j) : (j - i);
+         tmp = gsl_vector_get (s->arma->acf, lag);
+         gsl_matrix_set (gamma, i, j, tmp);
+         gsl_matrix_set (gamma, j, i, tmp);
+       }
+    }
+  for (i = 0; i < gamma->size1; i++)
+    {
+      fprintf (stderr, "%lf ", gsl_matrix_get (gamma, 0,i ));
+    }
+  fprintf (stderr, "\n");
+  for (i = 0; i < coef->size; i++)
+    {
+      fprintf (stderr, "%lf ", gsl_vector_get(coef, i));
+    }
+  fprintf (stderr, "\nmse=%lf\n", s->arma->mse);
+  rc = gsl_linalg_cholesky_decomp (gamma);
+
+  assert (rc != GSL_EDOM);
+  for (i = 0; i < scale->size; i++)
+    {
+      tmp = gsl_matrix_get (gamma, i, i);
+      gsl_vector_set (scale, i, tmp * tmp);
+      /* 
+        Scale the lower triangle and set the upper.
+      */
+      for (j = 0; j <= i; j++)
+       {
+         y = gsl_matrix_get (gamma, i, j) / fabs (tmp);
+         gsl_matrix_set (gamma, i, j, y);
+         gsl_matrix_set (gamma, j, i, y);
+       }
     }
-  r /= s->residuals->size;
-  r += log (s->arma->mse);
-  return r;
+  gsl_linalg_cholesky_solve (gamma, s->data, residuals);
+  for (i = 0; i < residuals->size; i++)
+    {
+      tmp = gsl_vector_get (residuals, i);
+      mse = gsl_vector_get (scale, i);
+      assert (mse > 0.0);
+      S += tmp * tmp / mse;
+      log_normalizer += log (mse);
+    }
+  s->arma->mse = S / residuals->size;
+  result = log (s->arma->mse) + log_normalizer / residuals->size;
+
+  gsl_matrix_free (gamma);
+  gsl_vector_free (scale);
+  gsl_vector_free (residuals);
+
+  return result;
 }
 
 void
-pspp_arma_d_loglikelihood (const gsl_vector *coef, void *parms, gsl_vector *g)
+pspp_arma_d_loglikelihood (const gsl_vector *coef, void *parms, gsl_vector 
*grad)
 {
-  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;
+  size_t k;
+  size_t lag;
+  size_t order;
+  int rc;
   struct pspp_arma_state *s;
+  double tmp;
+  double mse;
+  double y;
+  gsl_matrix *gamma; /* covariance */
+  gsl_vector *scale; /* mean-squared errors */
+  gsl_vector *residuals;
 
   s = (struct pspp_arma_state *) parms;
+  acf_vector (s, coef);
+  order = (s->arma->ar_order > s->arma->ma_order) ? s->arma->ar_order : 
s->arma->ma_order;
   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++)
+  gsl_vector_set_zero (grad);
+  residuals = gsl_vector_calloc (s->data->size);
+  gamma = gsl_matrix_calloc (residuals->size, residuals->size);
+  scale = gsl_vector_calloc (residuals->size);
+  for (i = 0; i < residuals->size; i++)
+    {
+      for (j = 0; j < residuals->size; j++)
+       {
+         lag = (i > j) ? (i - j) : (j - i);
+         tmp = gsl_vector_get (s->arma->acf, lag);
+         gsl_matrix_set (gamma, i, j, tmp);
+         gsl_matrix_set (gamma, j, i, tmp);
+       }
+    }
+  rc = gsl_linalg_cholesky_decomp (gamma);
+  assert (rc != GSL_EDOM);
+  for (i = 0; i < scale->size; i++)
+    {
+      tmp = gsl_matrix_get (gamma, i, i);
+      gsl_vector_set (scale, i, tmp * tmp);
+      /* 
+        Scale the lower triangle and set the upper.
+      */
+      for (j = 0; j <= i; j++)
        {
-         tmp += gsl_vector_get (&resid.vector, j) * gsl_vector_get 
(&data.vector, j) /
-           gsl_vector_get (&scale.vector, j);
+         y = gsl_matrix_get (gamma, i, j) / fabs (tmp);
+         gsl_matrix_set (gamma, i, j, y);
+         gsl_matrix_set (gamma, j, i, y);
        }
-      gsl_vector_set (g, i + s->arma->ar_order, -2.0 * tmp);
     }
+  gsl_linalg_cholesky_solve (gamma, s->data, residuals);
+  for (i = order; i < residuals->size; i++)
+    {
+      tmp = gsl_vector_get (residuals, i);
+      mse = gsl_vector_get (scale, i);
+      assert (mse > 0.0);
+      for (j = 0; j < s->arma->ar_order; j++)
+       {
+         gsl_vector_set (grad, j, 
+                         -2.0 * tmp * gsl_vector_get (s->data, i - j - 1) / 
mse);
+       }
+      for (j = s->arma->ar_order; j < (s->arma->ar_order + s->arma->ma_order); 
j++)
+       {
+         k = j - s->arma->ar_order;
+         gsl_vector_set (grad, j, -2.0 * tmp * 
+                         gsl_vector_get (residuals, i - k - 1) / mse);
+       }
+    }
+  gsl_matrix_free (gamma);
+  gsl_vector_free (scale);
+  gsl_vector_free (residuals);
 }
 /*
   Likelihood and its gradient together.
@@ -301,8 +548,100 @@
 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);
+  size_t i;
+  size_t j;
+  size_t k;
+  size_t lag;
+  size_t order;
+  int rc;
+  struct pspp_arma_state *s;
+  double tmp;
+  double S = 0.0;
+  double log_normalizer = 0.0;
+  double mse;
+  double y;
+  gsl_matrix *gamma; /* covariance */
+  gsl_vector *scale; /* mean-squared errors */
+  gsl_vector *residuals;
+
+  s = (struct pspp_arma_state *) parms;
+  acf_vector (s, coef);
+  order = (s->arma->ar_order > s->arma->ma_order) ? s->arma->ar_order : 
s->arma->ma_order;
+  coef_to_arma (coef, s->arma);
+  gsl_vector_set_zero (grad);
+  residuals = gsl_vector_calloc (s->data->size);
+  gamma = gsl_matrix_calloc (residuals->size, residuals->size);
+  scale = gsl_vector_calloc (residuals->size);
+  for (i = 0; i < residuals->size; i++)
+    {
+      for (j = 0; j < residuals->size; j++)
+       {
+         lag = (i > j) ? (i - j) : (j - i);
+         tmp = gsl_vector_get (s->arma->acf, lag);
+         gsl_matrix_set (gamma, i, j, tmp);
+         gsl_matrix_set (gamma, j, i, tmp);
+       }
+    }
+  for (i = 0; i < gamma->size1; i++)
+    {
+      fprintf (stderr, "%lf ", gsl_matrix_get (gamma, 0,i ));
+    }
+  fprintf (stderr, "\n");
+  for (i = 0; i < coef->size; i++)
+    {
+      fprintf (stderr, "%lf ", gsl_vector_get(coef, i));
+    }
+  fprintf (stderr, "\nmse=%lf\n", s->arma->mse);
+  rc = gsl_linalg_cholesky_decomp (gamma);
+
+  assert (rc != GSL_EDOM);
+  for (i = 0; i < scale->size; i++)
+    {
+      tmp = gsl_matrix_get (gamma, i, i);
+      gsl_vector_set (scale, i, tmp * tmp);
+      /* 
+        Scale the lower triangle and set the upper.
+      */
+      for (j = 0; j <= i; j++)
+       {
+         y = gsl_matrix_get (gamma, i, j) / fabs (tmp);
+         gsl_matrix_set (gamma, i, j, y);
+         gsl_matrix_set (gamma, j, i, y);
+       }
+    }
+  gsl_linalg_cholesky_solve (gamma, s->data, residuals);
+  for (i = 0; i < order; i++)
+    {
+      tmp = gsl_vector_get (residuals, i);
+      mse = gsl_vector_get (scale, i);
+      assert (mse > 0.0);
+      S += tmp * tmp / mse;
+      log_normalizer += log (mse);
+    }
+  for (i = order; i < residuals->size; i++)
+    {
+      tmp = gsl_vector_get (residuals, i);
+      mse = gsl_vector_get (scale, i);
+      assert (mse > 0.0);
+      for (j = 0; j < s->arma->ar_order; j++)
+       {
+         gsl_vector_set (grad, j, 
+                         -2.0 * tmp * gsl_vector_get (s->data, i - j - 1) / 
mse);
+       }
+      for (j = s->arma->ar_order; j < (s->arma->ar_order + s->arma->ma_order); 
j++)
+       {
+         k = j - s->arma->ar_order;
+         gsl_vector_set (grad, j, -2.0 * tmp * 
+                         gsl_vector_get (residuals, i - k - 1) / mse);
+       }
+      S += tmp * tmp / mse;
+      log_normalizer += log (mse);
+    }
+  s->arma->mse = S / residuals->size;
+  *like = log (s->arma->mse) + log_normalizer / residuals->size;
+  gsl_matrix_free (gamma);
+  gsl_vector_free (scale);
+  gsl_vector_free (residuals);
 }
 /*
   COEF stores the autoregressive coefficients, then the
@@ -335,7 +674,7 @@
       gsl_vector_set (coef, i, gsl_vector_get (a->ma_coeff, i - a->ar_order));
     }
 }
-
+#if 0
 /*
   This function assumes DATA is the reparameterized process.  This
   function will transform to give the residuals for the original
@@ -379,13 +718,14 @@
   get_mse (as);
   gsl_vector_free (tmp);
 }
+#endif
 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);
 
@@ -402,13 +742,20 @@
   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;  
 }
-
+static void
+demean (gsl_vector *x)
+{
+  double m = 0.0;
+  size_t i;
+  for (i = 0; i < x->size; i++)
+    {
+      m -= gsl_vector_get (x, i);
+    }
+  m /= x->size;
+  gsl_vector_add_constant (x, m);
+}
 /*
   Find the maximum likelihood estimates of an ARMA process. The
   function uses the EM algorithm with a gradient-descent.  The
@@ -431,17 +778,20 @@
   double epsabs = 1e-2;
   double m_stage_epsabs = 0.1;
   int rc;
+  int iter = 0;
 
+  demean (data);
   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);
+  acf_vector (as, coef);
+  s = gsl_multimin_fdfminimizer_alloc (t, order);
  
-  step_size = 0.1 * gsl_blas_dnrm2 ((const gsl_vector *) as->residuals);  
+  step_size = as->arma->mse;  
   rc = gsl_multimin_fdfminimizer_set (s, m, coef, step_size, tol);
   assert (rc != GSL_EBADLEN);
   
@@ -450,32 +800,23 @@
   /*
     EM algorithm.
    */
-  while (gsl_multimin_test_gradient (gradient, epsabs) ==
-        GSL_CONTINUE)
-    {
-      while (gsl_multimin_test_gradient (gradient, m_stage_epsabs) 
-            == GSL_CONTINUE)
+  while (iter < 20 && (gsl_multimin_test_gradient (gradient, epsabs) ==
+                      GSL_CONTINUE))
        {
-         /* Minimization */
          gsl_multimin_fdfminimizer_iterate (s);
          gradient = gsl_multimin_fdfminimizer_gradient (s);
+      iter++;
        }
-#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);
+  coef = gsl_multimin_fdfminimizer_x (s);
   gsl_vector_free (as->inf_ma_coeff);
   gsl_vector_free (as->scale);
+  printf ("%lf\t%lf\t%lf\n", gsl_vector_get (coef, 0), gsl_vector_get (coef, 
1), gsl_vector_get (coef,2));
   gsl_vector_free (coef);
   result = as->arma;
   free (as);
+/*   gsl_multimin_fdfminimizer_free (s); */
   free (m);
+
   return result;
 }
 bool