*
* @MCVARDoDraws(model=VAR model to analyze,other options)
*
* Computes impulse response functions for a VAR using Monte Carlo
* simulation. This assumes the model is an OLS VAR estimated using the
* ESTIMATE instruction. This uses a Cholesky factorization, though it
* can be modified fairly easily to do any just-identified factorization.
*
* Options:
* MODEL=model to analyze [required]
* STEP=number of response steps[48]
* DRAWS=number of Monte Carlo draws[1000]
* ACCUMULATE=list of variables (by position) to accumulate [none]
*
* FFUNCTION=FUNCTION[RECT](SYMM,MODEL) [not used]
* This must be a FUNCTION which returns an alternate factorization (rather than
* the Cholesky) using the draw for the covariance matrix (first argument in the
* function) and the model (second argument). Note that this is only technically
* correct if the factorization involves a just-identified model. An overidentified
* model requires more complicated methods.
*
* Variables Defined:
* %%responses is a <> vector of RECTANGULAR arrays, each of which
* has dimensions N^2 x <>. The responses of variable i to shock j
* at step s on draw d is at %%responses(d)((j-1)*N+i,s)
*
* Revision Schedule:
* 06/2009 Stripped from montevar.src procedure
* 07/2009 Revised to block rows of %%responses by shocks, not variables.
* 04/2011 Change to use FACTOR rather than FSIGMAD. Easier to explain
* how to modify.
* 08/2014 Added FFUNCTION option.
* 05/2017 Check for %DEFINED on ACCUM option. (Could carry over value from
* previous run that used it).
*
procedure MCVARDoDraws
*
option model model
option integer steps 48
option integer draws 1000
option vect[int] accum
option function[rect](symm,model) ffunction
*
local integer i j nvar
local integer wishdof draw
local rect fxx fwish fsigmad betaols betau betadraw factor
local symm sigmad
local vect ix
local rect[series] impulses
*
declare vect[rect] %%responses
*
if .not.%defined(model) {
disp "Syntax: @MCVARDoDraws(model=model to analyze,other options)"
return
}
compute nvar=%modelsize(model)
*
dim %%responses(draws) impulses(nvar,nvar)
*
compute fxx =%decomp(%xx)
compute fwish =%decomp(inv(%nobs*%sigma))
compute betaols=%modelgetcoeffs(model)
compute wishdof=%nobs-%rows(fxx)
*
* This uses antithetic acceleration to improve convergence of the Monte
* Carlo draws. On the odd values for draws, 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
* inv(S)~Wishart with T-p d.f. (p = number of explanatory variables per
* equation) and covariance matrix inv(T(S-hat)).
*
* A draw for the coefficients is computed by postmultiplying the
* Kronecker product of the decompositions of the sigma draw and X’X
* inverse by a vector of standard Normal draws. For even values of
* "draws", we take the reflection of the original draw around the OLS
* vector.
*
infobox(action=define,progress,lower=1,upper=draws) "Monte Carlo Integration"
do draw = 1,draws
if %clock(draw,2)==1 {
compute sigmad =%ranwisharti(fwish,wishdof)
compute fsigmad =%decomp(sigmad)
compute betau =%ranmvkron(fsigmad,fxx)
compute betadraw=betaols+betau
}
else
compute betadraw=betaols-betau
compute %modelsetcoeffs(model,betadraw)
*
* If FFUNCTION is defined, the function passed by the user is called
* with the just-computed covariance matrix and the model (in case the
* coefficients are needed). In returns, what is expected is a factor
* of the covariance matrix.
*
* Otherwise, the factor is the Cholesky factor.
*
if %defined(ffunction)
compute factor=ffunction(%outerxx(fsigmad),model)
else
compute factor=fsigmad
*
impulse(noprint,model=model,results=impulses,steps=steps,$
factor=factor)
*
* Accumulate the responses as requested
*
if %defined(accum) {
do i=1,%rows(accum)
do j=1,%cols(impulses)
acc impulses(accum(i),j) 1 steps
end do j
end do i
}
*
* Store the impulse responses
*
dim %%responses(draw)(%rows(impulses)*%cols(impulses),steps)
ewise %%responses(draw)(i,j)=ix=%vec(%xt(impulses,j)),ix(i)
infobox(current=draw)
end do draws
infobox(action=remove)
*
* Restore the original coefficients
*
compute %modelsetcoeffs(model,betaols)
end