To follow up on the first entry in this thread, I am posting below my coding attempt at replicating Anzuini, Lombardi and Pagano (2013). Sepcifically, the code below modifies MONTEVAR.RPF (After reading a number of Tom's responses on the forum) to obtain the responses to a
negative 100 basis point decrease in the Federal funds rate (FFR) (as the authors do) along with their Confidence Intervals. The MONTEVAR.RPF is modified using:
Code: Select all
compute fsigma2 = fsigma*inv(%diag(%xdiag(fsigma)))
to obtain a one unit shock (1%) to the FFR and
Code: Select all
compute flipper = ||-1.0,1.0,1.0,1.0,1.0||
as well as:
Code: Select all
compute fsigma3 = fsigma2*%diag(flipper)
to obtain the negative FFR shock. Given that I reconstructed the authors' data, the magnitudes of the responses are off (so is the shape for the CPI and the commodity price index). The responses I obtain are attached. I am posting this just in case it is useful for those of you trying to work with the paper's VAR model.
Code: Select all
*
* MONTEVAR.RPF for Anzuini, Lombardi and Pagano (2013)
* RATS Version 8, User's Guide, Example 16.2
* Monte Carlo Integration of Impulse Responses
*
compute lags=12 ;*Number of lags
compute nstep=50 ;*Number of response steps
compute ndraws=10000 ;*Number of keeper draws
*
calendar(m) 1970:01
allocate 2008:12
OPEN DATA "C:\Replications VAR Commodities\Anzuini, Lombardi and Pagano (2013)\ALP_Data.xlsx"
DATA(FORMAT=XLSX,ORG=COLUMNS) 1970:01 2008:12 CPI FFR IP M2 CRB_BLS Oil
*
**Obtain log levels of variables
set lm2 = log(m2)
set lcpi = log(cpi)
set lcrb = log(crb_bls)
set lip = log(ip)
*
system(model=varmodel)
variables ffr lip lm2 lcpi lcrb
lags 1 to lags
det constant
end(system)
******************************************************************
estimate(noprint)
compute nvar =%nvar
compute fxx =%decomp(%xx)
compute fwish =%decomp(inv(%nobs*%sigma))
compute wishdof=%nobs-%nreg
compute betaols=%modelgetcoeffs(varmodel)
*
*
dec frml[rect] afrml
nonlin g12 g15 g21 g23 g24 g34 g51 g52 g53 g54
frml afrml = || 1.0, g12, 0.0, 0.0, g15|$
g21, 1.0, g23, g24, 0.0|$
0.0, 0.0, 1.0, g34, 0.0|$
0.0, 0.0, 0.0, 1.0, 0.0|$
g51, g52, g53, g54, 1.0||
cvmodel(a=afrml,method=bfgs,iters=1000,factor=fsigma) %sigma
declare vect[rect] %%responses(ndraws)
declare rect[series] impulses(nvar,nvar)
infobox(action=define,progress,lower=1,upper=ndraws) "Monte Carlo Integration"
do draw=1,ndraws
*
* On the odd values for <<draw>>, a draw is made from the inverse Wishart
* distribution for the covariance matrix. This assumes use of the
* Jeffrey's prior |S|^-(n+1)/2 where n is the number of equations in
* the VAR. The posterior for S with that prior is inverse Wishart with
* T-p d.f. (p = number of explanatory variables per equation) and
* covariance matrix inv(T(S-hat)).
*
* Given the draw for S, a draw is made for the coefficients by adding
* the mean from the least squares estimate to a draw from a
* multivariate Normal with (factor of) covariance matrix as the
* Kroneker product of the factor of the draw for S and a factor of
* the X'X^-1 from OLS.
*
* On even draws, the S is kept at the previous value, and the
* coefficient draw is reflected through the mean.
*
if %clock(draw,2)==1 {
compute sigmad =%ranwisharti(fwish,wishdof)
compute fsigma =%decomp(sigmad)
compute betau =%ranmvkron(fsigma,fxx)
compute betadraw=betaols+betau
}
else
compute betadraw=betaols-betau
*
* Push the draw for the coefficient back into the model.
*
compute %modelsetcoeffs(varmodel,betadraw)
compute fsigma2 = fsigma*inv(%diag(%xdiag(fsigma)))
compute flipper = ||-1.0,1.0,1.0,1.0,1.0||
compute fsigma3 = fsigma2*%diag(flipper)
impulse(noprint,model=varmodel,factor=fsigma3,$
results=impulses,steps=nstep)
*
* Save the impulse responses
*
dim %%responses(draw)(nvar*nvar,nstep)
ewise %%responses(draw)(i,j)=ix=%vec(%xt(impulses,j)),ix(i)
infobox(current=draw)
end do draw
infobox(action=remove)
*
@mcgraphirf(model=varmodel,page=byshock)