The RATS Software Forum

Posted: **Thu Nov 02, 2017 9:33 am**

Dear Tom,

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

I safed all the GIRFs at the end of the code in the rectangular[series] and thats where I got stuck.

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

Posted: **Thu Nov 02, 2017 3:50 pm**

Have you looked at the Kilian-Vigfusson examples? Or any of the other GIRF examples? The bookkeeping on those is pretty similar---what's different is the method of generating the responses, not the method of collecting the information. There's no point in saving all the responses, since it's only the average that matters. So just add them up as you go,

Posted: **Thu Nov 02, 2017 4:48 pm**

Thank you for your reply.
Of course the average matters, but don't I have to save the draws to be able to calculate the average?

Can you recommend another GIRF example except Kilian and vigfusson or balke? Just to compare the way they do it?

Thank you for your great Support

Best

Jules

Posted: **Thu Nov 02, 2017 5:52 pm**

ONLY the average matters, which means only the sum matters. Look at the other GIRF code. In pseudocode

girf=0
do draws=1,ndraws
   girf=girf+(perturbed-base)
end do draws
girf=girf/ndraws

Posted: **Fri Nov 03, 2017 11:29 am**

Dear Tom,

for the above mentioned paper I tried to program the bootstrap for the GIRFs. The paper estimates an Interacted VAR (I-VAR), which analyses uncertainty shocks measured by the VIX during times of ZLB and during normal times.
The program below estimates the GIRFs, which should correspond to Steps 1 to 6 of the Appendix "Computation of the Generalized Impulse Response Funktions"
Step 7 is another bootstrap for confidence bands, but I left that out for a moment.

The code is the following, while the data are the same as above:

Code: Select all

***********************
*** Reading in Data ***
***********************

OPEN DATA "...\DataReplicationRats.xlsx"
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 = 500                                                          ;* #-of inner Bootstrap Draaws for GIRFs
compute nhist = 100                                                           ;* #-of outer draws for initial conditions
compute nsteps = 30                                                           ;* #-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

** reducedform
set u1 = u(1)
set u2 = u(2)
set u3 = u(3)
set u4 = u(4)
set u5 = u(5)
set u6 = u(6)

** structural
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 vect[series] girf(nvar)                                            ;*storage space for GIRF
declare vect[series] girfnormal(nvar)                                      ;*storgae space for GIRF during normal times
declare vect[series] girfzlb(nvar)                                         ;*storgae space for GIRF during ZLB Times

declare vect[strings] irflabels(neq)
compute irflabels =||"VIX","PCI","GDP","INV","CONS","FFR"||

******************
*** GIRF Inner ***
******************


infobox(action=define,progress,lower=1,upper=nvar) "Boot"
do i = 1, nvar
   infobox(current=i)

   set girf(i) hstart hend = 0.0
   set girfnormal(i) hstart hend = 0.0
   set girfzlb(i) hstart hend = 0.0

   compute nzlb = 0
   compute nnormal =0

   seed 123                                                          ;* To have for every variable the same bootstrap draw

   do hist = 1,nhist

      ********************************************************
      *** Draw block of length nlag for initial conditions ***
      ********************************************************

      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))

      if ientries(hstart-1)>= 182{
         compute nzlb = nzlb+1
      }
      else{
         compute nnormal = nnormal+1
      }

      do draw=1,ndraws

         *************************************************************
         *** Draw random entries for residuals for forecast period ***
         *************************************************************

         boot entries hstart hend bstart bend            
         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)     
         # 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

         ********************************************************************
         *** Compute GIRF for given draw and sum them up generate Average ***
         ********************************************************************

         set girf(i) hstart hend = girf(i) + (wshock(i)-base(i))

      end do draw

      set girf(i) hstart hend = girf(i)/ndraws

      if ientries(hstart-1)>= 182{                                              ;* 182 corresponds to 2008:4 where ZLB period starts, hstart-1 is the end period of the block bootstrap
         set girfzlb(i) hstart hend = girfzlb(i) + girf(i)
      }
      else{
         set girfnormal(i) hstart hend = girfnormal(i) + girf(i)
      }

   end do hist

   ****************************************************
   *** Average GIRF over ZLB State and Normal State ***
   ****************************************************

   set girfzlb(i) hstart hend = girfzlb(i)/nzlb
   set girfnormal(i) hstart hend = girfnormal(i)/nnormal

   display irflabels(i)
   display "ZLB Periods" nzlb
   display "normal Periods" nnormal


