stochastic forecasts

[ac_1]

Hi Tom,

I want to generate stochastic forecasts/simulations for (1) an AR model and (2) a VAR model, as follows:
(a) recursive static stochastic 1-step ahead forecasts via bootstrapping residuals and adding them to the static forecast.
(b) dynamic stochastic multiperiods-periods ahead forecasts via bootstrapping residuals and adding the residual initially to the static 1-step ahead forecast, and thereafter to the stochastic forecasts - as they are iterated.

The average of the RHS across EACH of the simulations being treated as the stochastic forecast.

The idea catering for the fact that a normal distribution for the forecast errors maybe an unreasonable assumption, further assumes future residuals are similar to past residuals, and related to https://otexts.com/fpp3/prediction-intervals.html.

static stochastic 1-step ahead
stoch_fcast(t+1) = static_fcast_one-step(t+1) + unknown_future_residual(t+1)

The unknown_future_residual(t+1) is bootstrap (10000 draws) sampled from the model "fitted" residuals and averaged.

and repeating,

stoch_fcast(t+2) = static_fcast_one-step(t+2) + unknown_future_residual(t+2)

where unknown_future_residual(t+2) is another bootstrap (10000 draws) sampled from the collection of residuals and averaged.

etc etc.

dynamic stochastic multiperiods-periods ahead
stoch_fcast(t+1) = static_fcast_one-step(t+1) + unknown_future_residual(t+1)
stoch_fcast(t+2) = stoch_fcast(t+1) + unknown_future_residual(t+2)
stoch_fcast(t+3) = stoch_fcast(t+2) + unknown_future_residual(t+3)

etc etc.


Here's an incorrect attempt at recursive static stochastic 1-step ahead forecasts for an AR model. I am unclear as to the RATS commands/syntax for the bootstrap?

Code: Select all

seed 13939
do regend = ibegin, iend
   boxjenk(constant,ar=1,define=eqmodel) y istart regend resids
   uforecast(equation=eqmodel,from=regend+1,steps=1) yhat_S_model
   do i = 1, 10000
      boot resids
      set ac_sample = resids(t)
      set stoch_fcasts = yhat_S_model + ac_sample
   end do i
   sstats(mean) stoch_fcasts>>stoch_fcasts_mean
end do regend
prin / stoch_fcasts_mean

many thanks,
Amarjit

[TomDoan]

You're correct that that doesn't work---all that is doing is adding noise to the point forecasts. Bootstrapping requires adding shocks at each horizon, which then feed through into later steps.

The general procedure for including bootstrapping in forecasts is FORECAST with PATHS as described in 16.9 of the User's Guide. However (as described there), for simple bootstrapping of a univariate model, you can just use the BOOTSTRAP option on UFORECAST, which draws shocks randomly (with replacement) from the residuals associated with the equation when it was estimated. So

do i=1,10000
uforecast(equation=eqmodel,from=regend+1,steps=1) stoch_fcasts
end do i

would do 10000 bootstrapped one step simulations. (Note that those really aren't forecasts---they're simulations). If you want multi-step bootstrapped simulations, you would just change the STEPS option.

However, it's not clear what you plan to do with the results. That's overwriting the same series 10000 times. You need to look at section 16.1 of the User's Guide and figure out what you want to do with the 10000 simulated values. If all you want is the mean (the point forecast), you just need the sum across simulations. (Not across time, which is what the SSTATS will do).

set stoch_fcasts = 0.0
do i=1,10000
uforecast(equation=eqmodel,from=regend+1,steps=1) stoch_sims
set stoch_fcasts = stoch_fcasts+stoch_sims
end do i
set stoch_fcasts = stoch_fcasts/10000

