How can I make one-step ahead predictions in a Bayesian state-space model? Specifically, for the model

\begin{align} y_t \sim& N(x_t, py) \\ x_t \sim& N(rx_{t-1}, px) \end{align}

where $y_t$ are the observations, and $x_t$ the states, I want to predict $\Pr(x_{t+1} \mid y_{1:t})$ or $\Pr(y_{t+1} \mid x_{1:t}, y_{1:t})$.

I have written some JAGS code below for a Gaussian state-space model but I do not know how to extend it to make predictions. I have tried to predict $x_{t+1}$ in the code below with xpred but it just seems to nudge the time-series along so I do not think it is implemented correctly. I do not know how to predict $\Pr(y_{t+1})$. I'd really welcome any help please (not necessarily complete code but some implementation hints to get me moving). Thank you.

y <- as.numeric(arima.sim(n = 100, list(ar = c(0.8)), sd=1))

# JAGS model
    x[1] ~ dnorm(0, 1)
    y[1] ~ dnorm(x[1], py)
    for (t in 2:N)  {
      x[t] ~ dnorm(r*x[t-1], px)
      y[t] ~ dnorm(x[t], py)
    # This is what I added for the predictions of x[t+1]
    for(i in 2:(N+1)){
      xpred[i] ~ dnorm(r*x[i-1], px)
    py ~ dgamma(1,1) 
    r ~ dnorm(0,0.1)
    px ~ dgamma(1,1)

# Sample
jm <- jags.model(textConnection(mod), data = list(y=y, N=length(y)), n.chain=2)
js <- coda.samples(jm, variable.names=c('x','xpred','px','py','r'), n.iter=1000)

A plot of the time-series observations $y_t$, the states $x_t$, and predictions $x_{t+1}$.

enter image description here


1 Answer 1


When you forecast, the object you ultimately need is the posterior predictive distribution:

$$p(Y_{t+1} | Y_{1:t})$$

In a state space model with state $X_t$ and parameters $\theta$, the specific structure of the model allows you to decompose this as follows:

$$p(Y_{t+1} | Y_{1:t}) = \int\int\int p(Y_{t+1}|X_{t+1}, \theta)p(X_{t+1}|X_t,\theta)p(X_t|Y_{1:t},\theta)p(\theta|Y_{1:t}) d\theta dX_t dX_{t+1}$$

This suggests the following algorithm:

  1. Obtain a sample $\theta^{(i)}$ from the posterior distribution of the parameters $p(\theta|Y_{1:t})$. This is the part you've already done (via JAGS), where $\theta^{(i)} = (r^{(i)}, px^{(i)}, py^{(i)})$.

  2. Obtain a sample $X_t^{(i)}$ for the final state from the filtering distribution $p(X_t|Y_{1:t},\theta)$, conditional on the same $\theta^{(i)}$ you've just sampled in 1. You can do this with the Kalman filter but you have also already obtained this from JAGS.

  3. Obtain a sample from $X_{t+1}^{(i)}$ from the state transition distribution $p(X_{t+1}|X_t,\theta)$, conditional on the same $\theta^{(i)}$ and the same $X_t^{(i)}$ from 1-2.

    In your specific example, that would be $X_{t+1}^{(i)} \sim \mathcal{N}(r^{(i)}X_t^{(i)}, px^{(i)})$. I don't use JAGS but what you have there looks correct. You could also just do it by hand afterwards. The plot doesn't look obviously wrong to me for an AR(1) model with mean zero.

  4. Obtain a sample from $Y_{t+1}^{(i)}$ from the measurement distribution $p(Y_{t+1}|X_{t+1}, \theta)$, conditional on the same $\theta^{(i)}$ and the same $X_{t+1}^{(i)}$ from 1-3. In your specific example, that's $Y_{t+1}^{(i)} \sim \mathcal{N}(X_{t+1}^{(i)}, py^{(i)})$.

  5. Repeat for $i = 1, ..., N$, the number of samples you have from JAGS. The $Y_{t+1}^{(i)}$ are then a sample from $p(Y_{t+1} | Y_{1:t})$.

The hardest part is getting the posterior for the parameters and the state (which you have from JAGS), everything after that is relatively simple: conditional on that, just sample a transition to the new state, then sample a measurement of the new state.

Edit to add: The above works correctly for an out of sample prediction, at the end of the data. It would also work at any point in time, aside from the fact that you do not have the parameter posterior $p(\theta| Y_{1:t})$ for every $t$, just for the last one. You would have to re-run JAGS on increasing subsamples or use a different algorithm to get the others. I just noticed that you have plotted x_pred at every time point so wanted to make that clear.

  • $\begingroup$ Thank you for your answer Chris, I appreciate you going through it step by step. Also thanks for the comment on the plot; I don't have the experience to judge the fit of these models and so from my untrained eye it seems as if some bias (probably not the correct word) was getting introduced as they seemed one step out. $\endgroup$ Jul 12, 2020 at 14:46
  • 1
    $\begingroup$ The mean of the transition distribution at $t+1$ is $r X_t$. Since $r = 0.8$, this is a number which is close to $X_t$, but shrunk slightly towards zero. That's why the time series of predictions looks like that. $\endgroup$
    – Chris Haug
    Jul 12, 2020 at 14:53
  • $\begingroup$ Chris thank you for your edit. It seems as if I have completely misunderstood what I was doing in JAGS. Do I understand you correctly that I cannot get one-step ahead predictions on the fly like this but would require to fit the model $N_{obs}$ times (or whatever chunk size of data I use) $\endgroup$ Jul 12, 2020 at 15:06
  • 1
    $\begingroup$ If you are trying to evaluate actual out of sample forecasting performance, then no, a single MCMC run like that won't work, because the parameter and state posterior then depends on "future" data. $\endgroup$
    – Chris Haug
    Jul 12, 2020 at 15:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.