Error Bands on IRFs in TVAR

[Jules89]

Dear Tom,

I have two questions on Threshold VAR models. Since I am an absolute beginner in this topic this might already have been dealt with:

1) TVARs switch coefficients when an estimated threshold value has been exceeded. Assume I want to analyse shocks under a zero lower bound regime on interest rats, then technically I don't need to estimate the threshold, right? I would just have to set the threshold value of the policy rate to zero.

2) To analyse the IRFs generalized IRFs have to be computed. In the Balke (2000) paper there have no error bands been computed. Is it possible to compute error bands for GIRFs for example by bootstrapping or monte carlo integration? A reference would be Galvao and Marcellino (2014) -"The effects of the monetary policy stance on the transmission mechanism", they use a residual-based bootstrap method to the estimated nonlinear model.

3) Is it possible to calculate regime depentent IRFs including error bands?



It would be very nice if you could suggest some example code related to the above posted questions.

Thank you very much in advance

Best

Jules

[TomDoan]

I have two questions on Threshold VAR models. Since I am an absolute beginner in this topic this might already have been dealt with:

1) TVARs switch coefficients when an estimated threshold value has been exceeded. Assume I want to analyse shocks under a zero lower bound regime on interest rats, then technically I don't need to estimate the threshold, right? I would just have to set the threshold value of the policy rate to zero.

If the threshold is known, then you don't have to estimate it. However, you have to have data points on both sides of the threshold to estimate the model, and, in fact, have to have enough data points in each regime to estimate a VAR. This isn't a theoretical model, it's empirical, so if something hasn't happened, or only rarely has happened, you can't use it.

2) To analyse the IRFs generalized IRFs have to be computed. In the Balke (2000) paper there have no error bands been computed. Is it possible to compute error bands for GIRFs for example by bootstrapping or monte carlo integration? A reference would be Galvao and Marcellino (2014) -"The effects of the monetary policy stance on the transmission mechanism", they use a residual-based bootstrap method to the estimated nonlinear model.

It's a bit odd that that paper doesn't reference Balke, given that other than the adding bootstrap over parameters, their model is exactly the same as Balke's. You can bootstrap over the model, as they do in that paper---it's just an added complication.

3) Is it possible to calculate regime depentent IRFs including error bands?

How is that different from question 2?

[Jules89]

Thank you for the quick reply. Is there any code example for a similar bootstrap ?

Best jules

[TomDoan]

No. But it wouldn't be a big issue. From your description, it sounds like you have a bigger problem with getting what you're doing into a TVAR framework.

[Jules89]

Dear Tom,

comming back to my question regarding the residual based bootstrap for the TVAR, which I would like to use to simulate data from the estimated model to generate confidence bands for the GIRFs.

I used the Balke(2000) code and tried to simlate data from it by using resampled residuals. The problem is that the simulation generates explosive time series. The code is given by:

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 =1                                                              ;* 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 = 20                                                        ;* #-of Bootstrap replications for CIs
compute thresh = 0.48667                                                   ;* Threshold value estimated before


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


****************************************************************************
**** Create a FRML to determine which regime holds in a simulation. (It ****
**** will be 1 for upper and 0 for lower)                               ****
****************************************************************************

*frml(identity) upperfrml dup = thrfrml{d}>thresh


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



***********************
**** Simulate Data ****
***********************


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


do rep = 1, 10

   boot entr rstart rend

   gset bootres rstart rend = stdu(entr(t))
   forecast(model=tvar,from=rstart,to=rend,results=resample)

   spgraph(hfields=4,vfields=1)
      do i = 1,4
         graph 2
         # data(i) rstart rend
         # resample(i) rstart rend
      end do i
   spgraph(done)

end do rep

Is this explosive behavor a property of the model, or did I do something wrong?

Thank you in advance

Best Jules

[TomDoan]

The dominant root of the lower partition is about 1.2. If (as seems to be the case), the lower regime doesn't necessarily drive the credit variable up towards the threshold, then eventually the process will diverge. That's not as obvious with the IRF's which only do 15 steps and which take the difference between the results with and without one single change to the shocks (thus big number-big number can be small number).

One thing to note is that computing error bands in a model which requires simulations just to compute the IRF's in the first place requires a staggering amount of calculation.

[Jules89]

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?
The Dickey-Fuller test results are:

