OPEN DATA "C:\Documents and Settings\n00605626\Desktop\An and Wang\swi1.xls"
CALENDAR(Q) 1973
ALL 2009:02
DATA(FORMAT=XLS,ORG=COLUMNS) 1973:01 2009:02 NER S PPI CPI IP IMP GDP C

seed 33
compute n1=500
compute n2=500
compute nvar=7
compute nstep=20
compute KMAX=2

set NER = -log(NER)
set GDP = log(GDP)
set IMP = log(IMP)
set PPI = log(PPI)
set CPI = log(CPI)
set C = log(C)


set trend = t
set trendsq = t**2
linreg GDP / gap
# constant trend trendsq



diff c / dc
diff ner / dner
diff imp / dimp
diff ppi / dppi
diff cpi /dcpi



system(model=varmodel)
variables  S GAP dNER dIMP dPPI dCPI dc
lags 1 to 6
det constant
end(system)
estimate(outsigma=vmat)


equation(identity, coeffs=||1.0||) seq s1
# s
equation(identity, coeffs=||1.0||) gapeq gap1
# gap
equation(identify,coeffs=||1.0,1.0||) nereq ner
# ner{1} dner

equation(identify,coeffs=||1.0,1.0||) impeq imp
# imp{1} dimp
equation(identify,coeffs=||1.0,1.0||) ppieq ppi
# ppi{1} dppi
equation(identify,coeffs=||1.0,1.0||) cpieq cpi
# cpi{1} dcpi
equation(identify,coeffs=||1.0,1.0||) Ceq c
# C{1} dc



group identmodel  seq gapeq nereq impeq ppieq cpieq Ceq
compute numident=7


*
dec vect[strings] vl(14)
compute vl=||'interest rate','output gap','nominal exchange rate difference','import price difference','producer price difference', 'consumer price difference',$
'consumption difference','interest rate','output gap','nominal exchange rate','import price', 'producer price', 'consumer price', 'consumption'||
*
* n1 is the number of draws from the posterior of the VAR
* n2 is the number of draws from the unit sphere for each draw for the VAR
* nvar is the number of variables
* nstep is the number of IRF steps to compute
* KMAX is the "K" value for the number of steps constrained

dec symm Y(nvar,nvar)
dec vect v1(nvar)
dec vect v2(nvar-1)
dec vect v3(nvar-2)
dec vect v4(nvar-3)
dec vect v(nvar)

*
* This is the standard setup for MC integration of an OLS VAR
*
compute sxx    =%decomp(%xx)
compute svt    =%decomp(inv(%nobs*%sigma))
compute betaols=%modelgetcoeffs(varmodel)
compute ncoef  =%rows(sxx)
compute wishdof=%nobs-ncoef
dec rect ranc(ncoef,nvar)
*
* Most draws are going to get rejected. We allow for up to 1000
* good ones. The variable accept will count the number of accepted
* draws. GOODRESP will be a RECT(nsteps,nvar) at each accepted
* draw.
*
declare vect[rect] goodrespER(5000) goodfevdER(5000)
declare vect[rect] goodrespM(5000) goodfevdM(5000)
declare vect[rect] goodrespD(5000) goodfevdD(5000)
declare vect[rect] goodrespPR(5000) goodfevdPR(5000)

declare vector ik a(nvar) ones(nvar)
declare series[rect] irfsquared
compute ones=%const(1.0)
source foredfactor.src
*
compute accept=0
infobox(action=define,progress,lower=1,upper=n1) 'Monte Carlo Integration'
do draws=1,n1
   *
   * Make a draw from the posterior for the VAR and compute its impulse
   * responses.
   *
   compute sigmad  =%ranwisharti(svt,wishdof)
   compute P   =%decomp(sigmad)
   compute ranc    =%ran(1.0)
   compute betau   =sxx*ranc*tr(P)
   compute betadraw=betaols+betau
   compute %modelsetcoeffs(varmodel,betadraw)
   impulse(noprint,model=varmodel+identmodel,decomp=%identity(nvar),results=impulses, steps=nstep)