end do i
infobox(action=remove)




spgraph(vfields=2,hfields=3)
   do i = 1, nvar

      graph(number = 0,header = irflabels(i), key=upright) 2
      # girfzlb(i)
      # girfnormal(i)

   end do i
spgraph(done)

I have two questions:

1) Does the code look correct to you? Especially, I am concerned with the use of the forecast instruction for the I-VAR to generate E(y|d,e)-E(y|e). I had trouble with that because the interaction variable is an exogenous regressor, therefore I used an AR model to close it in the forecasting part.

2)
I simply used structural residuals to calulate the the GIRFs as in Kilian and Vigfusson (2011). The paper mentions that they are doing it similar as in Kilian and Vigfusson (2011), but in step 2 and step 3 of the appendix they use the reduced form residuals and pertubate them with a stuctural shock. I do not understand why there is the need for that, when instead structural residuals can be used from the beginning on.

My point estimates look different from those in the paper, but I don't know whether I did anything wrong related to question (1) or (2) or both.

It would be great I you could have a look at that

Thank you very much in advance

Best

Jules

Posted: **Fri Nov 03, 2017 12:21 pm**

What is it that you think is exogenous? Read their step 1. It's a self-contained model, which has a non-linear interaction term in it.

Posted: **Fri Nov 03, 2017 12:29 pm**

Thanks for your quick reply

The wording "exogenous" was wrong since the interaction term includes endogenously modeled variable. But what I mean by exogenous is the way the coefficients of the non-linear interaction term are estimated in the model:

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

The interaction term in "interact" does not have an own equation, but shows up in every other equation.

Do you have suggestions regarding the questions above? I would really appreciate you help, as I am stuck in the programming

Thanky in advance

Best Jules

Posted: **Fri Nov 03, 2017 12:53 pm**

It's not exogenous, and you're confusing yourself by using that terminology. The Kilian-Vigfusson model is similar, in that the equations are linear in things that are observable (and thus can be estimated very simply), but one of those observables (XPLUS) is a function of the other observables, thus making it non-linear out-of-sample. You need to add to the model an identity which maps the left-side variables to the non-linear function used on the right-side (the XPLUSDEF in the K-V example).

Posted: **Fri Nov 03, 2017 1:09 pm**

You are right, the wording is simply wrong...
So if I understand you correctly it is something like:

frml(identity) intereq interact =VIX*FFR

That raises three questions:

1) Whats the whole purpose of adding those identities? Is it that when I have that identity in the model and do the forecasting that the value for the interact variable is calculated automatically?

2) how do I add that to the system which I defined by

system(model=IVAR)
variables VIX lnP lnGDP lnInv lnCons FFR
lags 1 to nlag
det constant interact{1 to nlag}
end(system)

3) do I still need to estimate the AR model for the forecasting, which was defined by

linreg(define=intereq) interact
# interact{1 to nlag}

Thank you

Best Jules

Posted: **Fri Nov 03, 2017 1:50 pm**

Jules89 wrote:You are right, the wording is simply wrong...
So if I understand you correctly it is something like:

frml(identity) intereq interact =VIX*FFR

That raises three questions:

1) Whats the whole purpose of adding those identities? Is it that when I have that identity in the model and do the forecasting that the value for the interact variable is calculated automatically?

Yes.

Jules89 wrote: 2) how do I add that to the system which I defined by

system(model=IVAR)
variables VIX lnP lnGDP lnInv lnCons FFR
lags 1 to nlag
det constant interact{1 to nlag}
end(system)

Use IVAR+INTEREQ for the model

Jules89 wrote: 3) do I still need to estimate the AR model for the forecasting, which was defined by

linreg(define=intereq) interact
# interact{1 to nlag}

Absolutely not. You asked me how to handle an exogenous variable. It isn't exogenous. The proper handling of what it actually is is shown above.

Posted: **Mon Nov 06, 2017 6:09 am**

Dear Tom,

