I have modified your Uhlig(2005)'s sign restriction code with historical decomposition for 4shocks
It works but since I'm a beginner, I might make some mistakes.
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=20
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) baseplusa(nvar)
dec vect[series] upathsaa(nvar) baseplusaa(nvar)
dec vect[series] upathsaaa(nvar) baseplusaaa(nvar)
clear upaths
clear upathsa
clear upathsaa
clear upathsaaa
*************************************************************************************
*
* 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 v3(4)
dec vect v4(3)
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 vect[rect] goodrespaa(1000) goodfevdaa(1000)
declare vect[rect] goodrespaaa(1000) goodfevdaaa(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) goodhistaa(1000) goodhistaaa(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)
*************************************************************************************
*
* 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 - 2nd Shock
*****************************************************
do k=1,KMAX+1
compute ik=%xt(impulses,k)*v2
if ik(1)<0.or.ik(2)<0
branch 105
end do k
@forcedfactor(force=column) sigmad i1~i2 f
compute v3=%ransphere(4)
compute i3=%xsubmat(f,1,6,3,6)*v3
compute v3=inv(p)*i3
*****************************************************
* third Set of Restrictions - 3rd Shock
*****************************************************
do k=1,KMAX+1
compute ik=%xt(impulses,k)*v3
if ik(1)<0
branch 105
end do k
@forcedfactor(force=column) sigmad i1~i2~i3 f
compute v4=%ransphere(3)
compute i4=%xsubmat(f,1,6,4,6)*v4
compute v4=inv(p)*i4
*****************************************************
* Forth Set of Restrictions - 4th Shock
*****************************************************
do k=1,KMAX+1
compute ik=%xt(impulses,k)*v4
if ik(1)>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)
dim goodrespaa(accept)(nstep,nvar) goodfevdaa(accept)(nstep,nvar)
dim goodrespaaa(accept)(nstep,nvar) goodfevdaaa(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)
ewise goodrespaa(accept)(i,j)=(ik=%xt(impulses,i)*v3),ik(j)
ewise goodfevdaa(accept)(i,j)=(ik=(irfsquared(i)*(v3.^2))./(irfsquared(i)*ones)),ik(j)
ewise goodrespaaa(accept)(i,j)=(ik=%xt(impulses,i)*v4),ik(j)
ewise goodfevdaaa(accept)(i,j)=(ik=(irfsquared(i)*(v4.^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) goodhistaa(accept)(nhor,nvar) goodhistaaa(accept)(nhor,nvar)
compute remap=p*v1*tr(v1)*inv(p)
do t=hstart,hend
compute %pt(upaths,t,remap*%xt(u,t))
end do t
compute remap=p*v2*tr(v2)*inv(p)
do t=hstart,hend
compute %pt(upathsa,t,remap*%xt(u,t))
end do t
compute remap=p*v3*tr(v3)*inv(p)
do t=hstart,hend
compute %pt(upathsaa,t,remap*%xt(u,t))
end do t
compute remap=p*v4*tr(v4)*inv(p)
do t=hstart,hend
compute %pt(upathsaaa,t,remap*%xt(u,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
forecast(model=varmodel,from=hstart,to=hend,results=baseplusaa,paths)
# upathsaa
forecast(model=varmodel,from=hstart,to=hend,results=baseplusaaa,paths)
# upathsaaa
*
* 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)-base(j)(i+hstart-1)
ewise goodhistaa(accept)(i,j)=baseplusaa(j)(i+hstart-1)-base(j)(i+hstart-1)
ewise goodhistaaa(accept)(i,j)=baseplusaaa(j)(i+hstart-1)-base(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=4,subhea='Impluse Response')
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
*
do i=1,nvar
compute minlower=maxupper=0.0
smpl 1 accept
do k=1,nstep
set work = goodrespaa(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 = goodrespaaa(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
*
*
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='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
*
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 = goodfevdaa(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='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 = goodfevdaaa(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 1st 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 2nd shocks to "+vl(i)) 3
# respH hstart hend
# upperH / 2
# lowerH / 2
end do i
*
spgraph(done)
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 = goodhistaa(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 3rd shocks to "+vl(i)) 3
# respH hstart hend
# upperH / 2
# lowerH / 2
end do i
*
spgraph(done)
clear upperH lowerH respH
spgraph(vfields=3,hfields=2,hlabel="Figure y. 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 = goodhistaaa(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 4th shocks to "+vl(i)) 3
# respH hstart hend
# upperH / 2
# lowerH / 2
end do i
*
spgraph(done)