gset irfsquared 1 1 = %xt(impulses,t).^2
gset irfsquared 2 nstep = irfsquared{1}+%xt(impulses,t).^2

   * Do the subdraws over the unit sphere. These give the weights on the
   * orthogonal components.
   *
   do subdraws=1,n2
      compute v1=%ransphere(nvar)
      compute i1=p*v1

 * This is to identify the exchange rate shock
      do k=1,KMAX+1
         compute ik=%xt(impulses,k)*v1
         if ik(8)<0.or.ik(10)<0
      branch 105
      end do k

** this is to identify the monetary shock

@forcedfactor(force=column) sigmad i1 f
compute v2=%ransphere(nvar-1)
compute i2=%xsubmat(f,1,nvar,2,nvar)*v2
compute v2=inv(p)*i2
 do k=1,KMAX+1
         compute ik=%xt(impulses,k)*v2
         if ik(8)>0.or.ik(10)<0
      branch 105
      end do k



** this is to identify the positive demand shock
@forcedfactor(force=column) sigmad i1~i2 f
compute v3=%ransphere(nvar-2)
compute i3=%xsubmat(f,1,nvar,3,nvar)*v3
compute v3=inv(p)*i3
 do k=1,KMAX+1
         compute ik=%xt(impulses,k)*v3
         if ik(8)<0.or.ik(10)>0.or.ik(14)<0
      branch 105
      end do k




      compute accept=accept+1
      dim goodrespER(accept)(nstep,nvar+numident) goodfevdER(accept)(nstep,nvar+numident)
      dim goodrespM(accept)(nstep,nvar+numident) goodfevdM(accept)(nstep,nvar+numident)
dim goodrespD(accept)(nstep,nvar+numident) goodfevdD(accept)(nstep,nvar+numident)
dim goodrespPR(accept)(nstep,nvar+numident) goodfevdPR(accept)(nstep,nvar+numident)
      ewise goodrespER(accept)(i,j)=ik=(%xt(impulses,i)*v1),ik(j)
      ewise goodfevdER(accept)(i,j)=(ik=(irfsquared(i)*(v1.^2))./(irfsquared(i)*ones)),ik(j)
     ewise goodrespM(accept)(i,j)=ik=(%xt(impulses,i)*v2),ik(j)
      ewise goodfevdM(accept)(i,j)=(ik=(irfsquared(i)*(v2.^2))./(irfsquared(i)*ones)),ik(j)
 ewise goodrespD(accept)(i,j)=ik=(%xt(impulses,i)*v3),ik(j)
      ewise goodfevdD(accept)(i,j)=(ik=(irfsquared(i)*(v3.^2))./(irfsquared(i)*ones)),ik(j)



      if accept>=5000
         break
   :105
   end do subdraws
   if accept>=5000
      break
   infobox(current=draws)
end do draws
infobox(action=remove)
*
* Post-processing. Graph the mean of the responses along with the 16% and 84%-iles
declare vector[series] IRFer IRFerL IRFerU IrFm IRFmU IRFmL IRFd IRFdU IRFdL IRFpr IRFprU IRFprL
dim IRFer(nvar+numident)
dim IRFerL(nvar+numident)
dim IRFerU(nvar+numident)
dim IRFm(nvar+numident)
dim IRFmL(nvar+numident)
dim IRFmU(nvar+numident)
dim IRFd(nvar+numident)
dim IRFdL(nvar+numident)
dim IRFdU(nvar+numident)
dim IRFpr(nvar+numident)
dim IRFprL(nvar+numident)
dim IRFprU(nvar+numident)


clear upper lower resp
spgraph(vfieds=4,hfields=2,hlabel='impulse response function to exchange rate shock')

do i=nvar+1,nvar+numident
   compute minlower=maxupper=minlowervd=maxuppervd=0.0
   smpl 1 accept
   do k=1,nstep
      set work = goodrespER(t)(k,i)
      compute frac=%fractiles(work,||.16,.84,.5||)
      compute lower(k)=frac(1)
      compute upper(k)=frac(2)
      compute resp(k)=frac(3)


   end do k


   smpl 1 nstep
   graph(ticks,number=0,picture='##.##',header=vl(i)) 3
   # resp
   # upper / 2
   # lower / 2
