The error bands for forecasts are relatively simple if you act as if the estimated coefficients are known rather than subject to sampling error.
FORECAST with the
STDERRS option takes care of most of the work. A couple of fairly simple examples are the Lutkepohl textbook example lutkp098.rpf and the Tsay textbook example tsayp407.rpf.
If you would like to take the uncertainty in the coefficients into account, there are, as with error bands for impulse responses, three basic ways to handle it: by Monte Carlo simulation methods, bootstrapping and the delta method. I don't have any examples applying the last two methods to forecasting. Bootstrapping requires a relatively simple change to a bootstrapping loop for impulse responses. The delta method is quite a bit harder--Lutkepohl derives the formula for the forecast variance in his book but doesn't have an application.
The following is an example of simulation methods, applied to a BVAR. This is an example from the
Bayesian Econometrics e-course. See
http://www.estima.com/courses_completed.shtml for details.
Code: Select all
*
* Bayesian Econometrics Workbook, Example 6.8
* Gibbs sampling on a Bayesian VAR with a "Minnesota" prior
*
compute lags=4 ;*Number of lags
compute nstep=16 ;*Number of response steps
compute nburn=500 ;*Number of burn-in draws
compute ndraws=2500 ;*Number of keeper draws
*
open data haversample.rat
cal(q) 1959
data(format=rats) 1959:1 2006:4 ftb3 gdph ih cbhm
*
set loggdp = log(gdph)
set loginv = log(ih)
set logc = log(cbhm)
*
* These are the controlling parameters for the symmetric "Minnesota"
* prior - the own lag tightness and the relative tightness on the other
* lags.
*
compute tight=.1
compute other=.5
*
* First, estimate the system with a prior via mixed estimation, to see if
* we come up with something similar via Gibbs sampling.
*
system(model=varmodel)
variables loggdp loginv logc ftb3
lags 1 to lags
specify(type=symmetric,tight=tight) other
det constant
end(system)
************************************************************************************
*
* Estimate with a prior
*
estimate
compute nvar=%nvar
*
source bvarbuildprior.src
@BVARBuildPriorMN(model=varmodel,tight=.1,other=.5) hbpriors hpriors
@BVARFinishPrior hbpriors hpriors bprior hprior
*
compute sigmad=%sigma
*
dec vect[series] forecast(nvar)
dec vect[series] forestderr(nvar)
*
* Range to forecast
*
compute fstart=%regend()+1
compute fend =fstart+nstep-1
*
dec series[rect] forecasts
gset forecasts fstart fend = %zeros(nvar,ndraws)
*
infobox(action=define,progress,lower=-nburn,upper=ndraws) "Gibbs Sampler"
do draw=-nburn,ndraws
infobox(current=draw)
*
* Draw b given sigma
*
compute bdraw =%ranmvkroncmom(%cmom,inv(sigmad),hprior,bprior)
compute rssmat=%sigmacmom(%cmom,bdraw)
*
* Draw sigma given b
*
compute sigmad=%ranwisharti(%decomp(inv(rssmat)),%nobs)
if draw<=0
next
compute %modelsetcoeffs(varmodel,bdraw)
simulate(model=varmodel,cv=sigmad,results=simresults,steps=nstep)
gset forecasts fstart fend = %psubmat(forecasts(t),1,draw,%xt(simresults,t))
end do draw
infobox(action=remove)
*
dec vect[series] faninput(7)
do i=1,7
set faninput(i) fstart fend = 0.0
end do i
*
dec vect fjt percent
do j=1,nvar
do time=fstart,fend
compute fjt=%xrow(forecasts(time),j)
compute percent=%fractiles(fjt,||.005,.05,.25,.50,.75,.95,.995||)
compute %pt(faninput,time,percent)
end do time
graph(style=line,overlay=fan,ovcount=6,ovsame,$
header="Forecasts of "+%l(%modeldepvars(varmodel)(j))) 7
# faninput(4)
# faninput(1)
# faninput(2)
# faninput(3)
# faninput(5)
# faninput(6)
# faninput(7)
end do j