' bootstrap structural VAR impulse responses
' from long-run restriction (Blanchard-Quah 1989)
' for version 4
' 11/28/2000 h

include sub_bootirf.prg

' create workfile
workfile blanquah q 1948:1 1987:4

' fetch data from database (already demeaned)
fetch(d="i:\hiroyuki\example files\var\data_svar") dy u

' change to percentage
dy = 100*dy

' estimate VAR (no constant)
var var1.ls 1 8 dy u @ 

' impose long-run restrictions
' @u1 = aggregate demand shock
' @u2 = aggregate supply shock
var1.cleartext(svar)
var1.append(svar) @lr1(@u1)=0

' estimate factorization
var1.svar(rtype=text,conv=1e-5)

' set response horizon
!hrz = 40

' store level and accumulated responses
do var1.impulse(!hrz,imp=struct,matbys=rsp_lvl) dy u @ 1
do var1.impulse(!hrz,imp=struct,a,matbys=rsp_acc) dy u @ 1

' get residuals to bootstrap
var1.makeresid(n=gres) res_dy res_u

' make model to generate bootstrap data
var1.makemodel(mod1)
' assign add factors for bootstrap residuals
mod1.addassign(i) @all

'set random number generator
rndseed(type=mt) 1234567
!reps = 100
'declare storage matrices
for !i=1 to @columns(rsp_lvl)
	matrix(!hrz,!reps) brsp_lvl{!i}
	matrix(!hrz,!reps) brsp_acc{!i}
next
'temporary working vector
vector vtmp
'count number of cases that converged
scalar cnt = 0

'bootstrap loop
for !i=1 to !reps
	'bootstrap residuals
	smpl @all if res_dy<>na
	gres.resample dy_a u_a
	'generate bootstrap data
	mod1.solve
	'estimate VAR (no constant) with bootstrap data
	smpl @all
	var var2.ls 1 8 dy_0 u_0 @ 
	'impose long-run restrictions
	var2.cleartext(svar)
	var2.append(svar) @lr1(@u1)=0
	'estimate factorization
	var2.svar(rtype=text,conv=1e-5,maxiter=100,nostop)
	'store results only if converged
	if (var2.@svarcvgtype=0) then
		'increment counter
		cnt = cnt + 1
		'get level and accumulated bootstrap responses
		do var2.impulse(!hrz,imp=struct,matbys=tmp_lvl) dy_0 u_0 @ 1
		do var2.impulse(!hrz,imp=struct,a,matbys=tmp_acc) dy_0 u_0 @ 1
		'and store in appropriate position
		for !j=1 to @columns(rsp_lvl)
			vtmp = @columnextract(tmp_lvl,!j)
			colplace(brsp_lvl{!j}, vtmp, cnt)
			vtmp = @columnextract(tmp_acc,!j)
			colplace(brsp_acc{!j}, vtmp, cnt)
		next
		delete tmp_lvl tmp_acc
	endif
next

'reshape results matrix if there were non-converged cases
if (cnt <> !reps) then
	for !i=1 to @columns(rsp_lvl)
		brsp_lvl{!i} = @subextract(brsp_lvl{!i}, 1, 1, !hrz, cnt)
		brsp_acc{!i} = @subextract(brsp_acc{!i}, 1, 1, !hrz, cnt)
	next
	table warning
	%msg = @str(!reps-cnt) + " case(s) out of " + @str(!reps) + " replications failed to converge"
	warning(1,1) = %msg
	show warning
endif

'find bootstrap percentiles for each horizon
matrix(!hrz,@columns(rsp_lvl)) lvl_upp
matrix(!hrz,@columns(rsp_lvl)) lvl_low
matrix(!hrz,@columns(rsp_lvl)) acc_upp
matrix(!hrz,@columns(rsp_lvl)) acc_low

rowvector rtmp
for !i=1 to !hrz
	for !j=1 to @columns(rsp_lvl)
		rtmp = @rowextract(brsp_lvl{!j},!i)
		lvl_upp(!i,!j) = @quantile(rtmp,0.833)
		lvl_low(!i,!j) = @quantile(rtmp,0.167)
		rtmp = @rowextract(brsp_acc{!j},!i)
		acc_upp(!i,!j) = @quantile(rtmp,0.833)
		acc_low(!i,!j) = @quantile(rtmp,0.167)
	next
next

'plot impulse responses with bootstrap error bounds
%gname = "fig_lvl"
call sub_bootirf(rsp_lvl, lvl_upp, lvl_low, %gname)
' add zero line to all graphs
{%gname}.draw(line,left,rgb(180,180,180)) 0
' add title to graph
{%gname}.addtext(0.1,-0.2) Bootstrap 66% Confidence Bounds (Levels)
' realign
{%gname}.align(2,1,0.5)
show {%gname}

%gname = "fig_acc"
call sub_bootirf(rsp_acc, acc_upp, acc_low, %gname)
' add zero line to all graphs
{%gname}.draw(line,left,rgb(180,180,180)) 0
' add title to graph
{%gname}.addtext(0.1,-0.2) Bootstrap 66% Confidence Bounds (Accumulated Responses)
' realign
{%gname}.align(2,1,0.5)
show {%gname}
