Thank you very much for your suggestions. I believe that I have successfully modified the code for two shocks. Could you please give me some feedbacks whether my code is correct?
Code: Select all
OPEN DATA uhligdata.xls
CALENDAR 1965 1 12
compute missc=1.0e+32
data(format=xlS,org=columns) 1965:1 2003:12 gdpc1 gdpdef cprindex totresns bognonbr fedfunds
set gdpc1 = log(gdpc1)*100.0
set gdpdef = log(gdpdef)*100.0
set cprindex = log(cprindex)*100.0
set totresns = log(totresns)*100.0
set bognonbr = log(bognonbr)*100.0
*
system(model=varmodel)
variables gdpc1 gdpdef cprindex fedfunds bognonbr totresns
lags 1 to 6
deterministic constant
end(system)
estimate(noprint,resid=resids)
dec vect[strings] vl(6)
compute vl=||'Consumer Prices','Real GDP','Interest Rates',$
'Market Index Real Stock Prices','Financial Sector Index Real Stock Prices','Real Exchange Rates'||
compute n1=1000
compute n2=1000
compute nvar=6
compute nstep=60
compute KMAX=5
*************************************************************************************
compute hstart = 1995:1
compute hend = 2003:12
compute nhor = hend-hstart+1
dec vect[series] upaths(nvar) u(nvar) onesteps(nvar) base(nvar) baseplus(nvar)
dec vect[series] upathsa(nvar) ua(nvar) onestepsa(nvar) basea(nvar) baseplusa(nvar)
clear upaths
clear upathsa
*************************************************************************************
*
* This is the standard setup for MC integration of an OLS VAR
*
dec symm s(6,6)
dec vect v1(6) ;* For the unit vector on the 1st draw
dec vect v2(5) ;* For the unit vector on the 2nd draw
dec vect v(6) ;* Working impulse vector
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] goodresp(1000) goodfevd(1000)
declare vect[rect] goodrespa(1000) goodfevda(1000)
declare vector ik a(nvar) ones(nvar)
declare series[rect] irfsquared
compute ones=%const(1.0)
source forcedfactor.src
*************************************************************************************
declare vector hk(nvar)
declare vect[rect] goodhist(1000) goodhista(1000)
*************************************************************************************
*
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)
*************************************************************************************
*
* Use static forecasts over the period for historical decomposition to get the
* residuals. Use dynamic forecasts over the same period to get the base forecast
* series for the historical decomposition.
*
forecast(model=varmodel,from=hstart,to=hend,results=onesteps,errors=u,static)
forecast(model=varmodel,from=hstart,to=hend,results=base)
forecast(model=varmodel,from=hstart,to=hend,results=onestepsa,errors=ua,static)
forecast(model=varmodel,from=hstart,to=hend,results=basea)
*************************************************************************************
*
* This is changed to unit shocks rather than orthogonalized shocks.
*
impulse(noprint,model=varmodel,decomp=p,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
************************************************
* First Set of Restrictions - Unique Impulse Vector
************************************************
compute v1=%ransphere(6)
compute i1=p*v1
do k=1,KMAX+1
compute ik=%xt(impulses,k)*v1
if ik(1)>0.or.ik(2)>0.or.ik(3)<0.or.ik(6)>0
branch 105
end do k
*
* Meets the first restriction
* Draw from the orthogonal complement of i1 (last five columns of the
* factor "f").
*
@forcedfactor(force=column) sigmad i1 f
compute v2=%ransphere(5)
compute i2=%xsubmat(f,1,6,2,6)*v2
compute v2=inv(p)*i2
*****************************************************
* Second Set of Restrictions - Demand Shock
*****************************************************
do k=1,KMAX+1
compute ik=%xt(impulses,k)*v2
if ik(1)<0.or.ik(2)<0.or.ik(3)<0.or.ik(6)<0
branch 105
end do k
*
* Meets both restrictions
*
compute accept=accept+1
dim goodrespa(accept)(nstep,nvar) goodfevda(accept)(nstep,nvar)
dim goodresp(accept)(nstep,nvar) goodfevd(accept)(nstep,nvar)
ewise goodresp(accept)(i,j)=(ik=%xt(impulses,i)*v1),ik(j)
ewise goodfevd(accept)(i,j)=(ik=(irfsquared(i)*(v1.^2))./(irfsquared(i)*ones)),ik(j)
ewise goodrespa(accept)(i,j)=(ik=%xt(impulses,i)*v2),ik(j)
ewise goodfevda(accept)(i,j)=(ik=(irfsquared(i)*(v2.^2))./(irfsquared(i)*ones)),ik(j)
*************************************************************************************
*
* Compute the historical paths of shocks to the system corresponding to the
* <<a>> weights on the Choleski factor.
*
dim goodhist(accept)(nhor,nvar) goodhista(accept)(nhor,nvar)
do t=hstart,hend
compute %pt(upaths,t,p*(i1.*(inv(p)*%xt(u,t))))
end do t
do t=hstart,hend
compute %pt(upathsa,t,p*(i2.*(inv(p)*%xt(ua,t))))
end do t
*
* Forecast over the decomposition period with those added shocks
*
forecast(model=varmodel,from=hstart,to=hend,results=baseplus,paths)
# upaths
forecast(model=varmodel,from=hstart,to=hend,results=baseplusa,paths)
# upathsa
*
* Save the difference between the forecasts with and without and added shocks
* (and shift them into the slots 1,...,nhor) of <<goodhist>>.
*
ewise goodhist(accept)(i,j)=baseplus(j)(i+hstart-1)-base(j)(i+hstart-1)
ewise goodhista(accept)(i,j)=baseplusa(j)(i+hstart-1)-basea(j)(i+hstart-1)
*************************************************************************************
if accept>=1000
break
:105
end do subdraws
if accept>=1000
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
*
clear upper lower resp
*
spgraph(vfields=6,hfields=2,subhea='Monetary Policy Shocks and Exchange Rate Shocks, Thailand')
do i=1,nvar
compute minlower=maxupper=0.0
smpl 1 accept
do k=1,nstep
set work = goodresp(t)(k,i)
compute frac=%fractiles(work,||.16,.84||)
compute lower(k)=frac(1)
compute upper(k)=frac(2)
compute resp(k)=%avg(work)
end do k
*
smpl 1 nstep
graph(pattern,ticks,number=0,picture='##.##',header=+vl(i)) 3
# resp
# upper / 2
# lower / 2
end do i
*
do i=1,nvar
compute minlower=maxupper=0.0
smpl 1 accept
do k=1,nstep
set work = goodrespa(t)(k,i)
compute frac=%fractiles(work,||.16,.84||)
compute lower(k)=frac(1)
compute upper(k)=frac(2)
compute resp(k)=%avg(work)
end do k
*
smpl 1 nstep
graph(pattern,ticks,number=0,picture='##.##',header=+vl(i)) 3
# resp
# upper / 2
# lower / 2
end do i
*
spgraph(done)
smpl
*
clear upper lower resp
*
spgraph(vfields=3,hfields=2,hlabel='Fraction of Variance Explained with Pure-Sign Approach',subhea='Subhea')
do i=1,nvar
compute minlower=maxupper=0.0
smpl 1 accept
do k=1,nstep
set work = goodfevd(t)(k,i)
compute frac=%fractiles(work,||.16,.50,.84||)
compute lower(k)=frac(1)
compute upper(k)=frac(3)
compute resp(k)=frac(2)
display resp(k)
end do k
*
smpl 1 nstep
graph(pattern,ticks,number=0,min=0.0,picture="##.##",header="Fraction Explained for "+vl(i)) 3
# resp
# upper / 2
# lower / 2
end do i
*
spgraph(done)
smpl
*
clear upper lower resp
*
spgraph(vfields=3,hfields=2,hlabel='SECOND Fraction of Variance Explained with Pure-Sign Approach',subhea='Subhea')
do i=1,nvar
compute minlower=maxupper=0.0
smpl 1 accept
do k=1,nstep
set work = goodfevda(t)(k,i)
compute frac=%fractiles(work,||.16,.50,.84||)
compute lower(k)=frac(1)
compute upper(k)=frac(3)
compute resp(k)=frac(2)
display resp(k)
end do k
*
smpl 1 nstep
graph(pattern,ticks,number=0,min=0.0,picture="##.##",header="Fraction Explained for "+vl(i)) 3
# resp
# upper / 2
# lower / 2
end do i
*
spgraph(done)
smpl
*************************************************************************************
declare series upperH lowerH respH
*
spgraph(vfields=3,hfields=2,hlabel="Figure X. Historical Decomposition with Pure-Sign Approach")
do i=1,nvar
compute minlower=maxupper=0.0
smpl 1 accept
clear upperH lowerH respH
do k=hstart,hend
set work 1 accept = goodhist(t)(k-hstart+1,i)
compute frac=%fractiles(work,||.16,.50,.84||)
compute lowerH(k)=frac(1)
compute upperH(k)=frac(3)
compute respH(k)=frac(2)
end do k
* compute maxupper=%max(maxupper,%maxvalue(upper)) ,min=minlower,max=maxupper
* compute minlower=%min(minlower,%minvalue(lower))
*
/*
display 'contributions to variable' +vl(i)
do h=1,nhor
display respH(h)
end do h
*/
smpl hstart hend
graph(pattern,ticks,dates,picture="##.##",header="Contribution of Monetary Policy shocks to "+vl(i)) 3
# respH hstart hend
# upperH / 2
# lowerH / 2
end do i
*
spgraph(done)
smpl
*
clear upperH lowerH respH
*
spgraph(vfields=3,hfields=2,hlabel="Figure X. Historical Decomposition with Pure-Sign Approach")
do i=1,nvar
compute minlower=maxupper=0.0
smpl 1 accept
clear upperH lowerH respH
do k=hstart,hend
set work 1 accept = goodhista(t)(k-hstart+1,i)
compute frac=%fractiles(work,||.16,.50,.84||)
compute lowerH(k)=frac(1)
compute upperH(k)=frac(3)
compute respH(k)=frac(2)
end do k
* compute maxupper=%max(maxupper,%maxvalue(upper)) ,min=minlower,max=maxupper
* compute minlower=%min(minlower,%minvalue(lower))
*
/*
display 'contributions to variable' +vl(i)
do h=1,nhor
display respH(h)
end do h
*/
smpl hstart hend
graph(pattern,ticks,dates,picture="##.##",header="Contribution of Exchange Rate shocks to "+vl(i)) 3
# respH hstart hend
# upperH / 2
# lowerH / 2
end do i
*
spgraph(done)