*
* Example of forecasts (no worked example in text)
* Based upon discussion on pp 578-9
*
open data tablef5-1[1].txt
calendar(q) 1950
all 2002:04
data(format=prn,org=columns) 1950:1 2000:4 year qtr realgdp realcons realinvs realgovt realdpi $
  cpi_u m1 tbilrate unemp pop infl realint
*
set logc = log(realcons)
set logy = log(realgdp)
*
* This creates a dummy variable for quarter 4. This series with leads 0 to 3
* (that is regressors seas{0 to -3}) gives you a full set of dummies.
*
seasonal seas
*
linreg(define=ardl) logc / resids
# logc{1 to 3} logy{0 to 3} seas{0 to -3}
*
* Two years worth of out-of-sample forecasts. We need an input path for the
* "exogenous" variable logy. We'll give it 3% annual growth. Since we're in logs,
* we need to increment by .0075 per period.
*
set logy 2001:1 2002:4 = logy{1}+.0075
*
* We're going to define some variables so it will be easier to modify the samples
* below.
*
compute regStart=%regStart(),regEnd=2000:4
compute foreStart=2001:1,foreEnd=2002:4,nfore=8
*
* Compute the base forecasts and the standard errors ignoring sampling error in
* the coefficient estimates
*
uforecast(equation=ardl,stderrs=basestderr,steps=nfore) basefore
*
* Bootstrapped forecast errors
* Get a copy of the estimated equation with logc replaced by the new variable
* "bootc". We'll use this for creating the bootstrapped data.
*
set bootc * 2000:4 = logc
modify ardl bootdl
vreplace logc with bootc
*
* This is the number of bootstrap draws
*
compute ndraws=1000
*
* We're going to do two sets of bootstraps. One will involve re-estimating the
* model with bootstrapped residuals throughout the sample, and then bootstrapped
* errors during the forecast period as well. This will give an estimate of the
* standard error of forecast taking into account everything but the uncertainty
* about the future exogenous variable. The other will use the original equation
* and will bootstrap only the out of sample errors. This should give us a result
* similar to the "basestderr" series generated above.
*
set ssq1 foreStart foreEnd = 0.0
set ssq2 foreStart foreEnd = 0.0
set boot1 foreStart foreEnd = 0.0
set boot2 foreStart foreEnd = 0.0
do draws=1,ndraws
*
*    This draws with replacement a set of entry numbers between regStart and regEnd
*    (the range for which we have residuals), one for each entry from regStart to
*    foreEnd (the range for which we need draws for the errors).
*
   boot entries regStart foreEnd regStart regEnd
*
*    And this takes the entries and converts them into a vector of resampled
*    residuals
*
   set reDraw regStart foreEnd = resids(entries(t))
*
*    FORECAST with PATHS generates the bootstrapped version of the logc series, by
*    simulating the model with the resampled shocks added as the errors. Note that
*    this is conditional on the entire y series and on the pre-regression sample
*    values of logc, that is, the initial lagged values of logc aren't altered. The
*    equation used is the original estimated equation modified only to work with
*    bootc and not logc.
*
   forecast(paths,from=regStart,to=regEnd)
   # bootdl bootc
   # reDraw
*
*    Reestimate the model with the bootstrapped data. Define a new equation from
*    this, since we need to keep bootdl as it was for bootstrapping the sample data
*    on other draws.
*
   linreg(equation=bootdl,define=workdl,noprint) bootc regStart regEnd
*
*    Now that the estimation is finished, we need to replace the bootstrapped data
*    with the original for doing the forecasts.
*
   set bootc regStart regEnd = logc
*
*    This does the first set of bootstrapped forecasts, using the model just
*    estimated, with the resampled residuals for the forecast period added in.
*
   forecast(paths,steps=nfore)
   # workdl boot1
   # reDraw
*
*    And this does the second set, which uses the original equation, adding in only
*    the resampled shocks.
*
   forecast(paths,steps=nfore)
   # ardl boot2
   # reDraw
*
*    Add up the squared differences between the bootstrapped forecast and the base
*    forecast
*
   set ssq1 foreStart foreEnd = ssq1+(boot1-basefore)**2
   set ssq2 foreStart foreEnd = ssq2+(boot2-basefore)**2
end do draws
*
* Convert the sums of squared differences to estimated standard errors by
* dividing by the number of draws to get the variance, then taking the square
* root.
*
set stderr1 foreStart foreEnd = sqrt(ssq1/ndraws)
set stderr2 foreStart foreEnd = sqrt(ssq2/ndraws)
print foreStart foreEnd basestderr stderr1 stderr2
*
* Method Two: Monte Carlo integration. Rather than drawing a new set of data and
* re-estimating the model, this draws a random set of coefficients from the
* posterior distribution. (?? take care of the sigma)
*
linreg(define=ardl) logc / resids
# logc{1 to 3} logy{0 to 3} seas{0 to -3}
compute xxbase=%xx,betabase=%beta
compute sxx=sqrt(%seesq)*%decomp(xxbase)
*
modify ardl montedl
set ssq1 foreStart foreEnd = 0.0
set ssq2 foreStart foreEnd = 0.0
set boot1 foreStart foreEnd = 0.0
set boot2 foreStart foreEnd = 0.0
do draws=1,ndraws
   compute beta=betabase+%ranmvnormal(sxx)
   compute %eqnsetcoeffs(montedl,beta)
   simulate(steps=nfore)
   # montedl boot1

   simulate(steps=nfore)
   # ardl boot2
*
*    Add up the squared differences between the Monte Carlo forecast and the base
*    forecast
*
   set ssq1 foreStart foreEnd = ssq1+(boot1-basefore)**2
   set ssq2 foreStart foreEnd = ssq2+(boot2-basefore)**2
end do draws
set stderr1 foreStart foreEnd = sqrt(ssq1/ndraws)
set stderr2 foreStart foreEnd = sqrt(ssq2/ndraws)
print foreStart foreEnd basestderr stderr1 stderr2