One-step ahead predicitons in a Bayesian state-space model 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.
library(rjags)
load.module("glm") 
set.seed(1)
y <- as.numeric(arima.sim(n = 100, list(ar = c(0.8)), sd=1))

# JAGS model
mod<-
  "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}$.

 A: 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:

*

*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)})$.


*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.


*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.


*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)})$.


*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.
