I am currently writing a programm to calculate GIRFs form a non-linear VAR based on http://economia.unipv.it/eco-pol/mbf201 ... mating.pdf .
Within the Bootstrap I calculate 100 draws of my GIRFs for six series. The draws are safed in a rect[series] girf(6,100). Each of the GIRFs has 20 steps. So GIRF(1,1)(2) would for example be the response of variable 1, of draw 1, for step 2 .
Now I want to do for each of the GIRFs montecarlo integration by calculating the horizon wise average. That is I want to calculate:
(GIRF(1,1)(1)+GIRF(1,2)(1)+...+GIRF(1,100)(1))/100 = ghat1(1)
(GIRF(1,1)(2)+GIRF(1,2)(2)+...+GIRF(1,100)(2))/100 = ghat1(2)
...
(GIRF(1,1)(20)+GIRF(1,2)(20)+...+GIRF(1,100)(20))/100 = ghat1(20)
To get for the first variable the GIRF in a vector (ghat1(1), ghat1(2), ..., ghat1(20))
My code (so far) is:
Code: Select all
***********************
*** Reading in Data ***
***********************
OPEN DATA "...."
CALENDAR(Q) 1962:7
DATA(FORMAT=XLSX,ORG=COLUMNS) 1963:03 2016:04 VIX lnP lnGDP lnInv lnCons FFR
**********************************************************
************ I-VAR with exogenous regressors ************
**********************************************************
compute nvar = 6 ;* #-of VAR equations
compute nlag = 3 ;* #-of VAR Lags
compute ndraws = 100 ;* #-of inner Bootstrap Draaws for GIRFs
compute nhist = 100 ;* #-of outer draws for initial conditions
compute nsteps = 20 ;* #-of GRIF steps
compute neq = 6 ;* #-of equations in VAR
set interact = VIX*FFR
system(model=IVAR)
variables VIX lnP lnGDP lnInv lnCons FFR
lags 1 to nlag
det constant interact{1 to nlag}
end(system)
estimate(noprint) * 2015:04
linreg(define=intereq) interact
# interact{1 to nlag}
********************************************
*** Generate Series of Structural Shocks ***
********************************************
forecast(model=IVAR,static,from=%regstart(),to=%regend(),errors=u)
@structresids(factor=%decomp(%sigma)) u %regstart() %regend() v ;* Generate sturctural residuals from Cholesky Decomposition
set v1 = v(1)
set v2 = v(2)
set v3 = v(3)
set v4 = v(4)
set v5 = v(5)
set v6 = v(6)
stats(noprint) v1
compute delta = sqrt(%variance) ;* Generate one-standard deviation uncertainty shock
*********************
*** Define Ranges ***
*********************
compute bstart=%regstart(),bend=%regend() ;* Range of residuals as source for bootstrapping
compute hstart=%regend()+1,hend=%regend()+nsteps ;* Range of simulated forecasts
compute dstart=%regstart()-nlag,dend=%regend() ;* Range of data for blocks of pre-sample values
*********************************************
*** Make copies of data for bootstrapping ***
*********************************************
set VIXraw dstart dend = VIX
set lnPraw dstart dend = lnP
set lnGDPraw dstart dend = lnGDP
set lnINVraw dstart dend = lnInv
set lnConsraw dstart dend = lnCons
set FFRraw dstart dend = FFR
set interactraw dstart dend = interact
************************************
*** Setup storage space for IRFs ***
************************************
declare rect[series] girf(neq,ndraws)
declare vect[strings] irflabels(neq)
compute irflabels =||"VIX","PCI","GDP","INV","CONS","FFR"||
declare rect[series] girf(neq,ndraws)
***********************************************
*** Bootstrap for one initial condition ***
***********************************************
boot(block=nlag) ientries hstart-nlag hstart-1 dstart dend
set VIX hstart-nlag hstart-1 = VIXraw(ientries(t))
set lnP hstart-nlag hstart-1 = lnPraw(ientries(t))
set lnGDP hstart-nlag hstart-1 = lnGDPraw(ientries(t))
set lnInv hstart-nlag hstart-1 = lnInvraw(ientries(t))
set lnCons hstart-nlag hstart-1 = lnConsRaw(ientries(t))
set FFR hstart-nlag hstart-1 = FFRraw(ientries(t))
set interact hstart-nlag hstart-1 = interactraw(ientries(t))
do draw=1,ndraws
*************************************************************
*** Draw random entries for residuals for forecast period ***
*************************************************************
boot entries hstart hend bstart bend ;*draws shocks from the range bstart to bend and puts them in the range hstart to hend
set x1path hstart hend = v1(entries(t))
set x2path hstart hend = v2(entries(t))
set x3path hstart hend = v3(entries(t))
set x4path hstart hend = v4(entries(t))
set x5path hstart hend = v5(entries(t))
set x6path hstart hend = v6(entries(t))
****************************************************************************************
*** Simulate the path y_(t+h) by randomly loading the VAR with the bootstrapt shocks ***
****************************************************************************************
forecast(model=IVAR+intereq,from=hstart,to=hend,results=base,paths,noprint) ;* Forecasts values from hstart to hend and loads the forecast with the shocks
# x1path x2path x3path x4path x5path x6path
****************************************************************************************************************
*** Simulate the path y_(t+h) by loading the VAR with shocks where the first period is the uncertainty shock ***
****************************************************************************************************************
set x1path hstart hend = %if(t==hstart,delta,v1(entries(t)))
forecast(model=IVAR+intereq,from=hstart,to=hend,results=wshock,paths,noprint)
# x1path x2path x3path x4path x5path x6path
do i= 1,6
set girf(i,draw) hstart hend = wshock(i)-base(i)
end do i
end do draw
So far I have just coded that up for one random history of data. in what follows I will do this in a loop for more histories, similar to the Kilian and Vigfusson (2011) replication file.
I would really apprecite your help
Best Jules