OPEN DATA "C:******\Daten.xlsx"
CALENDAR(M) 1992:1
Allocate 2017:02
DATA(FORMAT=XLSX,ORG=COLUMNS,SHEET="Data") 1992:01 2017:02 GER_GDP GER_CPI Euribor shadow GER_EPU CISS GER_SMV $
GER_Credit

* Estimation of threshold
*
* Control parameters
*
* Fraction of threshold values excluded at each end
*
compute pi=.20
*
* Number of VAR lags
*Moving average terms in smoothing threshold variable
* Threshold delay
compute maxlag=4
compute malength = 2
compute d = 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 = 1997:1
compute send   = 2016:10
*
set lGER_GDP = log(GER_GDP)
set lGER_CPI = log(GER_CPI)
set lGER_Credit = log(GER_Credit)
set lGER_EPU = log(GER_EPU)
set trend = t
system(model=varmodel)
variables lGER_GDP lGER_CPI shadow GER_EPU
lags 1 to maxlag
det constant trend
end(system)
*
compute rstart=sstart+maxlag,rend=send
estimate(noprint,resids=resids) rstart rend

compute nvar =%nvar




****RESIDUAL ANALYSE
**Autocorrelation
****LB-test
display ""
display "LB Autokorrelationstest"
dofor i =  1 to nvar
dofor lag = 1 2 3 4 5 6
@regcorrs(noprint,dfc=%narma,number=lag,qstats,nograph, method=burg,title="AR(4) model diagnostics;RGDP") resids(i)
disp  "Gleichung" i "Lag" lag  #.### %Qsignif
end dofor
end dofor

****WESTCHOTEST / Robust to Heteroskedasticity in the Residuals
display ""
display "West CHO Autokorrelationstest"
dofor i = 1 to nvar
dofor lag = 1 2 3 4 5 6
@WestChoTest(number=lag, noprint) resids(i)
disp  "Gleichung" i "Lag" lag  #.### %Qsignif
end dofor
end dofor

********Normalverteilung (JBTest)
dofor i = 1 to nvar
stats(noprint) resids(i)
disp  "Gleichung" i #.### %JBSIGNIF
end dofor


compute loglr=%logl
display %logl
filter(type=lagging,span=malength) lGER_EPU / lGER_EPUthr
*
* Get a sorted copy of the threshold variable
*
set copy rstart rend = lGER_EPUthr{d}
order copy rstart rend
*

compute piskip=fix(pi*%nobs)+(maxlag*nvar+1)
compute pistart=rstart+piskip,piend=rend-piskip
*
compute bestlogl=loglr
*
do pientry=pistart,piend
   compute thresh=copy(pientry)
   *
   * This allows for the covariance matrix to differ between regimes
   *
   sweep(var=hetero,group=lGER_EPUthr{d}<thresh) rstart rend
   # %modeldepvars(varmodel)
   # %rlfromeqn(%modeleqn(varmodel,1))
   *
   * If this is best so far, save it
   *
   if %logl>bestlogl
       compute bestlogl=%logl,bestthresh=thresh
end do pientry
disp "Best Threshold Value" bestthresh
disp "LR Statistic" 2.0*(bestlogl-loglr)
disp bestlogl

*
* The degrees of freedom of the test count both the number of VAR
* parameters and the number of elements in the second estimated
* covariance matrix.
disp "Degrees of freedom" %nregsystem/2+%nvar*(%nvar+1)/2

graph(style=lines,nodates,grid=(t==pistart).or.(t==piend),VGRID=bestthresh) 1
# copy

**set shading = thrfrml{d}>thresh

*graph(style=lines,shad=shading) 1
*# lGER_EPU

********************************************Schaetzung der IRF
* 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    =4
comp horizon =36
*
* Number of bootstrap replications used in computing GIRF's
*
comp nkrep   =50
*
* Change to upper=0 to get the lower 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.

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    =||lGER_GDP,lGER_CPI,shadow,lGER_EPU||
compute shortlabels=||"GDP","CPI","SR","PolUnc"||
compute longlabels =||"GDP","Consumer Price Index","Shadow Rate","Uncertainty"||
compute lab        =|| "Y", "CPI", "SR", "UNC" ||
*
*  set thresholds, mas, and delays
*
compute d = 1
*
* 15% window
*
dec vector macoeffs
compute malength = 2  ; comp thresh = 4.72038 ;* PolUnc, SR
*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
*
* This generates an equation (and from it a FRML) to compute the moving
* average of the data
* /
dim macoeffs(malength)
comp macoeffs = %const(1./malength)
equation(identity,coeffs=macoeffs) threqn lGER_EPUthr
# lGER_EPU{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
 4 4 4 4
 4 4 4 4
 4 4 4 4
 4 4 4 4
*************************************************************************
*
* Dummy variables for the two regimes
*
set d1 = thrfrml{d}>thresh
set d2 = 1.-d1
*
compute maxlag=d+(malength-1)
inquire(reglist) rstart<<sstart rend<<send
# lGER_EPUthr{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 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) 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"||
*
declare string patch
if upper = 1
compute patch = "High"
else
compute patch = "Low"
end if
display patch

spgraph(vfields=2,hfields=2)
compute j = 3
do sh=1,4
   do i=1,4
      set graphs(i) 1 horizon = girfs(i)(sh,j)(t+wstart-1)
   end do i
   graph(series=graphs,number=0,key=upright,klabels=klabels,$
     header= patch+" Uncertainty: Response of "+shortlabels(sh),subheader="Shock to "+shortlabels(j))
end do sh
spgraph(done)