*
* 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
*
comp nvar    =3
comp horizon =12
*
* Number of bootstrap replications used in computing GIRF's
*
comp nkrep   =500
*
* Change to upper=0 to get the lower regime
*
comp upper   =2
*
comp [vector] shocksizes=||1.5,3.0||
comp [vector] shocksigns=||1.0||
*
OPEN DATA "C:\Users\Juan Manuel\Desktop\STAR\pass_through12.rat"
CALENDAR(M) 1945:1
DATA(FORMAT=RATS) 1945:01 2014:07  TCE S IPC2 IPC2_D11 M0 M1_TODO M2_TODO M3_TODO



SET VE = LOG(TCE)-LOG(TCE{1})
SET INFLAIPC = LOG(IPC2)-LOG(IPC2{1})
SET INFLAIPCD = LOG(IPC2_D11)-LOG(IPC2_D11{1})
SET VSAL = LOG(S)-LOG(S{1})
SET BM = LOG(M0)-LOG(M0{1})
SET M1 = LOG(M1_TODO)-LOG(M1_TODO{1})
SET M2 = LOG(M2_TODO)-LOG(M2_TODO{1})
SET M3 = LOG(M3_TODO)-LOG(M3_TODO{1})

dummy(ao=1959:01) dummy5901
dummy(ao=1973:06) dummy7306
dummy(ao=1974:04) dummy7404
dummy(ao=1975:06) dummy7506
dummy(ao=1976:03) dummy7603
dummy(ao=1982:07) dummy8207
dummy(ao=1985:07) dummy8507
dummy(ao=1987:10) dummy8710
dummy(ao=1989:04) dummy8904
dummy(ao=1989:05) dummy8905
dummy(ao=1989:06) dummy8906
dummy(ao=1989:07) dummy8907
dummy(ao=1989:08) dummy8908
dummy(ao=1989:09) dummy8909
dummy(ao=1989:12) dummy8912
dummy(ao=1990:01) dummy9001
dummy(ao=1990:03) dummy9003
dummy(ao=1990:04) dummy9004
dummy(ao=1991:02) dummy9102
dummy(ao=1992:10) dummy9210
dummy(ao=2002:01) dummy0201
dummy(ao=2002:04) dummy0204
dummy(ao=2011:04) dummy1104
dummy(ao=2014:01) dummy1401

set x = VE*100
set y = INFLAIPCD*100
set y2 = (INFLAIPCD*100)^2
set y3 = (INFLAIPCD*100)^3
set w = VSAL*100
set m = M2*100

set dx = VE-VE{1}
set dy = INFLAIPCD-INFLAIPCD{1}
set dv = VSAL-VSAL{1}
*
* 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 = 1970:1
compute send   = 2011:7
*