set IRFer(i) = resp
set  IRFerL(i) = lower
set irferu(i) = upper

end do i
spgraph(done)

** draw the impulse response to the monetary shock

clear upper lower resp
spgraph(vfieds=4,hfields=2,hlabel='impulse response function to monetary shock')

do i=nvar+1,nvar+numident
   compute minlower=maxupper=minlowervd=maxuppervd=0.0
   smpl 1 accept
   do k=1,nstep
      set work = goodrespM(t)(k,i)
      compute frac=%fractiles(work,||.16,.84,.5||)
      compute lower(k)=frac(1)
      compute upper(k)=frac(2)
      compute resp(k)=frac(3)


   end do k


   smpl 1 nstep
   graph(ticks,number=0,picture='##.##',header=vl(i)) 3
   # resp
   # upper / 2
   # lower / 2
set IRFm(i) = resp
set  IRFmL(i) = lower
set irfmu(i) = upper

end do i
spgraph(done)

** draw the impulse response to the demand shock
clear upper lower resp
spgraph(vfieds=4,hfields=2,hlabel='impulse response function to demand shock')

do i=nvar+1,nvar+numident
   compute minlower=maxupper=minlowervd=maxuppervd=0.0
   smpl 1 accept
   do k=1,nstep
      set work = goodrespD(t)(k,i)
      compute frac=%fractiles(work,||.16,.84,.5||)
      compute lower(k)=frac(1)
      compute upper(k)=frac(2)
      compute resp(k)=frac(3)


   end do k


   smpl 1 nstep
   graph(ticks,number=0,picture='##.##',header=vl(i)) 3
   # resp
   # upper / 2
   # lower / 2
set IRFd(i) = resp
set  IRFdL(i) = lower
set irfdu(i) = upper

end do i
spgraph(done)


**get PT to ER shock
clear PPIPTer PPIPTerU PPIPTerL PPIPTera PPIPTeraU PPIPTeraL
set PPIPerT = IRFer(12)(t)/IRFer(10)(1)
set PPIPTerU = IRFerU(12)(t)/IRFerU(10)(1)
set PPIPTerL = IRFerL(12)(t)/IRFerL(10)(1)
set PPIPTera = IRFer(12)(t)/IRFer(10)(t)
set PPIPTeraU = IRFerU(12)(t)/IRFerU(10)(t)
set PPIPTeraL = IRFerL(12)(t)/IRFerL(10)(t)
clear CPIPTer CPIPTerU CPIPTerL CPIPTera CPIPTeraU CPIPTeraL
set CPIPTer = IRFer(13)(t)/IRFer(10)(1)
set CPIPTerU = IRFerU(13)(t)/IRFerU(10)(1)
set cPIPTerL = IRFerL(13)(t)/IRFerL(10)(1)
set CPIPTera = IRFer(13)(t)/IRFer(10)(t)
set CPIPTeraU = IRFerU(13)(t)/IRFerU(10)(t)
set CPIPTeraL = IRFerL(13)(t)/IRFerL(10)(t)
clear IMPPTer IMPPTerU IMPPPTerL IMPPTera IMPPTeraU IMPPTeraL
set IMPPTer = IRFer(11)(t)/IRFer(10)(1)
set IMPPTerU = IRFerU(11)(t)/IRFerU(10)(1)
set IMPPTerL = IRFerL(11)(t)/IRFerL(10)(1)
set IMPPTera = IRFer(11)(t)/IRFer(10)(t)
set IMPPTeraU = IRFerU(11)(t)/IRFerU(10)(t)
set IMPPTeraL = IRFerL(11)(t)/IRFerL(10)(t)

