*
* Replication file for Balke(2000), "Credit and Economic Activity:
* Credit Regimes and Nonlinear Propagation of Shocks," Review of
* Economics and Statistics, vol 82, 344-349.
*
* This program computes average IRF for output, conditional on being in
* the upper or lower regime, averaging across initial conditions and then,
* for each initial setting, across bootstrapped residuals.
*
* Control parameters
OPEN DATA "C:\Users\J. Jiang\Desktop\D\data\rawdata6.csv"
CALENDAR(M) 1968:4
DATA(FORMAT=PRN,ORG=COLUMNS,DATEFORM="m/d/y") 1968:04 2016:07 g_pm b_y cpi sp500 ffr brent_oil exch_cny m2 $
 ted vix baa WPU1594 cpigpctd3m


set logg  = log(g_pm)
set dlogg = logg- logg{1}
set logp  = log(cpi)
set pi    = logp - logp{1}
set nrg   = dlogg
set rrg   = dlogg - pi

set logo = log(brent_oil)
set nro = logo -logo{1}
set rro = nro - pi

set logjppi =log(WPU1594)
set dlnjppi = logjppi - logjppi{1}

comp nvar    =5
comp horizon =8


*
* Number of bootstrap replications used in computing GIRF's
*
comp nkrep   =100
*
* Change to upper=0 to get the WENDlower regime
*
comp upper   =1
*
comp [vector] shocksizes=||1.0,2.0||
comp [vector] shocksigns=||1.0,-1.0||
*


*
* These are the desired start and end of the sample to be used in the
* estimation. However, the transformations use the entire range.
*
compute sstart = 1987:7
compute send   = 2014:12
*
dec vector[series] res(nvar) vres(nvar)
dec vector v(nvar) resmat(nvar)
dec vect[int]      depvars(nvar)
dec rect[int]      laglengths(nvar,nvar)
dec vect[equation] eqn(nvar)
dec vect[series]   resup(nvar) resdn(nvar)
dec rect[frml]     fitud(nvar,2)
dec vect[frml]     tvarf(nvar) empty(nvar) tfrml(nvar) rfrml(nvar)
dec vect[series]   bootu(nvar) data(nvar)
dec series[vect]   bootuv bootres
*
dec vect[string]   shortlabels(nvar) longlabels(nvar)
declare vector[labels] lab(nvar)
*

compute depvars    =||rro,dlnjppi,pi,ted,rrg||
compute shortlabels=||"Real Return of Oil","Jewlry PPI Growth","Inflation","Ted Spread","Real Return of Gold"||
compute longlabels =||"Real Return of Oil","Jewlry PPI Growth","Inflation","Ted Spread","Real Return of Gold"||
compute lab        =|| "O", "J", "I", "Ted", "G" ||


*
*  set thresholds, mas, and delays
*
compute d = 0
*
* 15% window
*
dec vector macoeffs
compute malength = 1  ; comp thresh = 0.49 ;* cpbill, ff
*compute malength = 6 ; comp thresh = -0.00275 ;* dmix, ff
*compute malength = 6 ; comp thresh = -3.19 ;* dsd, ff
*
* This generates the actual moving average of the data
*
filter(type=lagging,span=malength) ted / tedthr
*
* This generates an equation (and from it a FRML) to compute the moving
* average of the data
*
dec vector macoeffs
dim macoeffs(malength)
comp macoeffs = %const(1./malength)
equation(identity,coeffs=macoeffs) threqn tedthr
# ted{0 to malength-1}
frml(equation=threqn,identity) thrfrml
*
* Lag lengths are allowed to differ among equations. These give the lag
* lengths with equation in a row and variable in a column.
*
input laglengths
  3 3 3 3 3
  3 3 3 3 3
  3 3 3 3 3
  3 3 3 3 3
  3 3 3 3 3
