Page 2 of 2
Re: GIRFs for I-VAR
Posted: Tue Nov 14, 2017 6:54 am
by Jules89
After going deeper into their Matlab code, I found the method they apply to discard GIRFs which are explosive.
They "repeat a draw if the residuals extracted make the GIRF explosive". So a GIRF is considered explosive if the last step ahead of the considered shock variable takes values which are extreme. In their example it is 10 standard deviations of the structural shock. I included this into the ZLB GIRF code. but I also made some additions/changes. I also discard values which make the first step ahead of the shock variable extremely large. My threshold value for this is 5 standrard deviations. Moreover since the Gauss-Seidel-Algorithm within the FORECAST instruction produces an error whenever there is an explosive draw, which would stop the entire loop, I use ENV TRAPERRORS before the Bootstrap. This should prevent the loop from breaking. After the loop continues the draw is repeated. Here is my code for the ZLB part:
Code: Select all
***********
*** ZLB ***
***********
do i = 1,nvar
set girf(i) hstart hend = 0.0
set girfzlb(i) hstart hend = 0.0
end do i
env traperrors
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
compute countrep = 0 ;* to count for repetitions in the draws
do draw=1,ndraws
*************************************************************
*** Draw random entries for residuals for forecast period ***
*************************************************************
:newtry
if countrep == ndraws{
compute nhist=nhist-1
break
}
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)
*********************************
*** Check for explosive IRFs ***
*********************************
compute nexpl = 0
compute test1 = wshock(1)(hend)-base(1)(hend)
compute test2 = wshock(1)(hstart)-base(1)(hstart)
if test1>=chol(1,1)*10*abs(delta).OR.test1<=(-1)*chol(1,1)*10*abs(delta).OR.test2>=chol(1,1)*5*abs(delta)
compute nexpl = nexpl+1
if nexpl>0{
compute countrep = countrep + 1
goto newtry
}
******************************************************
*** Accumulate GIRFs to calculate average later on ***
******************************************************
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
env notraperrors
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-1
# girfzlb(i) hstart hend-1
end do i
spgraph(done)
Together with the code for the normal part of the sample this produces exactly the mean GIRFS of figure 1 of the paper for any seed number I have tried so far.
Whats you opinion about the code? Do you have any further suggestions?
In general is there another way to force the loop to continue after an error occured, different from ENV TTRAPERRORS?
Thank you
best Jules
Re: GIRFs for I-VAR
Posted: Tue Nov 14, 2017 10:53 am
by TomDoan
I would be much happier if something like that were actually addressed in the paper rather than buried in the Matlab code.
It takes a lot for RATS to consider the forecast solution to be "explosive"---one of the simulated series has to have an absolute value larger than 1.e+100, which is hard to do in 25 steps. The fact that you're hitting that shows how unstable the model itself is when the FF goes negative. I'm surprised that they don't deal directly with that and reject draws which give a negative FF.
At any rate, yes, you have to do the TRAPERRORS to make it continue through the explosive solution error.
Re: GIRFs for I-VAR
Posted: Wed Nov 15, 2017 3:19 am
by Jules89
Sorry to bother you again with this, but ENV TRAPERRORS does not work as I implemented it above... Maybe I used it wrong. I added a progressbar to the loop and this clearly shows that after the error occured the loop breaks and does not go on:
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)
***********
*** ZLB ***
***********
do i = 1,nvar
set girf(i) hstart hend = 0.0
set girfzlb(i) hstart hend = 0.0
end do i
env traperrors
seed 1234
compute nhist = 0
infobox(action=define,progress,lower=zlbstart,upper=zlbend) "Boot zlb"
do hist = zlbstart, zlbend
infobox(current=hist)
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
compute countrep = 0 ;* to count for repetitions in the draws
do draw=1,ndraws
*************************************************************
*** Draw random entries for residuals for forecast period ***
*************************************************************
:newtry
if countrep == ndraws{
compute nhist=nhist-1
break
}
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)
*********************************
*** Check for explosive IRFs ***
*********************************
compute nexpl = 0
compute test1 = wshock(1)(hend)-base(1)(hend)
compute test2 = wshock(1)(hstart)-base(1)(hstart)
if test1>=chol(1,1)*10*abs(delta).OR.test1<=(-1)*chol(1,1)*10*abs(delta).OR.test2>=chol(1,1)*5*abs(delta)
compute nexpl = nexpl+1
if nexpl>0{
compute countrep = countrep + 1
goto newtry
}
******************************************************
*** Accumulate GIRFs to calculate average later on ***
******************************************************
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
infobox(action=remove)
env notraperrors
The seed 1234 generates a set of bootstrap draws that make the model explosive and therefore generates the error
"## FO7. Gauss-Seidel solution is explosive. Check model. Try DAMP option.
The Error Occurred At Location 1256, Line 63 of loop/block
36362784 Position 8220
" after the error occurs for the first time the loop breaks and does not continue to go on. The above code works fine for some values of SEED, but produces for others the explosive draws. Since I want to derive errorbands by bootstrapping within the next step, I really need this loop to work for any number of random draws. I would like to simply discard the explosive GIRFs and then rerun the iteration for another set of simulated residuals, but this requires RATS to be forced to continue after the above error occurs.
As far as I understand for now, the ENV TRAPERRORS just suppresses the print of the error message, but does not force the loop to continue. How would I use the ENV TRAPERRORS correctly? Or is there any other method to force the loop to continue after an error occurs?
Best
Jules
Re: GIRFs for I-VAR
Posted: Thu Nov 16, 2017 6:05 am
by Jules89
Dear Tom,
regarding my question above: Are there any other methods of forcing the loop to continue whenever errors occur in later versions of RATS? I am working on 9.0. Maybe it would be worth updating
Best Jules
Re: GIRFs for I-VAR
Posted: Thu Nov 16, 2017 12:14 pm
by TomDoan
You can add an extra equation which will force a downgraded error on the forecast:
frml(identity) checkf vixcheck = %if(vix<0.0,%na,vix)
group fullmodel %modeleqn(ivar,1) %modeleqn(ivar,2) %modeleqn(ivar,3) $
%modeleqn(ivar,4) %modeleqn(ivar,5) %modeleqn(ivar,6) intereq checkf
(Solving to %NA is a more common problem in forecasting than explosive values---the latter is generally associated with a non-linear model that has no solution, while the former is typically due to failing to initialize an exogenous variable).
Re: GIRFs for I-VAR
Posted: Fri Nov 17, 2017 3:38 am
by Jules89
Thank you,
as far as I understand that correctly an identity is added to the model that sets the value of VIXCHECK equal to %na at any point in time when the forecast for the VIX produces a negative value. I assume that you chose that condition, because around that point the IRFs are starting to become explosive. Having an identity in the model that produces %NAs produces an error in the forecast instruction, which is somewhat downgraded the error of explosive draws is no longer shown. Because the downgraded error does not break the loop (for what ever reason...it seems like there is some kind of ranking of errors, where some cause the program to stop while others do not), the nextiteration starts.
However that does not seem to work... the error which breaks the loop still occurs for example for SEED 1234, 1, 2 and many more. The code is posted below...
Code: Select all
frml(identity) checkf vixcheck = %if(vix<0.0,%na,vix)
group fullmodel %modeleqn(ivar,1) %modeleqn(ivar,2) %modeleqn(ivar,3) $
%modeleqn(ivar,4) %modeleqn(ivar,5) %modeleqn(ivar,6) intereq
do i = 1,nvar
set girf(i) hstart hend = 0.0
set girfzlb(i) hstart hend = 0.0
end do i
seed 1234
compute nhist = 0
env traperrors
infobox(action=define,progress,lower=zlbstart,upper=zlbend) "Boot zlb"
do hist = zlbstart, zlbend
infobox(current=hist)
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
compute countrep = 0 ;* to count for repetitions in the draws
do draw=1,ndraws
*************************************************************
*** Draw random entries for residuals for forecast period ***
*************************************************************
:newtry
if countrep == ndraws{
compute nhist=nhist-1
break
}
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)
*********************************
*** Check for explosive IRFs ***
*********************************
compute nexpl = 0
compute test1 = wshock(1)(hend)-base(1)(hend)
compute test2 = wshock(1)(hstart)-base(1)(hstart)
if test1>=chol(1,1)*10*abs(delta).OR.test1<=(-1)*chol(1,1)*10*abs(delta).OR.test2>=chol(1,1)*5*abs(delta)
compute nexpl = nexpl+1
if nexpl>0{
compute countrep = countrep + 1
goto newtry
}
******************************************************
*** Accumulate GIRFs to calculate average later on ***
******************************************************
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
infobox(action=remove)
do i = 1,nvar
set girfzlb(i) hstart hend = girfzlb(i)/nhist
end do i
I also tried the check function for FFR, as you mentioned earlier that the FFR going negative has something to do with the explosive behavior. But that doesn't work either. Is there maybe some alternative to the forecast instruction?
Or do I implemented your above suggestion wrong?
However when I reduce the umber of steps in the GIRF it seems to work fine. So the explosive behavio starts somwhere between forecast step 20 and 25. It would still be nice to fix this problem such that I can use any number of forecasting steps
Thanks
Best
Jules
Re: GIRFs for I-VAR
Posted: Wed Nov 22, 2017 9:49 am
by Jules89
I was again thinking about the above stated problem with the explosive draws within the forecast instruction.
The FORECAST instruction does not provide any way to discard Forecasts which are getting too large. However at some point the FORECAST instruction produces this error and the loop breaks.
Therefore I cannot apply the criterium the authors use to discard explosive draws, because I have to check for that after the FORECAST instruction, when the loop already stoped.
A possible solution to that problem would be to split the FORECAST into several parts. For example into 5 steps parts. Then for each part I can check whether the forecast gets too large and is eventually becoming explosive within the next iterations.
If it turns out that the forecast generates large values, I can apply the criterium of the authors to that and discard the forecast and jump to the next draw of the bootstrap. If the values are not too large I can go on with the next 5 interations.
Would you consider that as a possibility? If so how would I split the Forecast? It would be much easier is the forecast instruction would set the values to infinity (so a values thats larger than everything) instead of producing an error, maybe that is something for the suggestion box for later releases of RATS?
Best Jules
Re: GIRFs for I-VAR
Posted: Wed Nov 22, 2017 10:30 am
by TomDoan
The current trial version of RATS downgrades the explosive results from an error to a warning, and let's you test the %CONVERGED variable to see if the forecast calculation worked.
While this would meet your requirements, you need to know that the problem isn't so much with the software as with this model. Even in the non-ZBT initial conditions, the simulated values for both FFR and VIX often go negative, just not so badly that the whole model blows up. I'm not sure what the point is of "impulse response functions" which are averaging in nonsensical simulations. (Basically, you have nonsense-nonsense=IRF).
Re: GIRFs for I-VAR
Posted: Wed Nov 22, 2017 10:35 am
by Jules89
The point is that I want to discard all draws which do not converge or produce expIosive draws. If for a set of draws this happens too often, then I discard this inital history and go on with the next one. Therefore all histories which produce nonsense IRFs will not be included. I am runnung on RATS 9.20d. So there is no possibility to do the mentioned stuff with the %CONVERGED? Or is this version late enough? If not can I update the Forecast instruction without additional costs?
Re: GIRFs for I-VAR
Posted: Wed Nov 22, 2017 10:56 am
by TomDoan
The current (i.e. built over the weekend) trial has that feature. We've never seen a model with such behavior---usually explosive means the model itself is badly flawed by construction, not that it goes ballistic due to random shocks.
I'm sorry, there are plenty of nonsensical solutions of that model that don't trigger the overflow. Because the IRF calculation subtracts one value from another, the fact that both are absurd is hidden.
Re: GIRFs for I-VAR
Posted: Wed Nov 22, 2017 11:11 am
by Jules89
So there will be a release of a new RATS version( i.e. 9.20.e or so?), which has this new feature on the forecast instruction. When will it be released? I have a single user license 9.20d (std) how can I upgrate that?
Thanks
Jules
Re: GIRFs for I-VAR
Posted: Wed Nov 22, 2017 11:19 am
by TomDoan
You can download the trial which has it. We don't know when we will have the next build of the standard software.