Code:
Code: Select all
system(model=varmodel)
variables gdp gdpdef cprindex fedfunds bognonbr totresns X
lags 1 to 2
end(system)
estimate(noprint)
*
source uhligfuncs.src
*
dec vect[strings] vl(8)
compute vl=||"real GDP","GDP price defl","Comm. Price Ind.",$
"Fed Funds Rate","Nonborr. Reserv.","Total Reserves","X"||
*
* 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
*
compute n2=50000
compute nkeep=1000
compute nvar=7
compute nstep=20
compute KMIN=1
compute KMAX=2
*
* Most draws are going to get rejected. We allow for up to nkeep 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(nkeep)
declare vector ik a(nvar) ones(nvar)
*
impulse(noprint,model=varmodel,factor=%decomp(%sigma),results=impulses,steps=nstep)
compute accept=0
infobox(action=define,progress,lower=1,upper=n2) $
"Monte Carlo Integration"
do subdraws=1,n2
*
* Do the subdraws over the unit sphere. These give the weights on the
* orthogonal components.
*
compute a=%ransphere(nvar)
*
* Check that the responses have the correct signs for steps 1 to
* KMAX+1 (+1 because in the paper, the steps used 0-based
* subscripts, rather than the 1 based subscripts used by RATS).
*
if UhligAccept(a,KMIN,KMAX+1,||+4,-3,-2,-5,-1||)==0
goto reject
*
* This is an accepted draw. Copy the information out. If we have
* enough good ones, drop out of the overall loop.
*
compute accept=accept+1
dim goodresp(accept)(nstep,nvar)
ewise goodresp(accept)(i,j)=ik=%xt(impulses,i)*a,ik(j)
infobox(current=subdraws)
if accept>=nkeep
break
:reject
end do subdraws
infobox(action=remove)
*
* Compute Fry-Pagan minimizers. This uses the interquartile range of the
* responses at each horizon rather than standard error for rescaling the
* difference between a test response and the median.
*
dec vect[series] targets(nvar)
dec vect[series] rescale(nvar)
dec vect[series] lower(nvar) upper(nvar)
do i=1,nvar
set targets(i) 1 nstep = 0.0
set rescale(i) 1 nstep = 0.0
set lower(i) 1 nstep = 0.0
set upper(i) 1 nstep = 0.0
do k=1,nstep
set work 1 accept = goodresp(t)(k,i)
compute frac=%fractiles(work,||.25,.50,.75,.16,.84||)
compute targets(i)(k)=frac(2)
compute rescale(i)(k)=1.0/(frac(3)-frac(1))
compute lower(i)(k)=frac(4)
compute upper(i)(k)=frac(5)
end do k
end do i
*
* This does the optimization. It parameterizes the unit circle using the
* stereographic mapping. The function being minimized is the sum of
* squared standardized gaps between the impulse responses given the test
* rotation and median sign-restricted response where the standardization
* is 1/IQR.
*
dec vect g(nvar-1)
compute g=%fill(nvar-1,1,1.0)
nonlin g
dec real gap
find(trace) min gap
compute a=%stereo(g)
compute gap=0.0
do k=1,nstep
compute gap=gap+%normsqr(%xt(rescale,k).*(%xt(targets,k)-%xt(impulses,k)*a))
end do k
end find
*
* Graph with 16-84% bounds. Note that there is no reason the Fry-Pagan
* minimizing IRF has to be within the bounds at all horizons since it's
* chosen to minimize the sum across (in this case) 360 squared
* deviations from the median responses.
*
set resp 1 nstep = 0.0
spgraph(vfields=4,hfields=2,$
hlabel="Impulse Responses")
do i=1,nvar
set resp 1 nstep = ik=%xt(impulses,t)*a,ik(i)
graph(ticks,number=0,picture="##.##",header="Impulse Responses for "+vl(i)) 3
# resp 1 nstep
# lower(i) 1 nstep 2
# upper(i) 1 nstep 2
end do i
*
spgraph(done)