Functions to ensure stationarity in autoregression and invertibility in moving average parameters

[jrnold]

Summary:

Provide function to constrain a vector to the space of coefficients of stationary autoregressive functions and invertible moving-average functions.

Description:

When estimating ARMA time series models they are not stationary if the roots of the characteristic polynomial of the AR part lie outside the unit circle, or the roots of the characteristic polynomial of the MA part lie outside the unit circle. Since this space gets complicated once the order of those polynomials gets beyond 2, it would be useful to have functions that convert from the real line to this space, and the inverse. This issue is mentioned on p 95 of the v2.10.0 manual, and the constraint is manually imposed in the GARCH(1,1) model on p. 90. These functions would make it possible to ensure stationary of the coefficients of ARMA(p, q) and related time series models.

Below I provide Stan functions that implement the necessary transformations, described in Monahan (1984) and Jones (1987). I could reimplement them in C++ if this is of interest.

References

• Monahan, John F. 1984. "A Note on Enforcing Stationarity in Autoregressive-moving Average Models." Biometrika 71 (2) (August 1): 403-404
• Jones, R. H. (1980) Maximum likelihood fitting of ARMA models to time series with missing observations. Technometrics 22 389–395.
• Jones, 1987. Randomly choosing parameters from the stationarity and invertibility region of autoregressive-moving average models. Applied Statistics, 36.
• functions constrain_stationary_univariate and unconstrain_stationary_univariate in the Python package statsmodels.
• R function arima, with the transformations implemented in C with functions partrans and invpartrans

Expected Output:

Additional Information:

Stan functions implementing these transformations.

/*
Partial Autocorrelations to Autocorrelations

@param vector of coefficients of a partial autocorrelation function
@return vector of coefficients of an Autocorrelation function

*/
vector pacf_to_acf(vector x) {
  matrix[num_elements(x), num_elements(x)] y;
  int n;
  n <- num_elements(x);
  y <- rep_matrix(0.0, n, n);
  for (k in 1:n) {
    for (i in 1:(k - 1)) {
      y[k, i] <- y[k - 1, i] + x[k] * y[k - 1, k - i];
    }
    y[k, k] <- x[k];
    print(y);
  }
  return -y[n] ';
}

/*
Constrain vector of coefficients to the stationary and intertible region for AR or MA functions.

@param vector x An unconstrained vector in (-Inf, Inf)
@returns vector y A vector of coefficients for a stationary AR or inverible MA process.

*/
vector constrain_stationary(vector x) {
  vector[num_elements(x)] r;
  int n;
  n <- num_elements(x);
  // transform (-Inf, Inf) to (-1, 1)
  for (i in 1:n) {
    r[i] <- x[i] / (sqrt(1.0 + pow(x[i], 2)));
  }
  // Transform PACF to ACF
  return pacf_to_acf(r);
}

/*
Convert coefficients of an autocorrelation function to partial autocorrelations.

@param vector x Coeffcients of an autocorrelation function.
@returns vector y Coefficients of the corresponding partial autocorrelation function.

*/
vector acf_to_pacf(vector x) {
  matrix[num_elements(x), num_elements(x)] y;
  vector[num_elements(x)] r;
  int n;
  n <- num_elements(x);
  y <- rep_matrix(0.0, n, n);
  y[n] <- -x ';
  for (j in 0:(n - 1)) {
    int k;
    k <- n - j;
    for (i in 1:(k - 1)) {
      y[k - 1, i] <- (y[k, i] - y[k, k] * y[k, k - i]) / (1 - pow(y[k, k], 2));
    }
  }
  r <- diagonal(y);
  return r;
}

/*
Transform from stationary and invertible space to (-Inf, Inf)

@param vector x Coeffcients of an autocorrelation function.
@returns vector y Coefficients of the corresponding partial autocorrelation function.

*/
vector unconstrain_stationary(vector x) {
  matrix[num_elements(x), num_elements(x)] y;
  vector[num_elements(x)] r;
  vector[num_elements(x)] z;
  int n;
  n <- num_elements(x);
  // Transform ACF to PACF
  r <- acf_to_pacf(x);
  // Transform (-1, 1) to (-Inf, Inf)
  for (i in 1:n) {
    z[i] <- r[i] / (sqrt(1.0 - pow(r[i], 2)));
  }
  return z;
}

Provide any additional information here.

Current Version:

v2.9.0

[bgoodri]

Yeah, we have needed these for a long time.

[jrnold]

It would be awesome to have a stationary type, and the jacobians are available in (one of the many) sources of this transformation:

Monahan, John F. 1984. “A Note on Enforcing Stationarity in Autoregressive-moving Average Models.” Biometrika 71 (2) (August 1): 403-404.

Also, it is probably less in demand, but statsmodels has functions that extend it to the multivariate case.

[bob-carpenter]

If you can implement the transfrom from unconstrained to constrained and the log Jacobian, I can add the type to the language. We can speed up a lot of those functions above with tighter bounds and vectorization, but we can look at it in C++ rather than trying to make the Stan functions faster.

[bgoodri]

What should we call the type? Maybe stationary_vector or AR_vector or
ARMA_vector (this transformation also works for MA weights).

On Mon, Jul 4, 2016 at 4:34 PM, Bob Carpenter notifications@github.com
wrote:

