How to practically forecast a time series in a Bayesian framework

Question

can someone explains me how bayesian time series forecasting practically work?

Let me use the model on which I am working as example. It is a dynamic Poisson model or INGARCH(1,1) model \begin{equation} Y_t \sim Pois(\lambda_t) \\ \text{with} \qquad \lambda_t = \mu + \alpha Y_{t-1} + \beta \lambda_{t-1} \end{equation} In short, this model assumes that the data follows a Poisson whose mean/variance $\lambda_t$ changes over time in a way similar to an ARMA(1,1).

I have estimated the parameters of this model (and their distributions) using a MCMC algorithm since the posterior it is not tracatble.

What I want to know is how I can predict future distributions for multiple step-ahead. I would also like to know which is the so-called predictive distribution in a context like this.

I have tried it in the following way:

For time $t+1$:

1. I randomly choose values for $\mu, \alpha, \beta$ from the sample obtained from the MCMC algorithm
2. I use this values to compute $\lambda_{t+1}$ *
3. I draw from a $Pois(\lambda_{t+1})$
4. I repeat step 1-2-3 many times

For time $t+2$:

1. I randomly choose values for $\mu, \alpha, \beta$ from the sample obtained from the MCMC algorithm
2. I randomly choose $\lambda_{t+1}$ from the sample obtained for time $t+1$
3. I use this values to compute $\lambda_{t+2}$
4. I draw from a $Pois(\lambda_{t+2})$
5. I repeat step 1-2-3-4 many times

However, I am not confident on the correctness of this procedure.

Thanks a lot for your help.

• I recover $\lambda_{t}$ in this way:

for (i in 1:LastObservation) {
  
  if (i==1) {
    lambda[i] = mu+ alpha * InitialValueY + beta * InitialValueLambda
  } else {
    lambda[i] = mu + alpha * Y[i-1] + beta * lambda[i-1]
  }
}

#in this case lambda_t is lambda[LastObservation]

Answer 1

I will express your data as $Y_{0:T} = (Y_0,\ldots,Y_T)$ where $Y_0$ is the initial value. I'm also going to assume there is a prior distribution for $\lambda_0$. (I'm guessing you made up a values for $Y_0$ and $\lambda_0$. That doesn't affect my suggested solution.)

I'm not sure exactly how a sampler works for your model, but let's assume you have a valid sampler. I will represent the $R$ draws from the sampler as follows: \begin{equation} \{(\lambda_{0:T}^{(r)},\mu^{(r)},\alpha^{(r)},\beta^{(r)})\}_{r=1}^R . \end{equation}

The goal is to compute the following joint predictive distribution: \begin{equation} p(Y_{T+1:T+h}|Y_{0,T}) , \end{equation} where $Y_{T+1:T+h} = (Y_{T+1},\ldots,Y_{T+h})$ and $h$ is the number of steps ahead we wish to forecast. We do this by making draws of $Y_{T+1:T+h}^{(r)}$ for $r =1, \ldots, R$.

For each $r$ we can proceed as follows: First compute \begin{equation} \lambda_{T+1}^{(r)} = \mu^{(r)} + \alpha^{(r)}\,Y_T + \beta^{(r)}\,\lambda_T^{(r)} \end{equation} and draw \begin{equation} Y_{T+1}^{(r)} \sim \textsf{Poisson}(\lambda_{T+1}^{(r)}) . \end{equation} Next compute \begin{equation} \lambda_{T+2}^{(r)} = \mu^{(r)} + \alpha^{(r)}\,Y_{T+1}^{(r)} + \beta^{(r)}\,\lambda_{T+1}^{(r)} \end{equation} and draw \begin{equation} Y_{T+2}^{(r)} \sim \textsf{Poisson}(\lambda_{T+2}^{(r)}) . \end{equation} Repeat these steps until you have drawn $Y_{T+h}^{(r)}$.

This procedure produces $R$ draws from the desired distribution.