regarding adding an identity to the model, the instruction with frml(identity) intereq interact =VIX*FFR does not work in the code when I use IVAR+INTEREQ in the forecast instruction.
It produces the following error:

Can't Interpret MODEL + FRML[REAL]
## SX27. Illegal Combination of Data Types for Operation
>>>>model=IVAR+intereq,<<<<

So is there any possibility to add an identiy to a model defined by the system instruction, or do I have do estimate every equation of the VAR seperately and then group them together in a model?

Thank you in advance best

Jules

Posted: **Tue Nov 07, 2017 8:04 am**

I think I figured it out.
Instead of

Code: Select all

"frml(identity) intereq interact =VIX*FFR

I added

Code: Select all

equation(identity,coeffs=||1||) intereq interact
# VIX*FFR

Then the forecast instruction below does not produce an error:

Code: Select all

forecast(model=IVAR+intereq,from=hstart,to=hend,results=base,paths,noprint)   
# xpath(1) xpath(2) xpath(3) xpath(4) xpath(5) xpath(6)

Is this the correct way to add the interaction term identity to the model for the forecasting instruction?

Thank you

Best

Jules

Posted: **Tue Nov 07, 2017 1:53 pm**

No. VIX*FFR in the EQUATION is the product of the series handles, not the series. It looks like MODEL+FRML is on the wish list. You need to expand out the model on a GROUP instruction.

group fullmodel %modeleqn(ivar,1) %modeleqn(ivar,2) %modeleqn(ivar,3) %modeleqn(ivar,4) %modeleqn(ivar,5) %modeleqn(ivar,6) your_nonlinear_identity

and use fullmodel on the FORECAST.

Posted: **Mon Nov 13, 2017 8:20 am**

Thanks Tom it workes,

I am currently replicating the results for Caggiano et al. (2017). The paper is here http://www.sciencedirect.com/science/ar ... 2117301605
I also found the Matlab code of the paper on the personal page of one of the authors. See https://sites.google.com/site/giovannip ... e/research and https://onedrive.live.com/?authkey=!AAA ... ion=locate

Here is my code:

Code: Select all


***********************
*** Reading in Data ***
***********************

OPEN DATA "..."
CALENDAR(Q) 1962:3
DATA(FORMAT=XLSX,ORG=COLUMNS) 1962:03 2015: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 = 500                                                          ;* #-of inner Bootstrap Draaws for GIRFs
compute nsteps = 20+5                                                         ;* #-of GRIF steps
compute neq = 6                                                               ;* #-of equations in VAR

set time = t
set interact = VIX*FFR


frml(identity) intereq 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,ols) * 2015:04

group fullmodel %modeleqn(ivar,1) %modeleqn(ivar,2) %modeleqn(ivar,3) $
 %modeleqn(ivar,4) %modeleqn(ivar,5) %modeleqn(ivar,6) intereq


********************************************
*** 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

** reducedform
set u1 = u(1)
set u2 = u(2)
set u3 = u(3)
set u4 = u(4)
set u5 = u(5)
set u6 = u(6)

stats(noprint) u1 %regstart() 2008:3
compute numbnormal = %nobs

stats(noprint) u1 2008:4 %regend()
compute numbzlb =%nobs



** structural
set v1 = v(1)
set v2 = v(2)
set v3 = v(3)
set v4 = v(4)
set v5 = v(5)
set v6 = v(6)



*********************
*** Define Ranges ***
*********************


compute bstart=%regstart(),bend=%regend()                 ;* Range of residuals as source for Bootstrapping
compute normstart=%regstart(),normend=2008:3              ;* Range of residuals during normal period
compute zlbstart=normend+1,zlbend=%regend()               ;* Range of residuals during zlb period
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 extracting histories ***
****************************************************


set VIXraw dstart dend = VIX                             ;* Set up raw data to take the histories from
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 vect[series] girf(nvar)                                            ;*storage space for GIRF
declare vect[series] girfnormal(nvar)                                      ;*storgae space for GIRF during normal times
declare vect[series] girfzlb(nvar)                                         ;*storgae space for GIRF during ZLB Times
declare vect[series] xpath(nvar)
declare series[integer] hists