*************************************************************************
*
* Dummy variables for the two regimes
*
set d1 = thrfrml{d}>thresh
set d2 = 1.-d1

print / d1
*
compute maxlag=d+(malength-1)
inquire(reglist) rstart<<sstart rend<<send
# tedthr{d}
*
* Create an equation for each variable with the number of lags for each
* endogenous variable given by <<laglengths>>.
*
do i=1,nvar
   compute [vect[integer]] rl=||constant||
   do j=1,nvar
      compute rl=%rladdlaglist(rl,depvars(j),%seq(1,laglengths(i,j)))
      compute maxlag=%imax(maxlag,laglengths(i,j))
   end do j
   equation eqn(i) depvars(i)
   # rl
end do i
*
compute rstart=sstart+maxlag,rend=send
*display %datelabel(wstart) %datelabel(wend)
*display %datelabel(rstart) %datelabel(rend)
*display %datelabel(sstart) %datelabel(send)

do i=1,nvar
   linreg(equation=eqn(i)) * rstart rend
   frml empty(i) bootu(i) = 0.0
end do i
*
do i=1,nvar
   disp
   disp
   disp shortlabels(i)+" regression in upper regime"
   linreg(smpl=d1,equation=eqn(i),frml=fitud(i,1)) * rstart rend resup(i)
   disp
   disp
   disp shortlabels(i)+" regression in lower regime"
   linreg(smpl=d2,equation=eqn(i),frml=fitud(i,2)) * rstart rend resdn(i)
   set res(i) rstart rend = %if(d1,resup(i),resdn(i))
end do i
*
* 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.)
*
vcv(matrix=vup)
# resup
compute sup =%decomp(vup)
compute siup=inv(sup)

vcv(matrix=vdn)
# resdn
compute sdn =%decomp(vdn)
compute sidn=inv(sdn)
*
* Compute the joint covariance matrix and its factor.
*
vcv(matrix=vsigma)
# res
compute s=%decomp(vsigma)
*
* Compute (jointly) standardized residuals. Since we only use these as a
* VECTOR across variables at T (never as individual series), we define
* it as a SERIES[VECTOR].
*
dec series[vect] stdu
gset stdu rstart rend = %if(d1,siup*%xt(resup,t),sidn*%xt(resdn,t))
*
* Save the original data since we'll overwrite it as part of the
* bootstrap.
*
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,$
        fitud(&i,1)+%dot(%xrow(sup,&i),bootres),$
        fitud(&i,2)+%dot(%xrow(sdn,&i),bootres))
end do i
*
* Build the threshold var model using the switching equations and the
* definitional identity for the threshold variable.
*
*************************************************************************
*
* This is application specific. A different number of variables requires
* a different list of the tvarf's.
*
group tvar tvarf(1) tvarf(2) tvarf(3) tvarf(4) tvarf(5) thrfrml
*
*************************************************************************
*
* Figure out which set of entries we want. This is the simplest way
* to get this. d1{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=d1{0},id=values,identries=entries) d1 rstart rend
{
if upper==1.and.values(1)==1.or.upper==0.and.values(1)==0
   compute rdates=entries(1),nrep=%size(rdates)
else
   compute rdates=entries(2),nrep=%size(rdates)
}
*
* Set up target series. There is one for each combination of test sign
* and sizes on the shock, variable shocked and target variable.
*
declare real size sign
compute nexp=%size(shocksizes)*%size(shocksigns)
dec vect[rect[series]] girfs(nexp)
do i=1,nexp
   dim girfs(i)(nvar,nvar)
   clear(zeros) girfs(i)