*
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]   resupup(nvar) resupdn(nvar) resdnup(nvar) resdndn(nvar)
dec rect[frml]     fitud(nvar,2)
dec vect[frml]     tvarf(nvar) empty(nvar) tfrml1(nvar) tfrml2(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    =||y,w,m,x||
*compute shortlabels=||"Infla","Var Sal","Var M","Var TCN"||
*compute longlabels =||"Inflacion mensual","Variacion Salarial mensual","Variacion M","Var TCN"||
*compute lab        =|| "Y", "W", "M", "X"||

compute depvars    =||y,w,x||
compute shortlabels=||"Infla","Var Sal","Var TCN"||
compute longlabels =||"Inflacion mensual","Variacion Salarial mensual","Variacion TCN mensual"||
compute lab        =|| "Y", "W", "X"||
*

*
*  set thresholds, mas, and delays
*
compute d = 1
*
* 15% window
*
dec vector macoeffs1
dec vector macoeffs2
compute malength1 = 6  ; comp thresh1 = 2.55 ;* cpbill, ff
compute malength2 = 6 ; comp thresh2 = 1.25 ;* dmix, ff
*compute malength = 6 ; comp thresh = -3.19 ;* dsd, ff
*
* This generates the actual moving average of the data
*
filter(type=lagging,span=malength1) y / yma
filter(type=lagging,span=malength2) y / xma
*
* This generates an equation (and from it a FRML) to compute the moving
* average of the data
*
dim macoeffs1(malength1)
comp macoeffs1 = %const(1./malength1)
equation(identity,coeffs=macoeffs1) threqn1 yma
# y {0 to malength1-1}
frml(equation=threqn1,identity) thrfrml1

dim macoeffs2(malength2)
comp macoeffs2 = %const(1./malength2)
equation(identity,coeffs=macoeffs2) threqn2 xma
# y {0 to malength2-1}
frml(equation=threqn2,identity) thrfrml2
*
* 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
  6 6 6
  6 6 6
  6 6 6

*************************************************************************
*
* Dummy variables for the two regimes
*

set d1 = (thrfrml1{d}>thresh1)+2*(thrfrml2{d}>thresh2)

*
compute maxlag=d+(malength1-1)
inquire(reglist) rstart<<sstart rend<<send
# yma{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
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 high inflation"
   linreg(smpl=d1==3,equation=eqn(i),frml=fitud(i,1)) * rstart rend resupup(i)
   @lagpolyroots %eqnlagpoly(0,depvars(i))
   disp
   disp
   disp shortlabels(i)+" regression in intermediate inflation"
   linreg(smpl=d1==2,equation=eqn(i),frml=fitud(i,2)) * rstart rend resupdn(i)
   @lagpolyroots %eqnlagpoly(0,depvars(i))
   disp
   disp
   *disp shortlabels(i)+" regression in low inflation, high TCN depreciation"
   *linreg(smpl=d1==1,equation=eqn(i),frml=fitud(i,3)) * rstart rend resdnup(i)
   disp
   disp
   disp shortlabels(i)+" regression in low inflation"
   linreg(smpl=d1==0,equation=eqn(i),frml=fitud(i,4)) * rstart rend resdndn(i)
   @lagpolyroots %eqnlagpoly(0,depvars(i))
   set res(i) rstart rend = %if(d1==3,resupup(i),%if(d1==2,resupdn(i),%if(d1==1,resdnup(i),resdndn(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=vupup)
# resupup
compute supup =%decomp(vupup)
compute siupup=inv(supup)


vcv(matrix=vupdn)
# resupdn
compute supdn =%decomp(vupdn)
compute siupdn=inv(supdn)

*vcv(matrix=vdnup)
*# resdnup
*compute sdnup =%decomp(vdnup)
*compute sidnup=inv(sdnup)

vcv(matrix=vdndn)
# resdndn
compute sdndn =%decomp(vdndn)
compute sidndn=inv(sdndn)
*
* Compute the joint covariance matrix and its factor.
*
vcv(matrix=vsigma)
# res
compute f=%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==3,siupup*%xt(resupup,t),%if(d1==2,siupdn*%xt(resupdn,t),sidndn*%xt(resdndn,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((thrfrml1{d}>thresh1)+(thrfrml2{d}>thresh2),$
        fitud(&i,1)+%dot(%xrow(supup,&i),bootres),$
        %if((thrfrml1{d}<=thresh1)+(thrfrml2{d}>thresh2),$
        fitud(&i,2)+%dot(%xrow(supdn,&i),bootres),$
        fitud(&i,4)+%dot(%xrow(sdndn,&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) thrfrml1 thrfrml2
*
*************************************************************************
*
* 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==3.and.values(1)==3.or.upper==2.and.values(1)==2.or.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
         dofor size = 1.5 3.0
            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(2)
compute klabels=||"1.5 SD","3 SD"||
*
spgraph(vfields=1,hfields=1,footer=$
   "Figure 2. Response of Inflation to Shocks, Conditional on Regime")
do j=3,3
   do i=1,2
      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="Low Inflation Regime, <1.25% average: Response of Inflation",subheader="Shock to "+shortlabels(j))
end do j
spgraph(done)

*@LagPolyRoots(title="Raices estado 1") %eqnlagpoly(estado1,y)
*@LagPolyRoots(title="Raices estado 2") %eqnlagpoly(estado2,y)
*@LagPolyRoots(title="Raices estado 3") %eqnlagpoly(estado3,y)
*@LagPolyRoots(title="Raices estado 4") %eqnlagpoly(estado4,y)

do i =1,nvar
stats resupup(i)
stats resupdn(i)
stats resdndn(i)

corr(qstats,span=1,number=12,method=yule) resupup(i)
corr(qstats,span=1,number=12,method=yule) resupdn(i)
corr(qstats,span=1,number=12,method=yule) resdndn(i)

@archtest(lags=1,form=lm) resupup(i)
@archtest(lags=2,form=lm) resupup(i)
@archtest(lags=3,form=lm) resupup(i)
@archtest(lags=4,form=lm) resupup(i)
@archtest(lags=5,form=lm) resupup(i)
@archtest(lags=6,form=lm) resupup(i)
@archtest(lags=7,form=lm) resupup(i)
@archtest(lags=8,form=lm) resupup(i)
@archtest(lags=9,form=lm) resupup(i)
@archtest(lags=10,form=lm) resupup(i)
@archtest(lags=11,form=lm) resupup(i)
@archtest(lags=12,form=lm) resupup(i)

@archtest(lags=1,form=lm) resupdn(i)
@archtest(lags=2,form=lm) resupdn(i)
@archtest(lags=3,form=lm) resupdn(i)
@archtest(lags=4,form=lm) resupdn(i)
@archtest(lags=5,form=lm) resupdn(i)
@archtest(lags=6,form=lm) resupdn(i)
@archtest(lags=7,form=lm) resupdn(i)
@archtest(lags=8,form=lm) resupdn(i)
@archtest(lags=9,form=lm) resupdn(i)
@archtest(lags=10,form=lm) resupdn(i)
@archtest(lags=11,form=lm) resupdn(i)
@archtest(lags=12,form=lm) resupdn(i)

@archtest(lags=1,form=lm) resdndn(i)
@archtest(lags=2,form=lm) resdndn(i)
@archtest(lags=3,form=lm) resdndn(i)
@archtest(lags=4,form=lm) resdndn(i)
@archtest(lags=5,form=lm) resdndn(i)
@archtest(lags=6,form=lm) resdndn(i)
@archtest(lags=7,form=lm) resdndn(i)
@archtest(lags=8,form=lm) resdndn(i)
@archtest(lags=9,form=lm) resdndn(i)
@archtest(lags=10,form=lm) resdndn(i)
@archtest(lags=11,form=lm) resdndn(i)
@archtest(lags=12,form=lm) resdndn(i)


end do i