May the reason for this be that the variable "money" and "D1p" aren't stationary? Because in a TVAR the variables have to be stationary right? Otherwise shocks might drive you persistently in a regime, where the process stays extremely lone or evene forever, right?
Anyway, below is a code for the bootstrapped confidence bands. I am not quite sure whether this is correct, as the results seem odd to me. I use the median of all possible GIRFs which are estimated from the data simulated by the model, as the GIRFs central tendency. The confidence bands are derived by using upper and lower quantiles. Could you have a look on the code?
Code: Select all
*******************
**** Load Data ****
*******************
cal(q) 1947
allocate 1997:4
open data flofund2.dat
data(unit=data,org=obs) 1947:1 1997:4 $
date nfbfla nfbstk nfbfln nfcfla nfcstk nfcfln cmpapfln cmpapfla cmpapstk rgdpcw pgdpcw fyff m2 cp46 t6
****************************
**** Control parameters ****
****************************
compute nvar =4 ;* #-Dependent variables
compute maxlag=4 ;* Lag-length
comp upper =0 ;* Regime Selector. Change 0 for other regime
comp horizon =12 ;* IRF horizon
comp nboot =100 ;* #-of IRF Bootstrap Replications
compute shock= +1.0 ;* size and signs of shocks in terms of standard deviations
compute sstart = 1959:1 ;* Estimation start
compute send = 1997:3 ;* Estimation end
compute d = 1 ;* Threshold Lag
compute ciboot = 500 ;* #-of Bootstrap replications for CIs
compute thresh = 0.48667 ;* Threshold value estimated before
compute pi=.15 ;* Fraction of threshold values excluded at each end
******************************
**** Data Transformations ****
******************************
set mix = (nfbstk+nfcstk)/(cmpapstk+nfbstk+nfcstk)
set gdp = rgdpcw
set pgdp = pgdpcw
set d1mix = (mix-mix{1})
set cpbill1 = cp46-t6
set d1gdp = 400*log(gdp/gdp{1})
set d1pgdp = 400*log(pgdp/pgdp{1})
set ly = log(gdp)
set d1y = d1gdp
set d1p = d1pgdp
set money = fyff
set credit = cpbill1
***********************
**** Storage Space ****
***********************
dec vect[frml] tvarf(nvar) ;* Vector of formulas of T-VAR equations, which will include the switching
dec vect[series] bootu(nvar) data(nvar)
dec series[vect] bootuv bootres
dec vect[string] shortlabels(nvar) longlabels(nvar)
*****************************************
**** Set up the series for the T-VAR ****
*****************************************
dec vect[int] depvars(nvar) ;* Vector of Indicators of Dependent Variables
compute depvars =||d1y,d1p,money,credit||
***********************************************
**** set thresholds and Threhold Formulas ****
***********************************************
dec vector macoeffs
compute malength = 2
filter(type=lagging,span=malength) credit / credthr ;* Moving Average of the Threshold Variable
compute macoeffs=%fill(malength,1,1.0/malength)
equation(identity,coeffs=macoeffs) threqn credthr ;* This generates an equation (and from it a FRML) to compute the moving average
# credit{0 to malength-1}
frml(equation=threqn,identity) thrfrml
set dup = thrfrml{d}>thresh ;* Dummy for the upper regime
set dlo = 1.-dup ;* Dummy for the lower regime
***********************************************
**** Estimation of the Model coefficients *****
***********************************************
system(model=basevar)
variables depvars
lags 1 to maxlag
det constant
end(system)
compute rstart=sstart+maxlag
compute rend=send
dec vect[series] resup(nvar) reslo(nvar) ;* Vector of Series of Upper and Lower regime residuals splitted up
dec vect[frml] fitup(nvar) fitlo(nvar) ;* Vector of formulas of fitted VAR equations (includes the estimated coefficients)
*** Upper regime estimation
estimate(noprint,smpl=dup,resids=resup,title="Upper Regime") rstart rend ;* Estimate under the upper regime
do i=1,nvar
frml(equation=%modeleqn(basevar,i)) fitup(i) ;* Define regime-specific FRML’s
end do i
compute vup=%sigma ;* Upper Regime Covariance Matrix
*** Lower regime estimation
estimate(noprint,smpl=dlo,resids=reslo,title="Lower Regime") rstart rend ;* Estimate under the lower regime
do i=1,nvar
frml(equation=%modeleqn(basevar,i)) fitlo(i) ;* Define regime-specific FRML’s
end do i
compute vlo=%sigma ;* Lower Regime Covariance Matrix
********************************************************************************
**** Compute the upper and lower branch covariance matrices, their Choleski ****
**** factors and the inverse factor. (The inverse is for standardizing the ****
**** residuals and the factor for mapping standardized residuals back to ****
**** the true levels.) ****
********************************************************************************
compute sup =%decomp(vup) ;* Upper Regime Cholesky Decomp of Covariance Matrix
compute siup=inv(sup) ;* Upper Regime inverse Cholesky Decomp of Covariance Matrix
compute slo =%decomp(vlo) ;* Lower Regime Cholesky Decomp of Covariance Matrix
compute silo=inv(slo) ;* Lower Regime inverse Cholesky Decomp of Covariance Matrix
**************************************************
**** Compute (jointly) standardized residuals ****
**************************************************
dec series[vect] stdu
gset stdu rstart rend = %if(dup,siup*%xt(resup,t),silo*%xt(reslo,t)) ;* siup/silo are vect[series] %xt(resup,1) takes for t=1 the vector of realisations of resup
********************************
**** Save the original data ****
********************************
dec vect[series] data(nvar)
do i=1,nvar
set data(i) = depvars(i){0}
end do i
********************************************************************************
**** Define the switching formula. Because the bootstrap requires taking ****
**** the standardized residuals and reflating them using regime-specific ****
**** covariance matrices, the reflation part has to be incorporated into ****
**** the formula. Thus TVARF(i) is (in effect) an identity given the (still ****
**** to be created) bootres series. ****
********************************************************************************
do i=1,nvar
frml tvarf(i) depvars(i) = %if(thrfrml{d}>thresh,$
fitup(&i)+%dot(%xrow(sup,&i),bootres),$
fitlo(&i)+%dot(%xrow(slo,&i),bootres))
end do i
*****************************************************************************
**** Build the threshold var model using the switching equations and the ****
**** definitional identity for the threshold variable. ****
*****************************************************************************
group tvar tvarf(1) tvarf(2) tvarf(3) tvarf(4) thrfrml
*************************************************************************
***********************************************************************************
**** Figure out which set of entries we want. This is the simplest way ****
**** to get this. dup{0} takes two values (0 or 1) but the order isn't known ****
**** in advance since it will depend upon which value is seen first in the ****
**** <<rstart>> to <<rend>> range. So we have to figure out which of values(1) ****
**** and values(2) (and the corresponding entries(1) and entries(2)) is the ****
**** one that we need, given the choice for <<upper>>. ****
***********************************************************************************
panel(group=dup{0},id=values,identries=entries) dup rstart rend
{
if upper==1.and.values(1)==1.or.upper==0.and.values(1)==0
compute rdates=entries(1)
else
compute rdates=entries(2)
}
******************************************************************************
**** Set up target series. There is one for each combination of test sign ****
**** and sizes on the shock, variable shocked and target variable. ****
******************************************************************************
compute wstart=rend+1,wend=wstart+horizon-1 ;* This is the working range for calculating the GIRF's
compute ninitial=%size(rdates) ;* #-of periods in the specified regime
dec rect[series] nlirfs(nvar,nvar)
do i=1,nvar
do j=1,nvar
set nlirfs(i,j) wstart wend = 0.0
end do j
end do i
*************************************
**** Storage Space for Bootstrap ****
*************************************
dec series[vect] stdub
dec vect[series] resupb(nvar) reslob(nvar) ;* Vector of Series of Upper and Lower regime residuals splitted up
dec vect[frml] fitupb(nvar) fitlob(nvar) ;* Vector of formulas of fitted VAR equations (includes the estimated coefficients)
dec vect[frml] tvarfb(nvar) ;* Vector of formulas of T-VAR equations, which will include the switching
dec vect[series] bootub(nvar)
dec series[vect] bootuvb bootresb
dec vect[series] datab(nvar)
dec rect[real] matgirf ;* For TRansformimng series into matrix
declare vect[series[vect]] girf(nvar*nvar)
do i = 1, nvar*nvar
gset girf(i) 1 ciboot = %zeros(horizon,1)
end do i
*************************************************
**** Additional Configurations for Bootstrap ****
*************************************************
dec vect[series] resample(nvar+1) ;* Set up storage space for the simulated Data
do i = 1,nvar
set resample(i) = %modeldepvars(tvar)(i){0}
end do i
set resample(5) = credthr
compute piskip=fix(pi*%nobs)+(maxlag*nvar+1) ;* To prevent overfitting, exclude the lower and upper pi% and number of regressors per equation from the estimation range
compute pistart=rstart+piskip
compute piend=rend-piskip
************************************************************************************************************************
**************************************** GIRFs Bootstrap within Bootstrap Start ****************************************
************************************************************************************************************************
infobox(action=define,lower=1,upper=ciboot,progress) "Residual Based Bootstrap"
do r = 1,ciboot
infobox(current=r)
*************************************
**** Simulate Data from the TVAR ****
*************************************
boot entr rstart rend
gset bootres rstart rend = stdu(entr(t))
forecast(model=tvar,from=rstart,to=rend,results=resample)
**********************************************
**** Estimate the Model with the new Data ****
**********************************************
filter(type=lagging,span=malength) resample(4) / credthrre ;* Moving Average of the Threshold Variable
compute macoeffs=%fill(malength,1,1.0/malength)
equation(identity,coeffs=macoeffs) threqnb credthrre ;* This generates an equation (and from it a FRML) to compute the moving average
# resample(4){0 to malength-1}
frml(equation=threqnb,identity) thrfrmlb
******************************************
****** Estimate the Threshold Value ******
******************************************
* Threshold Value is estimated by the value that maximizes the Log likelihood
* of the estimated model for a given threshold value
system(model=threshestim)
variables resample(1) resample(2) resample(3) resample(4)
lags 1 to maxlag
det constant
end(system)
estimate(noprint,model=threshestim) rstart rend
compute loglr=%logl
set copy rstart rend = credthrre{d} ;* Get a sorted copy of the threshold variable
order copy rstart rend
compute bestlogl=loglr
set lrstats pistart piend = 0.0
do pientry=pistart,piend
compute thresh=copy(pientry)
sweep(var=hetero,group=credthrre{d}<thresh) rstart rend ;* This allows for the covariance matrix to differ between regimes
# %modeldepvars(threshestim)
# %rlfromeqn(%modeleqn(threshestim,1))
if %logl>bestlogl ;* If this is best so far, save it
compute bestlogl=%logl,bestthresh=thresh
end do pientry
compute thresh = bestthresh
*********************************************************
**** given the new threhold value estimate the Model ****
*********************************************************
set dupb = thrfrmlb{d}>thresh ;* Dummy for the upper regime
set dlob = 1.-dupb
system(model=basevarb)
variables resample(1) resample(2) resample(3) resample(4)
lags 1 to maxlag
det constant
end(system)
*** Upper & Lower regime estimation
estimate(noprint,model=basevarb,smpl=dupb,resids=resupb) rstart rend ;* Estimate under the upper regime
estimate(noprint,model=basevarb,smpl=dlob,resids=reslob) rstart rend ;* Estimate under the lower regime
do i=1,nvar
frml(equation=%modeleqn(basevarb,i)) fitupb(i) ;* Define regime-specific FRML’s
end do i
compute vupb=%sigma ;* Upper Regime Covariance Matrix
*** Lower regime estimation
do i=1,nvar
frml(equation=%modeleqn(basevarb,i)) fitlob(i) ;* Define regime-specific FRML’s
end do i
compute vlob=%sigma ;* Lower Regime Covariance Matrix
compute supb =%decomp(vupb) ;* Upper Regime Cholesky Decomp of Covariance Matrix
compute siupb=inv(supb) ;* Upper Regime inverse Cholesky Decomp of Covariance Matrix
compute slob =%decomp(vlob) ;* Lower Regime Cholesky Decomp of Covariance Matrix
compute silob=inv(slob) ;* Lower Regime inverse Cholesky Decomp of Covariance Matrix
gset stdub rstart rend = %if(dup,siupb*%xt(resupb,t),silob*%xt(reslob,t)) ;* siup/silo are vect[series] %xt(resup,1) takes for t=1 the vector of realisations of resup
do i=1,nvar
frml tvarfb(i) resample(i) = %if(thrfrml{d}>thresh,$
fitupb(&i)+%dot(%xrow(supb,&i),bootresb),$
fitlob(&i)+%dot(%xrow(slob,&i),bootresb))
end do i
group tvarb tvarfb(1) tvarfb(2) tvarfb(3) tvarfb(4) thrfrmlb
panel(group=dupb{0},id=values,identries=entries) dupb rstart rend
{
if upper==1.and.values(1)==1.or.upper==0.and.values(1)==0
compute rdates=entries(1)
else
compute rdates=entries(2)
}
compute ninitial=%size(rdates)
do i=1,nvar
set datab(i) = resample(i){0}
end do i
************************
**** GIRF BOOTSTRAP ****
************************
do jrep=1,ninitial
********************************************************
**** Copy History into the dependent Variable Solts ****
********************************************************
compute basedate=rdates(jrep)
do i=1,nvar
move datab(i) basedate-maxlag basedate-1 resample(i) wstart-maxlag
end do i
******************************************
**** Loop over bootstrap replications ****
******************************************
do boot=1,nboot
****************************************
**** Generate the bootstrap shuffle ****
****************************************
boot rentries wstart wend rstart rend
***********************************************************************
**** Generate the base set of shocks by premultiplying the ****
**** bootstrapped standardized shocks by the factor of the overall ****
**** covariance matrix. ****
***********************************************************************
gset bootuvb wstart wend = %if(%ranflip(.5),+1,-1)*stdub(rentries(t))
gset bootresb wstart wend = bootuvb
forecast(model=tvarb,from=wstart,to=wend,results=base)
do jshock=1,nvar
**************************************************************
**** Patch over the component for which are computing the ****
**** response with the selected size and sign. ****
**************************************************************
compute bootresb(wstart)=bootuvb(wstart)
compute bootresb(wstart)(jshock)=shock
forecast(model=tvarb,from=wstart,to=wend,results=withshock)
do i=1,nvar
set nlirfs(i,jshock) wstart wend = nlirfs(i,jshock)+withshock(i)-base(i)
end do i
end do jshock
end do boot
end do jrep
do j=1,nvar
do i=1,nvar
set nlirfs(i,j) wstart wend = nlirfs(i,j)/(nboot*ninitial)
end do i
end do j
*********************
**** Bookkeeping ****
*********************
*****************************************************************************************************
** Converts the series of IRFs for different shocks into matrices. Gives you k (#-shocks) **
** matrices (i,j), where i is the ith response step of each variable. **
** j is equal to nvar(=4)*nvar(4), coulumns 1 to 4 are the responses of all variables to a shock **
** in variable 1, coulums 5 to 8 are the responses of all variables to a shock in variable 2, and **
** so on. **
*****************************************************************************************************
make matgirf wstart wend
# nlirfs
*****************************************************************************************************
** Plugs the columns of the above matrices in the storage space vectors for j = 1,...,nvar*nvar. **
** GIRF vectors of series of vectors: j=1,..4 are responses of veriable all variables **
** to a shock in variable 1 for draw r, j=5,...,8 are the responses of all variables **
** to a shock in variable 2 and so on **
*****************************************************************************************************
do j =1,nvar*nvar
do i =1,horizon
compute girf(j)(r)(i) = matgirf(i,j)
end do i
end do j
end do r
infobox(action=remove)
*********************************
**** Monte Carlo Integration ****
*********************************
declare vect[vect[vect]] Q(nvar*nvar)
compute lowerci = 0.16
compute middle = 0.50
compute upperci = 0.84
do i = 1,nvar*nvar
@mcmcpostproc(ndraws=ciboot,percentiles=||lowerci,middle,upperci||,quantiles=Q(i)) girf(i)
end do i
declare vect[series] ciup(nvar*nvar)
declare vect[series] center(nvar*nvar)
declare vect[series] cilow(nvar*nvar)
do i= 1, nvar*nvar
set cilow(i) 1 horizon = Q(i)(t)(1)
set center(i) 1 horizon = Q(i)(t)(2)
set ciup(i) 1 horizon = Q(i)(t)(3)
end do i
declare rect[series] cilomat(nvar,nvar)
declare rect[series] centermat(nvar,nvar)
declare rect[series] ciupmat(nvar,nvar)
compute k = 1
do i = 1,4
do j= 1,4
set cilomat(j,i) = cilow(k)
set centermat(j,i) = center(k)
set ciupmat(j,i) = ciup(k)
compute k=k+1
end do j
end do i
compute shortlabels=||"Output","Inflation","Money","Credit"|| ;* Labels Short
spgraph(hfields=4,vfields=4,header="Shocks",ylabels=shortlabels,xlabels=shortlabels)
do i = 1,nvar
do j = 1,nvar
graph(nodates) 3
# cilomat(j,i) / 2
# centermat(j,i) / 1
# ciupmat(j,i) / 2
end do j
end do i
spgraph(done)