declare vect[strings] irflabels(neq)
compute irflabels =||"VIX","CPI","GDP","INV","CONS","FFR"||

compute chol = %decomp(%sigma)
compute invchol = inv(chol)


****************************************
*** All normal and all ZLB Histories ***
****************************************

**************
*** Normal ***
**************



do i = 1,nvar

   set girf(i) hstart hend = 0.0
   set girfnormal(i) hstart hend = 0.0

end do i
   

compute nhist = 0
do hist = normstart, normend


   do k=0, nlag-1

      compute VIX(hstart-nlag+k) = VIXraw(hist-nlag+k)          ;* Pulls the history and sets its values into the following range: Lag periods before hstart to one period before hstart, such that hstart is the first forecasting period
      compute lnP(hstart-nlag+k) = lnPraw(hist-nlag+k)
      compute lnGDP(hstart-nlag+k) = lnGDPraw(hist-nlag+k)
      compute lnINV(hstart-nlag+k) = lnINVraw(hist-nlag+k)
      compute lnCons(hstart-nlag+k) = lnConsraw(hist-nlag+k)
      compute FFR(hstart-nlag+k) = FFRraw(hist-nlag+k)
      compute interact(hstart-nlag+k) = interactraw(hist-nlag+k)

   end do k


   do draw=1,ndraws

      *************************************************************
      *** Draw random entries for residuals for forecast period ***
      *************************************************************

      boot entries hstart hend bstart bend            ;*draws residuals from the whole range and puts them in the range hstart to hend for loading the forecast

      do z=1,nvar

         set xpath(z) hstart hend = u(z)(entries(t))

      end do z


      ***********************************************************************************
      *** Simulate the path y_(t+h) by loading the VAR with the bootstraped residuals ***
      ***********************************************************************************

      forecast(model=fullmodel,from=hstart,to=hend,results=base,paths,noprint)     ;* Forecasts values from hstart to hend and loads the forecast with the shocks
      # xpath(1) xpath(2) xpath(3) xpath(4) xpath(5) xpath(6)

      ***************************************************************************************************************************
      *** Simulate the path y_(t+h) by loading the VAR with shocks where the first period is the structural uncertainty shock ***
      ***************************************************************************************************************************

      @structresids(factor=chol) xpath hstart hstart vimp
      *compute delta = 1/(chol(1,1))
      compute delta = 1

      compute vimp(1)(hstart) = vimp(1)(hstart)+delta
      @structresids(factor=invchol) vimp hstart hstart xpathpert

      do j =1, nvar
         compute xpath(j)(hstart) = xpathpert(j)(hstart)
      end do j

      forecast(model=fullmodel,from=hstart,to=hend,results=wshock,paths,noprint)
      # xpath(1) xpath(2) xpath(3) xpath(4) xpath(5) xpath(6)

      ********************************************************************
      *** Compute GIRF for given draw and sum them up generate Average ***
      ********************************************************************
      

      do i = 1, nvar

         set girf(i) hstart hend = girf(i) + (wshock(i)-base(i))

      end do i

   end do draw

      
   do i = 1, nvar
      set girf(i) hstart hend = girf(i)/ndraws
      set girfnormal(i) hstart hend = girfnormal(i) + girf(i)
   end do i

   compute nhist = nhist+1

end do hist


do i = 1,nvar

   set girfnormal(i) hstart hend = girfnormal(i)/nhist

end do i



***********
*** ZLB ***
***********


do i = 1,nvar

   set girf(i) hstart hend = 0.0
   set girfzlb(i) hstart hend = 0.0

end do i
   