Code: Select all

Dickey-Fuller Unit Root Test, Series D1P
Regression Run From 1959:04 to 1997:03
Observations        153
With intercept
With 2 lags chosen from 4 by GTOS/t-tests(0.050)

Sig Level    Crit Value
1%(**)         -3.47396
5%(*)          -2.88035
10%            -2.57669

T-Statistic    -1.93483


Dickey-Fuller Unit Root Test, Series MONEY
Regression Run From 1960:01 to 1997:03
Observations        152
With intercept
With 3 lags chosen from 4 by GTOS/t-tests(0.050)

Sig Level    Crit Value
1%(**)         -3.47423
5%(*)          -2.88047
10%            -2.57675

T-Statistic    -2.42320

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?

Best

Jules

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)

[TomDoan]

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?

That's not true. In the example below, this generates a non-explosive process even though one branch is non-stationary and one is explosive. For TAR models, the behavior of the combined process depends upon all sorts of things---for instance, the variance of the shocks matters (a big variance here does create a problem because it overwhelms the -1.0 in the upper branch that's driving the process back towards the threshold).

Code: Select all

dec series y
frml yfrml y = %if(y{1}<=1.0,y{1},-1.0+1.1*y{1})
set y 1 5000 = 0.0
simulate(model=yfrml,from=2,to=5000,cv=||1.0||,results=test)

[TomDoan]

Regarding your program, inside the loop, you're redefining the THRESH value that's needed to do the outer bootstrap.

Given the behavior of the bootstrapped data, I'm not surprised that you're getting odd results. Once those get off track, they very quickly turn into y(t)=1.2y(t-1) for each of the series (basically, you see only the dominant root and nothing else)

[Jules89]

Thanks Tom,

I have three questions:

1)when I get this right, then stationarity of the variables seems not to be a major issue in TVAR models. A summary paper by Hubrich and Teräsvirta (2013) with the title "Thresholds and Smooth Transitions in Vector Autoregressive Models (link: https://pure.au.dk/ws/files/54473123/rp13_18.pdf), defines the TVAR for stationary time series (see page 7). Basicly all TVAR papers I read use stationary time series (or variables that should be stationary). Would you recommend to transform variables to stationarity if they arn't or would you go with the standard behavior in linear VARs and leav everything in levels?

2) Right I have to assign a different name for it within the reestimation (maybe call it threshb or so, such that in the outer bootstrap I have the "old" threshvariable). I wanted to estimate a new threshold for the simulated data. Should I stick to the already estimated threshold and do not reestimate it for every "draw" of new data?

3) So the bad news is that for this data set bootstrapping confidence bands is not a good idea, right?

Thank you very much

Best

Jules

[TomDoan]

Thanks Tom,

I have three questions:

1)when I get this right, then stationarity of the variables seems not to be a major issue in TVAR models. A summary paper by Hubrich and Teräsvirta (2013) with the title "Thresholds and Smooth Transitions in Vector Autoregressive Models (link: https://pure.au.dk/ws/files/54473123/rp13_18.pdf), defines the TVAR for stationary time series (see page 7). Basicly all TVAR papers I read use stationary time series (or variables that should be stationary). Would you recommend to transform variables to stationarity if they arn't or would you go with the standard behavior in linear VARs and leav everything in levels?

I'm not really sure what you're saying. The 2013 paper keeps working with "stationary" models because you actually can do proofs of various conditions for them. They rather explicitly allow for non-stationary models as well, they just don't do much with them.

The model is what the model is. If the threshold is seen as depending upon the non-stationary credit variable, then you can't replace that with a stationary alternative.

2) Right I have to assign a different name for it within the reestimation (maybe call it threshb or so, such that in the outer bootstrap I have the "old" threshvariable). I wanted to estimate a new threshold for the simulated data. Should I stick to the already estimated threshold and do not reestimate it for every "draw" of new data?

You're doing a parametric bootstrap treating the original estimate as fixed, but part of the uncertainty in the model is that that's not known so you should be re-estimating it.

3) So the bad news is that for this data set bootstrapping confidence bands is not a good idea, right?

Correct. You probably would have an easier time doing MCMC (which treats the data as given)---the only thing non-linear is the threshold value; given it, everything else can be drawn by standard VAR techniques.

[Jules89]