end do i
*
* This is the working range for calculating the GIRF's
*
compute wstart=rstart,wend=rstart+horizon-1
*
* Outer loop is over the initial conditions, which walk through the data
* points in the desired regime. The residuals are bootstrapped by taking
* the standardized residuals (across the entire range), shuffling them,
* and reflating them based upon the current threshold value in the
* generated series.
*
infobox(action=define,lower=1,upper=nrep,progress) "Bootstrapping Across Initial Values"
do jrep=1,nrep
   infobox(current=jrep)
   *
   * Copy observed data into depvar slots
   *
   compute basedate=rdates(jrep)
   do i=1,nvar
      set depvars(i) wstart-maxlag wstart-1 = data(i)(t-wstart+basedate)
   end do i
   *
   * Loop over bootstrap replications
   *
   do krep=1,nkrep
      *
      * 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 bootuv wstart wend  = %if(%ranflip(.5),+1,-1)*stdu(rentries(t))
      *
      gset bootres wstart wend = bootuv
      *
      forecast(model=tvar,from=wstart,to=wend,results=base)
      compute ifill=0
      dofor sign = 1 -1
         dofor size = 1 2
            compute ifill=ifill+1
            do jshock=1,nvar
               *
               * Patch over the component for which are computing the
               * response with the selected size and sign.
               *
               compute bootres(wstart)=bootuv(wstart)
               compute bootres(wstart)(jshock)=sign*size
               *
               forecast(model=tvar,from=wstart,to=wend,results=withshock)
               do i=1,nvar
                  set girfs(ifill)(i,jshock) wstart wend = girfs(ifill)(i,jshock)+withshock(i)-base(i)
               end do i
            end do jshock
         end do sign
      end do size
   end do krep
end do jrep
infobox(action=remove)
*
do k=1,nexp
   do i=1,nvar
      do j=1,nvar
         set girfs(k)(i,j) wstart wend = girfs(k)(i,j)/(nkrep*nrep)
      end do j
   end do i
end do k
*
* Move the data back
*
do i=1,nvar
   set depvars(i) = data(i)
end do i
*
*************************************************************************
*
dec vect[series] graphs(4)
compute klabels=||"+1 SD","+2 SD","-1 SD","-2 SD"||
*
spgraph(vfields=nvar,hfields=1,footer=$
   "Figure 2.1 Response of Stock Return to Shocks, Conditional on Regime")
do j=1,nvar
   do i=1,4
      set graphs(i) 1 horizon = girfs(i)(1,j)(t+wstart-1)
   end do i
   graph(series=graphs,number=0,key=upright,klabels=klabels,$
     header="High Regime: Response of Stock",subheader="Shock to "+shortlabels(j))
end do j
spgraph(done)

*************************************************************************

spgraph(vfields=nvar,hfields=1,footer=$
   "Figure 2.2 Response of Bond Return to Shocks, Conditional on Regime")
do j=1,nvar
   do i=1,4
      set graphs(i) 1 horizon = girfs(i)(2,j)(t+wstart-1)
   end do i
   graph(series=graphs,number=0,key=upright,klabels=klabels,$
     header="High Regime: Response of Bond",subheader="Shock to "+shortlabels(j))
end do j
spgraph(done)

**************************************************************************

spgraph(vfields=nvar,hfields=1,footer=$
   "Figure 2.3 Response of Gold Return to Shocks, Conditional on Regime")
do j=1,nvar
   do i=1,4
      set graphs(i) 1 horizon = girfs(i)(3,j)(t+wstart-1)
   end do i
   graph(series=graphs,number=0,key=upright,klabels=klabels,$
     header="High Regime: Response of Gold",subheader="Shock to "+shortlabels(j))
end do j
spgraph(done)

**************************************************************************
spgraph(vfields=nvar,hfields=1,footer=$
   "Figure 2.4 Response of Inflation to Shocks, Conditional on Regime")
do j=1,nvar
   do i=1,4
      set graphs(i) 1 horizon = girfs(i)(4,j)(t+wstart-1)
   end do i
   graph(series=graphs,number=0,key=upright,klabels=klabels,$
     header="High Regime: Response of Inflation",subheader="Shock to "+shortlabels(j))
end do j
spgraph(done)