If you can implement the transfrom from unconstrained to constrained and
the log Jacobian, I can add the type to the language. We can speed up a lot
of those functions above with tighter bounds and vectorization, but we can
look at it in C++ rather than trying to make the Stan functions faster.

—
You are receiving this because you commented.
Reply to this email directly, view it on GitHub
#309 (comment), or mute
the thread
https://github.com/notifications/unsubscribe/ADOrqqM0rzwwrztUlNnhpXMWgE7ZsnU4ks5qSW5ZgaJpZM4I-dK1
.

[bob-carpenter]

How about stationary_ar with the understanding that it's
a vector? It'd be like simplex or positive_ordered that way.

• Bob

On Jul 4, 2016, at 1:38 PM, bgoodri notifications@github.com wrote:

What should we call the type? Maybe stationary_vector or AR_vector or
ARMA_vector (this transformation also works for MA weights).

On Mon, Jul 4, 2016 at 4:34 PM, Bob Carpenter notifications@github.com
wrote:

If you can implement the transfrom from unconstrained to constrained and
the log Jacobian, I can add the type to the language. We can speed up a lot
of those functions above with tighter bounds and vectorization, but we can
look at it in C++ rather than trying to make the Stan functions faster.

—
You are receiving this because you commented.
Reply to this email directly, view it on GitHub
#309 (comment), or mute
the thread
https://github.com/notifications/unsubscribe/ADOrqqM0rzwwrztUlNnhpXMWgE7ZsnU4ks5qSW5ZgaJpZM4I-dK1
.

—
You are receiving this because you were assigned.
Reply to this email directly, view it on GitHub, or mute the thread.

[bgoodri]

Maybe but stationary_ar is not my fav because y is stationary if these
coefficients generate it, rather than the coefficients themselves being
stationary. Also, we need to implement the informative prior on this vector
because the usual beta informative prior is put on the partial
correlations, which the user doesn't have access to.

On Mon, Jul 4, 2016 at 4:41 PM, Bob Carpenter notifications@github.com
wrote:

How about stationary_ar with the understanding that it's
a vector? It'd be like simplex or positive_ordered that way.

• Bob

On Jul 4, 2016, at 1:38 PM, bgoodri notifications@github.com wrote:

What should we call the type? Maybe stationary_vector or AR_vector or
ARMA_vector (this transformation also works for MA weights).

On Mon, Jul 4, 2016 at 4:34 PM, Bob Carpenter notifications@github.com
wrote:

If you can implement the transfrom from unconstrained to constrained
and
the log Jacobian, I can add the type to the language. We can speed up
a lot
of those functions above with tighter bounds and vectorization, but we
can
look at it in C++ rather than trying to make the Stan functions faster.

—
You are receiving this because you commented.
Reply to this email directly, view it on GitHub
#309 (comment),
or mute
the thread
<
https://github.com/notifications/unsubscribe/ADOrqqM0rzwwrztUlNnhpXMWgE7ZsnU4ks5qSW5ZgaJpZM4I-dK1

.

—
You are receiving this because you were assigned.

Reply to this email directly, view it on GitHub, or mute the thread.

—
You are receiving this because you commented.
Reply to this email directly, view it on GitHub
#309 (comment), or mute
the thread
https://github.com/notifications/unsubscribe/ADOrqooDuyijLM8PxVJ0sdhq3S4n3Oeeks5qSW_pgaJpZM4I-dK1
.

[jrnold]

I can't think of a good name for the type either, I had a hard enough time trying to write the title of this issue. Maybe stationary_arcoef of stationary_acf ?

For the informative prior, we could create a atanh-beta distribution, where X ~ atanh_beta(a, b) = atanh((X + 1) / 2) ~ beta(a, b).

[bob-carpenter]

I'd be OK with stationary_ar or stationary_ar_coef.

[betanalpha]

stationary_basis?

Having to write out “stationary” is a bit burdensome, so what aobut
sar_basis ro sar_coef?

On Jul 5, 2016, at 2:49 AM, Bob Carpenter notifications@github.com wrote:

I'd be OK with stationary_ar or stationary_ar_coef.

—
You are receiving this because you are subscribed to this thread.
Reply to this email directly, view it on GitHub, or mute the thread.

[jrnold]

But that's what autocomplete is for

[jrnold]

It'll take me a little bit of time to get re-familiarized with the C++ library since it's been so long since I looked at it, so I'll have a few questions.

[bgoodri]

See also https://arxiv.org/abs/1406.4584 for the vector case

[elbamos]

I'm curious whatever happened to this? It would definitely be nice to have.

[elbamos]

If anyone else comes across this page and attempts to use the code - I'm pretty sure the pacf_to_acf function should be:

  row_vector pacf_to_acf(vector x) {
    int n = num_elements(x);
    matrix[n, n] y = diag_matrix(x);

    for (k in 2:n) {
      for (i in 1:(k - 1)) {
        y[k, i] = y[k - 1, i] - x[k] * y[k - 1, k - i];
      }
    }
    return y[n];
  }

and you don't need the inverse or the log abs det jacobian if you set your priors on the auto-correlations and constrain them to [-1, 1].

[bob-carpenter]

@elbamos --- thanks for th clarification.