******Monte Carlo Integration******
*
OPEN DATA "C:\Users\jfw43\Dissertation\Eviews and Data Files\Five Variable VAR.RAT"
DATA(FORMAT=RATS) 1959:01 2014:03 D2LRGDP DLCONS DR DS INFLATION
*
compute nvar = 5
compute nlags = 1
compute nsteps = 100
compute ndraws = 10000
*
compute [vect[strings]] shocklabels=||"Output","Savings","Inflation","Interest Rate","Consumption"||
compute [vect[strings]] varlabels=||"Real GDP","Savings Rate","Inflation","Nominal Interest Rate","Consumption"||
*
system(model=fivevar)
var d2lrgdp ds inflation dr dlcons
lags 1 to nlags
det
end(system)
estimate(noprint,resids=resids)
compute vsigma=%sigma
*
dec rect lr(5,5) sr(5,5)
input lr
. 0 . . .
. . . . .
. . . . .
. . . . .
. . . . .
input sr
. . 0 0 0
. . 0 0 0
. . . 0 0
. . . . 0
. . . . .
dec vect[strings] shocklabels varlabels
compute shocklabels=||"Output","Savings","Inflation","Interest rate","Consumption"||
compute varlabels=||"Real GDP","Savings Rate","CPI Inflation","Nominal Interest rate","Consumption"||
*
@ShortAndLong(lr=lr,sr=sr,masum=inv(%varlagsums)) %sigma f
compute fd=%xdiag(f)
ewise f(i,j)=f(i,j)/fd(j)
*
@varirf(noprint,model=fivevar,steps=nsteps,factor=f,page=byshock,accum=||1,2,4,5||,$
  shocks=shocklabels,varlabels=varlabels)
*
procedure SRLRDoDraws
*
option model     model
option integer   ndraws   10000
option integer   nsteps   100
option vect[int] accum   ||1,2,4,5||
option rect      lr
option rect      sr
*
local integer nvar
local rect    fxx fwish fsigmad factor
local integer wishdof
local rect    betaols betau betadraw
local vect    ix
local symm    sigmad
local integer i j draw
*
if .not.%defined(fivevar) {
   disp "###SRLRDoDraws(MODEL=fivevar,other options)"
   return
}
compute nvar=%modelsize(fivevar)

declare vect[rect] %%responses
dim %%responses(ndraws)
*
compute fxx    =%decomp(%xx)
compute fwish  =%decomp(inv(%nobs*%sigma))
compute wishdof=%nobs-%nreg
compute betaols=%modelgetcoeffs(fivevar)
*
local rect[series] impulses(nvar,nvar)
*
*
infobox(action=define, progress, lower=1, upper=ndraws) "Monte Carlo Integration"
do draw=1, ndraws
compute fxx    =%decomp(%xx)
compute fwish  =%decomp(inv(%nobs*%sigma))
compute wishdof=%nobs-%nreg
compute betaols=%modelgetcoeffs(fivevar)
*
local rect[series] impulses(nvar,nvar)
declare vect[rect]   %%responses(ndraws)
*
infobox(action=define,progress,lower=1,upper=ndraws) "Monte Carlo Integration"
do draw=1,ndraws
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)
*
@ShortAndLong(lr=lr,sr=sr,masum=inv(%modellagsums(fivevar))) sigmad factor
impulse(noprint,model=fivevar,results=impulses,steps=nsteps)
	*
	do i=1, %rows(accum)
    	do j=1,nvar
        acc impulses(accum(i),j) 1 nsteps
    	end do j
	end do i
	*
	dim %%responses(draw)(nvar*nvar,nsteps)
	ewise %%responses(draw)(i,j)=ix=%vec(%xt(impulses,j)),ix(i)
	infobox (current=draw)
end do draw
infobox(action=remove)
*
@SRLRDoDraws(model=fivevar,accum=||1,2,4,5||,steps=nsteps)
*@MCProcessIRF(model=fivevar,lower=lower,upper=upper,irf=irf)
@MCGraphIRF(model=fivevar,shocklabels=shocklabels,varlabels=varlabels,page=byshock)