'bootstrap VAR forecasts
'with nonlinear transformation (log difference)
'12/15/2000 h

'load workfile
%data = @runpath + "lut1"
load "{%data}"

'rename original series
y1 = inv
y2 = inc
y3 = cs

'set sample objects
sample smpl_est	 1960:1 1978:4	'estimation sample
sample smpl_fcst 1979:1 1980:4	'forecast sample

'estimate VAR
smpl smpl_est
var1.ls 1 2 dlog(y1) dlog(y2) dlog(y3) @ c

'get residuals to bootstrap
var1.makeresid(n=gres) res1 res2 res3

'make model out of estimated VAR
var1.makemodel(mod1)

'solve model to obtain dynamic forecasts
smpl smpl_fcst
mod1.solve
'rename forecasts
smpl smpl_est
for !i=1 to var1.@neqn
	rename y{!i}_0 y{!i}f
	y{!i}f = na
next

'assign add factors for bootstrap residuals
mod1.addassign(i) @all

'set monte carlo parameters
!reps = 500						'bootstrap replications
!hrz = smpl_fcst.@obs			'forecast horizon

'allocate storage matrix for entire forecast distribution
for !i=1 to var1.@neqn
	matrix(!hrz,!reps) mboot{!i}
next
vector vtmp						'temporary working vector

'set random number generator
rndseed(type=mt) 1234567

'bootstrap loop
for !i=1 to !reps
	'bootstrap residuals
	smpl @all if res1<>na
	gres.resample y1_a y2_a y3_a
	'generate bootstrap data
	mod1.solve
	'estimate VAR with bootstrap data
	smpl smpl_est
	var var2.ls 1 2 dlog(y1_0) dlog(y2_0) dlog(y3_0) @ c
	'make model from bootstrap estimates
	var2.makemodel(mod2)
	'solve to bootstrap forecasts
	smpl smpl_fcst
	mod2.solve
	'store bootstrap forecasts
	for !j=1 to var2.@neqn
		stom(y{!j}_0_0, vtmp, smpl_fcst)
		colplace(mboot{!j}, vtmp, !i)
	next
	'clean up
	delete mod2
next

'find bootstrap percentiles for each horizon
matrix(!hrz,var1.@neqn) mupp	'upper percentile
matrix(!hrz,var1.@neqn) mlow	'lower percentile
rowvector rtmp					'temporary working vector

for !i=1 to !hrz
	for !j=1 to var1.@neqn
		rtmp = @rowextract(mboot{!j},!i)
		mupp(!i,!j) = @quantile(rtmp,0.975)
		mlow(!i,!j) = @quantile(rtmp,0.025)
	next
next

'plot actual and forecasts
smpl 1978:1 1980:4
for !i=1 to var1.@neqn
	'convert bounds into series
	series y{!i}_upp
	vtmp = @columnextract(mupp,!i)
	mtos(vtmp, y{!i}_upp, smpl_fcst)

	series y{!i}_low
	vtmp = @columnextract(mlow,!i)
	mtos(vtmp, y{!i}_low, smpl_fcst)
	'plot graph	
	group gtmp y{!i} y{!i}f y{!i}_upp y{!i}_low
	freeze(gra{!i}) gtmp.line
	gra{!i}.setelem(1) legend(actual)
	gra{!i}.setelem(2) legend(forecast)
	gra{!i}.setelem(3) legend(bootstrap bounds)
	gra{!i}.setelem(4) legend()
	
	%gname = %gname + "gra" + @str(!i) + " "
next

' merge all graphs into one
freeze(gfcst) {%gname}
gfcst.options size(3,4)
gfcst.align(3, 0.7, 0)
gfcst.legend columns(1) -inbox position(0,0)
'gfcst.scale(left) +zeroline
'gfcst.draw( dashline,left,rgb(155,155,155) ) 0.0
gfcst.option linepat	' need to set linepat
gfcst.elem(1) lcolor(blue) lpat(solid)
gfcst.elem(2) lcolor(red) lpat(solid)
gfcst.elem(3) lcolor(red) lpat(dash1)
gfcst.elem(4) lcolor(red) lpat(dash1)
gfcst.addtext(0.1,-0.2) Bootstrap Forecast Bounds

show gfcst