Thank you very much,

the last point sounds very interesting.
How would I estimate that with MCMC? As far as I understand you mean that I estimate the threshold in the usual way (as in Balke 2000).
Since the data are the same and I am not generating any new data the threshold keeps being the same. From the samplesplit implied by th eestimated threshold I get two data sets, for which I set up two likelihoods and two priors. Then I draw for the two resulting VARs the coefficients and run the GIRF Bootstrap.

That would mean I have model uncertainty with respect to the VAR parameters but not with respect to the threshold, since this is given by the data, correct?

how would I get the sample split data to set up the two gibbs samplers?


Thank you very much


best Jules

[TomDoan]

Given the threshold value, you just have two separate VAR's. Sampling from those will take care of the covariance matrices and VAR coefficients. Those are "one-off" calculations since the sample changes, so it's not the super-efficient calculation used by @MONTEVAR which can do everything for each draw using the least squares statistics, but it's no more complicated than having to re-estimate with bootstrapped data. The likelihood in the threshold value is discontinuous (since it only changes at the observed data points of the threshold variable), but you just need to come up with a proposal density for it, which could be continuous, to do Independence Chain M-H.

[Jules89]

Dear Tom,

thank you for your help, but I may have missed an argument. Within this bayesian setup I don't understand what you mean by "Given the threshold value, you just have two separate VAR's. Sampling from those will take care of the covariance matrices and VAR coefficients."?
I have two VARs and I draw the coefficients of those from the respective posteriors. Given the threshold value, I would draw for each of the two VARs the parameters similar to the @MONTEVAR procedure. Then I set up the TVAR as before and run the GIRF bootstrap.

What do you mean by " Those are "one-off" calculations since the sample changes"? Given the two likelihoods for the two VARs and the respective posteriors, where does the sample change and why?

I assume that you mean I should sample the threshold value as well, since you said "The likelihood in the threshold value is discontinuous (since it only changes at the observed data points of the threshold variable), but you just need to come up with a proposal density for it, which could be continuous, to do Independence Chain M-H."
But how would such a proposal density look like? I have never seen a paper doing something like that for TVARs and it sounds very challenging... Could you give me a reference or tell me a bit more detailed how the sampling steps would look like?

Thank you in advance

Best

Jules

[TomDoan]

Dear Tom,

thank you for your help, but I may have missed an argument. Within this bayesian setup I don't understand what you mean by "Given the threshold value, you just have two separate VAR's. Sampling from those will take care of the covariance matrices and VAR coefficients."?
I have two VARs and I draw the coefficients of those from the respective posteriors. Given the threshold value, I would draw for each of the two VARs the parameters similar to the @MONTEVAR procedure. Then I set up the TVAR as before and run the GIRF bootstrap.

What do you mean by " Those are "one-off" calculations since the sample changes"? Given the two likelihoods for the two VARs and the respective posteriors, where does the sample change and why?

One VAR is for the lower regime, one for the upper. Change the threshold, change the sample covered by each. Draws for a VAR with the standard non-informative prior only require the OLS coefficients, covariance matrix and X'X matrix for a particular sample which is why you can compute those sufficient statistics once for a fixed sample VAR and generate draws very quickly. In this case, the samples (likely) change with each sweep, so you can't use the old information.

I assume that you mean I should sample the threshold value as well, since you said "The likelihood in the threshold value is discontinuous (since it only changes at the observed data points of the threshold variable), but you just need to come up with a proposal density for it, which could be continuous, to do Independence Chain M-H."
But how would such a proposal density look like? I have never seen a paper doing something like that for TVARs and it sounds very challenging... Could you give me a reference or tell me a bit more detailed how the sampling steps would look like?

It would be Gibbs sampling (or Metropolis within Gibbs). You do everything but the threshold given the threshold (again, that's just two VAR's), and then the threshold given everything else, which is not standard, but the conditional likelihood depends only upon the residuals from the two VAR's and the two covariance matrices. As to what might be an appropriate proposal for the threshold, that would depend entirely upon the application. In some cases, it might be rather flat; in others it might need to be much more peaked. However, it's just a single parameter, so the sampler's likely to be fairly robust.

As I'm sure you're aware, there's very little done with error bands on TVAR's. It just seems that this is a case where parametric bootstrapping is fixing a bit too much at sample values to be workable.