OPEN DATA "C:\Users\icbc\Desktop\新建文件夹 (2)\换成实际有效汇率指数\199602开始差分x12.csv"
CALENDAR(M) 1996:2
DATA(FORMAT=PRN,ORG=COLUMNS) 1996:02 2016:04 brent wti stockreturn exchangerate monetary output cpi t ect


comp nvar    =5
comp horizon =20

comp nkrep   =300

comp upper   =2
comp [vector] shocksizes=||1.0,2.0||
comp [vector] shocksigns=||1.0,-1.0||

compute sstart = 1996:2
compute send   = 2016:4

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) resin(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    =||wti,exchangerate,monetary,output,cpi||
compute shortlabels=||"wti","exchangerate","monetary","output","cpi"||
compute longlabels =||"oil wti","exchange rate","monetary supply","industrial output","inflation cpi"||
compute lab        =|| "W", "E", "M","O","C"||

compute d = 1
dec vector macoeffs
compute thresh1 =0.01679
compute thresh2 =-0.04757
set dff1 = wti{1}
set dff2 = wti{1}
frml(identity) threshold1 dff1 = wti{1}
frml(identity) threshold2 dff2 = wti{1}
input laglengths
  1 1 1 1 1
  1 1 1 1 1
  1 1 1 1 1
  1 1 1 1 1
  1 1 1 1 1
set d1 = (threshold2>thresh2)+(threshold1>thresh1)


compute maxlag=d

equation eqn(1) depvars(1)
   # constant wti{1} exchangerate{1} monetary{1} output{1} cpi{1} 
equation eqn(2) depvars(2)
   # constant wti{1} exchangerate{1} monetary{1} output{1} cpi{1} 
equation eqn(3) depvars(3)
   # constant wti{1} exchangerate{1} monetary{1} output{1} cpi{1} 
equation eqn(5) depvars(5)
   # constant wti{1} exchangerate{1} monetary{1} output{1} cpi{1} 
 equation eqn(4) depvars(4)
   # constant wti{1} exchangerate{1} monetary{1} output{1} cpi{1} ect 
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==2,equation=eqn(i),frml=fitud(i,1)) * rstart rend resup(i)
   disp
   disp
   disp shortlabels(i)+" regression in lower regime"
   linreg(smpl=d1==0,equation=eqn(i),frml=fitud(i,2)) * rstart rend resdn(i)
   disp
   disp shortlabels(i)+" regression in intermediate regime"
   linreg(smpl=d1==1,equation=eqn(i),frml=fitud(i,3)) * rstart rend resin(i)
set res(i) rstart rend = %if(d1==2,resup(i),%if(d1==0,resdn(i),resin(i)))

end do i

vcv(matrix=vup)
# resup
compute sup =%decomp(vup)
compute siup=inv(sup)

vcv(matrix=vdn)
# resdn
compute sdn =%decomp(vdn)
compute sidn=inv(sdn)
vcv(matrix=vin)
# resin
compute sen =%decomp(vin)
compute seen=inv(sen)

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==2,siup*%xt(resup,t),%if(d1==0,sidn*%xt(resdn,t),seen*%xt(resin,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(threshold1>thresh1,$
        fitud(&i,1)+%dot(%xrow(sup,&i),bootres),$
        %if((threshold1<=thresh1)+(threshold2>=thresh2),fitud(&i,2)+%dot(%xrow(sdn,&i),bootres),$
        fitud(i,3)+%dot(%xrow(sen,&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) threshold1 threshold2
*************************************************************************
*
panel(group=d1{0},id=values,identries=entries) d1 rstart rend
{
if 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 -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
*
*************************************************************************