** get PT ratio to monetary shock
clear PPIPTm PPIPTmU PPIPTmL PPIPTma PPIPTmaU PPIPTmaL
set PPIPTm = IRFm(12)(t)/IRFm(10)(1)
set PPIPTmU = IRFmU(12)(t)/IRFmU(10)(1)
set PPIPTmL = IRFmL(12)(t)/IRFmL(10)(1)
set PPIPTma = IRFm(12)(t)/IRFm(10)(t)
set PPIPTmaU = IRFmU(12)(t)/IRFmU(10)(t)
set PPIPTmaL = IRFmL(12)(t)/IRFmL(10)(t)
clear CPIPTm CPIPTmU CPIPTmL CPIPTma CPIPTmaU CPIPTmaL
set CPIPTm = IRFm(13)(t)/IRFm(10)(1)
set CPIPTmU = IRFmU(13)(t)/IRFmU(10)(1)
set cPIPTmL = IRFmL(13)(t)/IRFmL(10)(1)
set CPIPTma = IRFm(13)(t)/IRFm(10)(t)
set CPIPTmaU = IRFmU(13)(t)/IRFmU(10)(t)
set CPIPTmaL = IRFmL(13)(t)/IRFmL(10)(t)
clear IMPPT IMPPTU IMPPPTL IMPPTa IMPPTaU IMPPTaL
set IMPPTm = IRFm(11)(t)/IRFm(10)(1)
set IMPPTmU = IRFmU(11)(t)/IRFmU(10)(1)
set IMPPTmL = IRFmL(11)(t)/IRFmL(10)(1)
set IMPPTma = IRFm(11)(t)/IRFm(10)(t)
set IMPPTmaU = IRFmU(11)(t)/IRFmU(10)(t)
set IMPPTmaL = IRFmL(11)(t)/IRFmL(10)(t)

** PT ratio to demand shock
clear PPIPTd PPIPTdU PPIPTdL PPIPTda PPIPTdaU PPIPTdaL
set PPIPTd = IRFd(12)(t)/IRFd(10)(1)
set PPIPTdU = IRFdU(12)(t)/IRFdU(10)(1)
set PPIPTdL = IRFdL(12)(t)/IRFdL(10)(1)
set PPIPTda = IRFd(12)(t)/IRFd(10)(t)
set PPIPTdaU = IRFdU(12)(t)/IRFdU(10)(t)
set PPIPTdaL = IRFdL(12)(t)/IRFdL(10)(t)
clear CPIPTd CPIPTdU CPIPTdL CPIPTda CPIPTdaU CPIPTdaL
set CPIPTd = IRFd(13)(t)/IRFd(10)(1)
set CPIPTdU = IRFdU(13)(t)/IRFdU(10)(1)
set cPIPTdL = IRFdL(13)(t)/IRFdL(10)(1)
set CPIPTda = IRFd(13)(t)/IRFd(10)(t)
set CPIPTdaU = IRFdU(13)(t)/IRFdU(10)(t)
set CPIPTdaL = IRFdL(13)(t)/IRFdL(10)(t)
clear IMPPTd IMPPTdU IMPPPTdL IMPPTda IMPPTdaU IMPPTdaL
set IMPPTd = IRFd(11)(t)/IRFd(10)(1)
set IMPPTdU = IRFdU(11)(t)/IRFdU(10)(1)
set IMPPTdL = IRFdL(11)(t)/IRFdL(10)(1)
set IMPPTda = IRFd(11)(t)/IRFd(10)(t)
set IMPPTdaU = IRFdU(11)(t)/IRFdU(10)(t)
set IMPPTdaL = IRFdL(11)(t)/IRFdL(10)(t)




smpl 1 nstep
spgraph(vfieds=6,hfields=4,hlabel='PT ratios')
 graph(ticks, number=0, picture='*.##', header='IMP PT to ER shock') 3
 # IMPPTer
 # IMPPTerU / 2
 # IMPpterl / 2


graph(ticks, number=0, picture='*.##', header='IMP PT to ER shock ALT') 3
 # IMPPTera
 # IMPPTeraU / 2
 # IMPpteral / 2




graph(ticks, number=0, picture='*.##', header='PPI PT to ER shock') 3
 # PPIPTer
 # PPIPTerU / 2
 # PPIpterl / 2

graph(ticks, number=0, picture='*.##', header='PPI PT to ER shock ALT') 3
 # PPIPTera
 # PPIPTeraU / 2
 # PPIpteral / 2


