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

jrnold commented 8 years ago

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 commented 8 years ago

Yeah, we have needed these for a long time.

jrnold commented 8 years ago

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 commented 8 years ago

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 commented 8 years ago

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 https://github.com/stan-dev/math/issues/309#issuecomment-230351232, or mute the thread https://github.com/notifications/unsubscribe/ADOrqqM0rzwwrztUlNnhpXMWgE7ZsnU4ks5qSW5ZgaJpZM4I-dK1 .

bob-carpenter commented 8 years ago

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 https://github.com/stan-dev/math/issues/309#issuecomment-230351232, 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 commented 8 years ago

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 https://github.com/stan-dev/math/issues/309#issuecomment-230351232, 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 https://github.com/stan-dev/math/issues/309#issuecomment-230351938, or mute the thread https://github.com/notifications/unsubscribe/ADOrqooDuyijLM8PxVJ0sdhq3S4n3Oeeks5qSW_pgaJpZM4I-dK1 .

jrnold commented 8 years ago

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 commented 8 years ago

I'd be OK with stationary_ar or stationary_ar_coef.

betanalpha commented 8 years ago

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 commented 8 years ago

But that's what autocomplete is for

jrnold commented 8 years ago

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 commented 8 years ago

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

elbamos commented 5 years ago

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

elbamos commented 5 years ago

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 commented 5 years ago

@elbamos --- thanks for th clarification.

stan-dev / math