Faster representation of AR(1) model for stochastic volatility?

Hello!

Context

The Stan manual gives excellent code for a stochastic volatility model here: 2.5 Stochastic volatility models | Stan User’s Guide

transformed parameters {
  vector[T] h = h_std * sigma;  // now h ~ normal(0, sigma)
  h[1] /= sqrt(1 - phi * phi);  // rescale h[1]
  h += mu;
  for (t in 2:T) {
    h[t] += phi * (h[t - 1] - mu);
  }
}

I wonder if there isn’t a way to avoid the loop here, thereby making the process amenable to vectorisation.

Conditional moments for AR(1) model

Suppose I have an AR(1) model given by:

X_t - \mu = \phi(X_{t-1} - \mu) + \epsilon_t.

Knowing the value at some value, the first mean and variance are given by the following

\mathrm{E}(X_t | X_{t-L}) = \mu(1 - \phi^L) + \phi^L X_{t-L} , and
\mathrm{Var}(X_t | X_{t-L}) = \frac{\sigma^2}{1-\phi^2}\left(1 - \phi^{2L-1}\right).

See eg https://mfe.baruch.cuny.edu/wp-content/uploads/2014/12/TS_Lecture1_2019.pdf equation 16.

This result is obtained by expanding the AR(1) model recursively, then using the standard result for the finite sum of a geometric series.

Faster in Stan?

Given an initial value, the results above give a direct formula for the conditional mean and standard deviation at any point in the AR(1) model. Can you validly use this to “directly” specify the mean and standard deviation at each time step, thereby avoiding the loop?

I guess this assumes that you don’t have knowledge of the intermediate points, which is not true in general, because actually you know X_{t-1} etc.

This also assumes the process is AR(1), so I’m not sure what impact model misspecification have.

Fallback

At the very least, there’s a common sub-expression that can be pulled out:

transformed parameters {
  vector[T] h = h_std * sigma;  // now h ~ normal(0, sigma)
  real phi_times_mu = phi*mu;
  h[1] /= sqrt(1 - phi * phi);  // rescale h[1]
  h += mu;
  for (t in 2:T) {
    //h[t] += phi * (h[t - 1] - mu);
   h[t] += phi * h[t - 1] - phi_times_mu;
  }
}

which in turn can be further simplified to:

transformed parameters {
  vector[T] h = h_std * sigma;  // now h ~ normal(0, sigma)
  real phi_times_mu = phi*mu;
  h[1] /= sqrt(1 - phi * phi);  // rescale h[1]
  h += mu - phi_times_mu;
  for (t in 2:T) {
    h[t] += phi * h[t - 1];
  }
}

It still feels like there should be a way to make this more efficient and vectorise it by relying on geometric series results.

1 Like

The entire h vector is also multivariate normal with a tridiagonal covariance matrix. So you can also write it in that form if you choose to sample in terms of the h instead of the v. Tridiagonal covariance matrices are super fast to work with since you can get their determinants and inverses in O(n) time.

Although writing it the way the user’s guide does is convenient for expanding the model e.g. using student-t distributed innovations instead of gaussian.

EDIT: I said tridiagonal covariance matrix, but it’s actually tridiagonal precision matrix.

1 Like

One problem with this stochastic volatility approach in Stan is that the number of things you need to estimate in the model will increase with the number of observations. Both GARCH and SV models have this problem in Stan, but ARCH-like models tend not to. What the GARCH/SV models do better with is capturing the partial autocorrelation of squared errors. So what you can do is have an ARCH(N) model where N is small and then include some additional function that attempts to capture the partial autocorrelation of squared errors better in a parsimonious kind of way.

Thanks - this is very insightful.