graph(ticks, number=0, picture='*.##', header='CPI PT to ER shock') 3
 # PPIPTer
 # PPIPTeraU / 2
 # PPIpteral / 2


 graph(ticks, number=0, picture='*.##', header='CPI PT to ER shock ALT') 3
 # CPIPTer
 # CPIPTeraU / 2
 # CPIpteral / 2

 graph(ticks, number=0, picture='*.##', header='IMP PT to Monetary shock') 3
 # IMPPTm
 # IMPPTmU / 2
 # IMPptml / 2


graph(ticks, number=0, picture='*.##', header='IMP PT to Monetary shock ALT') 3
 # IMPPTma
 # IMPPTmaU / 2
 # IMPptmal / 2




graph(ticks, number=0, picture='*.##', header='PPI PT to Monetary shock') 3
 # PPIPTm
 # PPIPTmU / 2
 # PPIptml / 2

graph(ticks, number=0, picture='*.##', header='PPI PT to Monetary shock ALT') 3
 # PPIPTma
 # PPIPTmaU / 2
 # PPIptmal / 2


graph(ticks, number=0, picture='*.##', header='CPI PT to Monetary shock') 3
 # PPIPTm
 # PPIPTmU / 2
 # PPIptml / 2


 graph(ticks, number=0, picture='*.##', header='CPI PT to Monetary shock ALT') 3
 # CPIPTm
 # CPIPTmaU / 2
 # CPIptmal / 2

graph(ticks, number=0, picture='*.##', header='IMP PT to Demand shock') 3
 # IMPPTd
 # IMPPTdU / 2
 # IMPptdl / 2


graph(ticks, number=0, picture='*.##', header='IMP PT to Demand shock ALT') 3
 # IMPPTda
 # IMPPTdaU / 2
 # IMPptdal / 2


graph(ticks, number=0, picture='*.##', header='PPI PT to Demand shock') 3
 # PPIPTd
 # PPIPTdU / 2
 # PPIptdl / 2

graph(ticks, number=0, picture='*.##', header='PPI PT to Demand shock ALT') 3
 # PPIPTda
 # PPIPTdaU / 2
 # PPIptdal / 2


graph(ticks, number=0, picture='*.##', header='CPI PT to Demand shock') 3
 # CPIPTd
 # CPIPTdU / 2
 # CPIptdl / 2


 graph(ticks, number=0, picture='*.##', header='CPI PT to Demand shock ALT') 3
 # CPIPTd
 # CPIPTdaU / 2
 # CPIptdal / 2



spgraph(done)


report(action=define, hlabels=||'lags', 'IMPPT', 'PPIPT', 'CPIPT'||)
do i=1, 16
report(row=new, atcol=1) i IMPPTer(i) PPIPTer(i) CPIPTer(i)
end do i
report(action=show, window='pass-through ratio under ER shock')

report(action=define, hlabels=||'lags', 'IMPPT', 'PPIPT', 'CPIPT'||)
do i=1, 16
report(row=new, atcol=1) i IMPPTera(i) PPIPTera(i) CPIPTera(i)
end do i
report(action=show, window='pass-through ratio to ER shock ALT')

report(action=define, hlabels=||'lags', 'IMPPT', 'PPIPT', 'CPIPT'||)
do i=1, 16
report(row=new, atcol=1) i IMPPTm(i) PPIPTm(i) CPIPT(mi)
end do i
report(action=show, window='pass-through ratio to monetary shock')

report(action=define, hlabels=||'lags', 'IMPPT', 'PPIPT', 'CPIPT'||)
do i=1, 16
report(row=new, atcol=1) i IMPPTma(i) PPIPTma(i) CPIPTma(i)
end do i
report(action=show, window='pass-through ratio to monetary shock ALT')

report(action=define, hlabels=||'lags', 'IMPPT', 'PPIPT', 'CPIPT'||)
do i=1, 16
report(row=new, atcol=1) i IMPPTd(i) PPIPTd(i) CPIPTd(i)
end do i
report(action=show, window='pass-through ratio to demand shock')

report(action=define, hlabels=||'lags', 'IMPPT', 'PPIPT', 'CPIPT'||)
do i=1, 16
report(row=new, atcol=1) i IMPPTda(i) PPIPTda(i) CPIPTda(i)
end do i
report(action=show, window='pass-through ratio to demand shock ALT')






display accept