compute nhist = 0
do hist = zlbstart, zlbend


   do k=0, nlag-1

      compute VIX(hstart-nlag+k) = VIXraw(hist-nlag+k)          ;* Pulls the history and sets its values into the following range: Lag periods before hstart to one period before hstart, such that hstart is the first forecasting period
      compute lnP(hstart-nlag+k) = lnPraw(hist-nlag+k)
      compute lnGDP(hstart-nlag+k) = lnGDPraw(hist-nlag+k)
      compute lnINV(hstart-nlag+k) = lnINVraw(hist-nlag+k)
      compute lnCons(hstart-nlag+k) = lnConsraw(hist-nlag+k)
      compute FFR(hstart-nlag+k) = FFRraw(hist-nlag+k)
      compute interact(hstart-nlag+k) = interactraw(hist-nlag+k)

   end do k


   do draw=1,ndraws

      *************************************************************
      *** Draw random entries for residuals for forecast period ***
      *************************************************************

      boot entries hstart hend bstart bend            ;*draws residuals from the whole range and puts them in the range hstart to hend for loading the forecast

      do z=1,nvar

         set xpath(z) hstart hend = u(z)(entries(t))

      end do z


      ***********************************************************************************
      *** Simulate the path y_(t+h) by loading the VAR with the bootstraped residuals ***
      ***********************************************************************************

      forecast(model=fullmodel,from=hstart,to=hend,results=base,paths,noprint)     ;* Forecasts values from hstart to hend and loads the forecast with the shocks
      # xpath(1) xpath(2) xpath(3) xpath(4) xpath(5) xpath(6)

      ***************************************************************************************************************************
      *** Simulate the path y_(t+h) by loading the VAR with shocks where the first period is the structural uncertainty shock ***
      ***************************************************************************************************************************

      @structresids(factor=chol) xpath hstart hstart vimp
      *compute delta = 1/(chol(1,1))
      compute delta = 1

      compute vimp(1)(hstart) = vimp(1)(hstart)+delta
      @structresids(factor=invchol) vimp hstart hstart xpathpert

      do j =1, nvar
         compute xpath(j)(hstart) = xpathpert(j)(hstart)
      end do j

      forecast(model=fullmodel,from=hstart,to=hend,results=wshock,paths,noprint)
      # xpath(1) xpath(2) xpath(3) xpath(4) xpath(5) xpath(6)

      ********************************************************************
      *** Compute GIRF for given draw and sum them up generate Average ***
      ********************************************************************
      
      do i = 1, nvar

         set girf(i) hstart hend = girf(i) + (wshock(i)-base(i))

      end do i

   end do draw

      
   do i = 1, nvar
      set girf(i) hstart hend = girf(i)/ndraws
      set girfzlb(i) hstart hend = girfzlb(i) + girf(i)
   end do i

   compute nhist = nhist+1

end do hist


do i = 1,nvar

   set girfzlb(i) hstart hend = girfzlb(i)/nhist

end do i




spgraph(vfields=2,hfields=3)
   do i = 1, nvar


      graph(number = 0,header = irflabels(i), key=loright) 2
      # girfnormal(i) hstart hend-5
      # girfzlb(i) hstart hend-5


   end do i
spgraph(done)

The code estimates the GIRFs without the bootstrap, that crates the confidence bands.
I can perfectly replicate the GIRFs for the normal times (Figure 1 in the paper). But for the ZLB times I get the following error message:

## FO7. Gauss-Seidel solution is explosive. Check model. Try DAMP option.
The Error Occurred At Location 930, Line 44 of loop/block
3725856 Position 7137

I Think it has something to do with explosive draws. However when I add for example "seed 12" before all the bootstrapping I do not get the problem. Therefor it depends on some explosive draws. Is there any possibility to discard explosive draws and simly jump to the next iteration? The authors write that they discard explosive draws.

Thank you

Best

Jules

Posted: **Mon Nov 13, 2017 11:29 am**

It has nothing to do with explosive "draws". The model itself has potentially explosive non-linear behavior. Starting with the initial conditions for 2009:04, FFR goes immediately negative, which causes all kinds of bad things to happen. You would have to look carefully through the paper to see what they do to prevent FFR from going negative during simulations.

The RATS Software Forum

GIRFs for I-VAR

GIRFs for I-VAR

Re: GIRFs for I-VAR

Re: GIRFs for I-VAR

Re: GIRFs for I-VAR

Re: GIRFs for I-VAR

Re: GIRFs for I-VAR

Re: GIRFs for I-VAR

Re: GIRFs for I-VAR

Re: GIRFs for I-VAR

Re: GIRFs for I-VAR

Re: GIRFs for I-VAR

Re: GIRFs for I-VAR

Re: GIRFs for I-VAR

Re: GIRFs for I-VAR

Re: GIRFs for I-VAR