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.