# How to practically forecast a time series in a Bayesian framework

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 $$$$Y_t \sim Pois(\lambda_t) \\ \text{with} \qquad \lambda_t = \mu + \alpha Y_{t-1} + \beta \lambda_{t-1}$$$$ 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]

• It would appear you also need a distribution for $\lambda_t$ assuming it is not observed.
– mef
Commented May 18, 2021 at 16:01
• As far as I understand, once obtained $\mu, \alpha, \beta$, $\lambda_t$ is completely determined since $Y_{t-1}$ and $\lambda_{t-1}$ are known given the information up to time t (I have added in the text how I revover $\lambda_{t}$) Commented May 19, 2021 at 16:04

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: $$$$\{(\lambda_{0:T}^{(r)},\mu^{(r)},\alpha^{(r)},\beta^{(r)})\}_{r=1}^R .$$$$
The goal is to compute the following joint predictive distribution: $$$$p(Y_{T+1:T+h}|Y_{0,T}) ,$$$$ 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 $$$$\lambda_{T+1}^{(r)} = \mu^{(r)} + \alpha^{(r)}\,Y_T + \beta^{(r)}\,\lambda_T^{(r)}$$$$ and draw $$$$Y_{T+1}^{(r)} \sim \textsf{Poisson}(\lambda_{T+1}^{(r)}) .$$$$ Next compute $$$$\lambda_{T+2}^{(r)} = \mu^{(r)} + \alpha^{(r)}\,Y_{T+1}^{(r)} + \beta^{(r)}\,\lambda_{T+1}^{(r)}$$$$ and draw $$$$Y_{T+2}^{(r)} \sim \textsf{Poisson}(\lambda_{T+2}^{(r)}) .$$$$ Repeat these steps until you have drawn $$Y_{T+h}^{(r)}$$.
This procedure produces $$R$$ draws from the desired distribution.
• Thanks a lot, it works, it confirms the forecast computed automatically by JAGS and by the frequentist version of the model. Sorry if I did not understand your comment, it was probably already clear but I have realized only now that adding a prior for $\lambda_0$ only, makes all the $\lambda_{t=1:T}$ stochastic. Commented May 20, 2021 at 16:45