would give you the point estimates of the forecasts. Note that those are likely to be very close to non-stochastic forecasts. (In taking means across a large number of simulations, the lack of Gaussianity doesn't really matter much).

[ac_1]

Here's an attempt at:
(a) recursive static stochastic 1-step ahead forecasts
(b) dynamic stochastic multistep ahead iterated forecasts
for an AR(1) model.

Code: Select all

* recursive static 1-step ahead forecasts
clear stoch_fcasts
*
do regend = ibegin, iend
   boxjenk(noprint,constant,ar=1,define=eqmodel) y istart regend resids
   uforecast(equation=eqmodel,from=regend+1,steps=1) stoch_fcasts
end do regend
prin / stoch_fcasts
*
* recursive static stochastic 1-step ahead forecasts
seed 13939
clear stoch_sims
clear stoch_fcasts_mean
*
do regend = ibegin, iend
   boxjenk(noprint,constant,ar=1,define=eqmodel) y istart regend
   uforecast(equation=eqmodel,from=regend+1,steps=1) stoch_fcasts
   do i=1,ndraws
    * draws "shocks" randomly (with replacement) from the residuals associated with the equation when it was estimated
      uforecast(bootstrap,equation=eqmodel,from=regend+1,steps=1) stoch_sims
      set stoch_fcasts regend+1 regend+1 = stoch_fcasts+stoch_sims; * stochastic simulations
   end do i
   set stoch_fcasts_mean regend+1 regend+1 = stoch_fcasts/ndraws; * point forecasts
end do regend

Code: Select all

seed 13939
clear stoch_fcasts
clear stoch_fcasts_mean
*
* dynamic multistep ahead iterated forecasts
boxjenk(noprint,constant,ar=1,define=eqmodel) y istart iend resids
uforecast(equation=eqmodel,from=iend+1,steps=nsteps) stoch_fcasts
prin / resids stoch_fcasts
*
* dynamic stochastic multistep ahead iterated forecasts
do i=1,ndraws
   uforecast(bootstrap,equation=eqmodel,from=iend+1,steps=nsteps) stoch_sims
   set stoch_fcasts iend+1 iend+nsteps = stoch_fcasts+stoch_sims; * stochastic simulations
end do i
set stoch_fcasts_mean iend+1 iend+nsteps = stoch_fcasts/ndraws; * point forecasts

Alright?


To be able to generate stochastic forecasts for VAR's I need to use FORECAST.
Have attempted for an equivalent result to the AR(1) dynamic stochastic multistep ahead iterated forecasts - but I get NA's for stoch_fcasts_mean?

Code: Select all

seed 13939
clear stoch_fcasts
clear stoch_fcasts_mean
*
* dynamic stochastic multistep ahead iterated forecasts
do i=1,ndraws
   boot entries / istart iend; * BOOT fills the ENTRIES series with random integers between istart and iend
   set path1 iend+1 iend+nsteps = resids(entries); * Choose random ENTRIES from RESIDS
   forecast(paths,from=iend+1,steps=nsteps) 1
   # eqmodel stoch_fcasts; * stochastic simulations
   # path1
end do i
set stoch_fcasts_mean iend+1 iend+nsteps = stoch_fcasts/ndraws; * point forecasts

[TomDoan]

For the last, you need

boot entries iend+1 iend+nsteps istart iend; * BOOT fills the ENTRIES series with random integers between istart and iend

ENTRIES needs to be defined over the range where it's needed---the simulation range.

[ac_1]

Am still having problems with stochastic simulations for an AR(1) with FORECAST - I get NA's for stoch_fcasts & stoch_fcasts_mean.

Code: Select all

clear(len=iend+nsteps,zeros) stoch_fcasts
clear(len=iend+nsteps,zeros) stoch_sims
clear(len=iend+nsteps,zeros) stoch_fcasts_mean
clear(len=iend+nsteps,zeros) path1
*
* dynamic multistep ahead iterated forecasts
boxjenk(noprint,constant,ar=1,define=eqmodel) y istart iend resids
uforecast(equation=eqmodel,from=iend+1,steps=nsteps) stoch_fcasts
prin / resids stoch_fcasts
*
* dynamic stochastic multistep ahead iterated forecasts
do i=1,ndraws
   boot entries iend+1 iend+nsteps istart iend
   set path1 iend+1 iend+nsteps = resids(entries)
   forecast(paths,from=iend+1,to=iend+nsteps) 1
   # eqmodel stoch_sims
   # path1
   set stoch_fcasts iend+1 iend+nsteps = stoch_fcasts+stoch_sims
end do i
set stoch_fcasts_mean iend+1 iend+nsteps = stoch_fcasts/ndraws

Similarly for a VAR.

I have estimated a VAR(1), 3 variables, and generated dynamic multistep iterated forecasts in stoch_fcasts. Then run the following to generate dynamic stochastic multistep ahead iterated forecasts, using FORECAST:

Code: Select all

* dynamic stochastic multistep ahead iterated forecasts
dec vec[series] path(3)
dec vec[series] stoch_sims(3)
dec vec[series] dfore_m(3)

clear(len=iend+nsteps,zeros) stoch_sims(1)
clear(len=iend+nsteps,zeros) stoch_sims(2)
clear(len=iend+nsteps,zeros) stoch_sims(3)
clear(len=iend+nsteps,zeros) path(1)
clear(len=iend+nsteps,zeros) path(2)
clear(len=iend+nsteps,zeros) path(3)
clear(len=iend+nsteps,zeros) dfore_m(1)
clear(len=iend+nsteps,zeros) dfore_m(2)
clear(len=iend+nsteps,zeros) dfore_m(3)

do draw=1,ndraws
   boot entries iend+1 iend+nsteps istart iend
   set path(1) iend+1 iend+nsteps = resids(1)(entries)
   set path(2) iend+1 iend+nsteps = resids(2)(entries)
   set path(3) iend+1 iend+nsteps = resids(3)(entries)
   forecast(paths,model=varmodel,results=stoch_sims,from=iend+1,to=iend+nsteps,steps=nsteps) 3
   # path(1) path(2) path(3)
   set stoch_fcasts(1) iend+1 iend+nsteps = stoch_fcasts(1)+stoch_sims(1)
   set stoch_fcasts(2) iend+1 iend+nsteps = stoch_fcasts(2)+stoch_sims(2)
   set stoch_fcasts(3) iend+1 iend+nsteps = stoch_fcasts(3)+stoch_sims(3)
end do draw

set dfore_m(1) iend+1 iend+nsteps = stoch_fcasts(1)/ndraws
set dfore_m(2) iend+1 iend+nsteps = stoch_fcasts(2)/ndraws
set dfore_m(3) iend+1 iend+nsteps = stoch_fcasts(3)/ndraws

(i) With ndraws only =10, sometimes I get all forecasts for the true out of sample period iend+1 iend+nsteps, sometimes a few NA's, sometimes all NA's.
(ii) If I increase ndraws to a reasonable amount from 10; to 100, 1000, 10000, I get all NA's.
(iii) Likewise for static stochastic 1 step ahead: there are on occasion NA's for out-of-sample forecasts, but all NA's for ndraws=100, or 1000, etc.

[TomDoan]

Works fine for me. Based upon your description, your estimation period must have missing values. That makes the residuals have (a few) NA's which will randomly get included in the bootstrapped simulations. If you do low number of NDRAWS, you're less likely to hit an NA; if you do a lot, you're almost certain to.

[ac_1]

Thanks. Yes - I changed istart and they work! I will have questions regarding prediction intervals from the simulations...

[ac_1]

You need to look at section 16.1 of the User's Guide and figure out what you want to do with the 10000 simulated values.

I have calculated stochastic SE's and generated stochastic PI's.

I have had a further read of section 16.1. From the perspective of forecasts, what else can be calculated from the simulated values?

[TomDoan]

You need to look at section 16.1 of the User's Guide and figure out what you want to do with the 10000 simulated values.

I have calculated stochastic SE's and generated stochastic PI's.

I have had a further read of section 16.1. From the perspective of forecasts, what else can be calculated from the simulated values?

You can estimate any statistic about the time path of the data through the forecast period, such as P(max value > K).

[ac_1]

You can estimate any statistic about the time path of the data through the forecast period, such as P(max value > K).

How about computing a density function for the forecasts e.g. superimposing the simulations vs the actuals, similar to the graphs in https://estima.com/ratshelp/index.html? ... bsrpf.html, but vs. the actuals?

Although there, within the Gibbs Sampler (for the parameter estimates) there is code like:

compute inter(draw)=beta(1)

but here there is:

set S_fore_m(1) begin+1 end+1 = S_stoch_fcasts(1)/ndraws; * point forecasts

So am uncertain how to plot?


What does P(max value > K) mean?

[TomDoan]

You can. I'm not sure what you think that would provide. Density graphs are most useful when the distribution is highly non-Gaussian (such as multiple peaks) and your simulations won't have that. There is no simple way to display those when you have a separate density for each data point. Some people have done 3-D graphs where each is arranged as a 2-D hill but those aren't particularly interesting, and in practice provide no information that isn't more easily seen from a fan graph.

P that the max of a series exceeds a particular value.

[ac_1]

Hi Tom,

I want to plot the PI's from a bootstrapped VAR for non-Gaussian distributions.

I have generated bootstrap forecasts as

Code: Select all

* dynamic multistep ahead iterated forecasts
dec vec[series] D_stoch_fcasts(2)
dec vec[series] D_stoch_stderrs(2)
dec vec[series] D_stoch_fcasts_sq(2)

system(model=varmodel)
variables x y
lags 1 2 3 4 
det constant
end(system)
estimate(print,resids=D_resids,noftests) istart iend
*
forecast(model=varmodel,results=D_stoch_fcasts,stderrs=D_stoch_stderrs,from=iend+1,steps=(iend+1-ibegin))

prin / D_stoch_fcasts D_stoch_stderrs

set D_stoch_fcasts_sq(1) iend+1 iend+1+nsteps = D_stoch_fcasts(1)
set D_stoch_fcasts_sq(2) iend+1 iend+1+nsteps = D_stoch_fcasts(2)

prin / D_stoch_fcasts_sq


* dynamic stochastic multistep ahead iterated forecasts
dec vec[series] D_path(2)
dec vec[series] D_stoch_sims(2)
dec vec[series] D_fore_m(2)

do draw=1,ndraws
   boot entries iend+1 iend+1+nsteps istart iend
   set D_path(1) iend+1 iend+1+nsteps = D_resids(1)(entries)
   set D_path(2) iend+1 iend+1+nsteps = D_resids(2)(entries)
   forecast(paths,model=varmodel,results=D_stoch_sims,from=iend+1,to=iend+1+nsteps)
   # D_path(1) D_path(2)

   set D_stoch_fcasts(1) iend+1 iend+1+nsteps = D_stoch_fcasts(1)+D_stoch_sims(1)
   set D_stoch_fcasts(2) iend+1 iend+1+nsteps = D_stoch_fcasts(2)+D_stoch_sims(2)
end do draw

set D_fore_m(1) iend+1 iend+1+nsteps = D_stoch_fcasts(1)/ndraws; * point forecasts
set D_fore_m(2) iend+1 iend+1+nsteps = D_stoch_fcasts(2)/ndraws; * point forecasts

But I need EACH draw to calculate something like, (this being for a single equation):

Code: Select all

set fcastMIN_D   iend+1 iend+nsteps = work=%xt(fcast_D,t),%fractiles(work,||.00||)(1); * 100%
set fcast95L_D   iend+1 iend+nsteps = work=%xt(fcast_D,t),%fractiles(work,||.025||)(1); * 95%

set fcastmean_D  iend+1 iend+nsteps = work=%xt(fcast_D,t),%avg(work); * mean
set fcastmedian_D    iend+1 iend+nsteps = work=%xt(fcast_D,t),%fractiles(work,||.50||)(1); * median OR middle quantile, i.e. 50th percentile, or 0.5 quantile

set fcast95U_D   iend+1 iend+nsteps = work=%xt(fcast_D,t),%fractiles(work,||.975||)(1); * 95%
set fcastMAX_D   iend+1 iend+nsteps = work=%xt(fcast_D,t),%fractiles(work,||1.00||)(1); * 100%

How do I extract EACH draw and calculate EACH of these fcastxxx series, per variable in a multivariate setting?


thanks,
Amarjit

[TomDoan]

Isn't it just the same thing, but with a separate SERIES[VECT} for each variable in the model?

[ac_1]

Isn't it just the same thing, but with a separate SERIES[VECT} for each variable in the model?

From the RATS manual UG-245 to UG-246, I now have something like:

Code: Select all

dec series[vect] D_fcasts_1 D_fcasts_2
gset D_fcasts_1 iend+1 iend+nsteps = %zeros(ndraws,1)
gset D_fcasts_2 iend+1 iend+nsteps = %zeros(ndraws,1)

do draw=1,ndraws
   boot entries iend+1 iend+nsteps istart iend
   set D_path(1) iend+1 iend+nsteps = D_resids(1)(entries)
   set D_path(2) iend+1 iend+nsteps = D_resids(2)(entries)
   forecast(paths,model=varmodel,results=D_stoch_sims,from=iend+1,to=iend+nsteps)
   # D_path(1) D_path(2)

   set D_stoch_fcasts(1) iend+1 iend+nsteps = D_stoch_fcasts(1)+D_stoch_sims(1)
   set D_stoch_fcasts(2) iend+1 iend+nsteps = D_stoch_fcasts(2)+D_stoch_sims(2)

   do t=iend+1,iend+nsteps
      compute D_fcasts_1(t)(draw)=D_stoch_sims(1)(t)
      compute D_fcasts_2(t)(draw)=D_stoch_sims(2)(t)
   end do t
end do draw


set D_fore_m(1) iend+1 iend+nsteps = D_stoch_fcasts(1)/ndraws; * point forecasts
set D_fore_m(2) iend+1 iend+nsteps = D_stoch_fcasts(2)/ndraws; * point forecasts

* equivalently as
set D_1_fmean iend+1 iend+nsteps = %avg(D_fcasts_1(t)); * mean
set D_2_fmean iend+1 iend+nsteps = %avg(D_fcasts_2(t)); * mean

and have computed fractiles, simiilarly.


AIM: I want to plot forecast densities, and plot PI's within a cloud of simulations in a multivariate setting.


I have done this in a single equation setting, where I have something like:

declare vector[series] fcasts_D(nsteps)
do ng = 1, nsteps
make array_fc_D begin+ng begin+ng
# fcast_D
set fcasts_D(ng) 1 ndraws = array_fc_D(1,t)
end do ng

fcast_D is VECTOR[SERIES[REAL]], size (ndraws*1)

fcasts_D is VECTOR[SERIES[REAL]], and there are nsteps of them for EACH timestep ahead, EACH of size ndraws.

Using fcasts_D, I have plotted the densities, and the PI's within a cloud of simulations.


If I try in a multivariate setting (D_fcasts_1 is defined with GSET in the code from above):

declare vector[series] D_fcasts_1_vs(nsteps)
do ng = 1, nsteps
make D_array_fc_1 begin+ng begin+ng
# D_fcasts_1
set D_fcasts_1_vs(ng) 1 ndraws = D_array_fc_1(1,t)
end do ng


I get the following error

## SX22. Expected Type INTEGER, Got SERIES[VECTOR[REAL]] Instead
The Error Occurred At Location 119, Line 3 of loop/block

[TomDoan]

I have no idea what the MAKE instructions are trying to accomplish, but why don't you just do what works and repeat for two series rather than one?