The RATS Software Forum

Dear RATS users,

In the program below, I’m estimating the error bands
of the impulse responses in a Near-SVAR model. The
MONTESUR.PRG (the most recent one) estimates the
Near-VAR model, it but uses the Choleski
decomposition, not the Structural decomposition. On
the other hand, the IMPORT.PRG estimates the SVAR
model, which is not a Near-SVAR model. In my program,
I’m trying to combine the MONTESUR.PRG and the
IMPORT.PRG to get the impulse responses in a Near-SVAR
model.

I think my problem is here-- I estimate 'betav' and
'beta', which respectively are the
covariance matrix of coefficients and the vector of
coefficients of the SUR model, out
side the IF-ELSE condition, and use these
outside-estimated coefficients to compute 'surdraw'.
Now If I want to compute "betav" in side the If-ELSE
condition, replacing "betav' by "%decomp(%xx), then
the resulting matrix doesn’t become the covariance
matrix of coefficients of the SUR model. Rather it
becomes a 5 x5 matirx, which comes from the number of
free parameters of the A matrix in the SVAR model.
Similarly, I need to be able to compute "beta", which
will mean the vector of coefficients of the SUR Model,
in the equation "compute surdraw=beta+betau".

Any ideas from anyone to fix it will be highly
appreciated. I’m using version 6.35.

Thanks,
Rokon

*Near-SVAR Model
/*Bwteen here and the next line of ****, I do the
following:
a) estimate the SUR model.
b) create the rectangle with dimension NREG X NEQN to
be used in %modelsetcoeffs
c) compute the covariance matrix of coefficients of
the Sur Model "betav=%decomp(%xx)", which is used
later to "compute betau=betav*ransur"
d) compute the vector of ceofficients of the Sur Model
"beta=%beta", which is used later to "compute
surdraw=beta+betau"
*/
compute lags=3
compute nvar=4
compute nstep=40
compute ncoef=nvar*lags+1
compute ndraws=1000

cal 1993 1 12
allocate 2006:12
open data svari.xls
data(format=xls,org=columns)

equation eq1 o
# constant o{1 to lags} r{1 to lags} inf{1 to lags}
ff{1 to lags}
equation eq2 r
# constant o{1 to lags} r{1 to lags} inf{1 to lags}
ff{1 to lags}
equation eq3 inf
# constant o{1 to lags} r{1 to lags} inf{1 to lags}
ff{1 to lags}
equation eq4 ff
# constant ff{1 to lags}

SUR(outsigma = vmat) 4
# eq1
# eq2
# eq3
# eq4

Group surmodel eq1 eq2 eq3 eq4

compute maxsize=0
do i=1,%modelsize(surmodel)
compute j=%eqnsize(%modeleqn(surmodel,i))
compute maxsize=%imax(maxsize,j)
end do j
dec rect fullcoeff(maxsize,%modelsize(surmodel))
dec rect[int] slots(2,%nreg)
compute fill=1
do i=1,%modelsize(surmodel)
do j=1,%eqnsize(%modeleqn(surmodel,i))
compute
slots(2,fill)=i,slots(1,fill)=j,fill=fill+1
end do j
end do i

dec vector ransur(%nreg) surdraw(%nreg)
declare vect[rect] responses(ndraws)
declare rect[series] impulses(nvar,nvar)
**
compute betav=%decomp(%xx)
compute beta=%beta
***********************************************************************************************
/* Between here and the next line of ******, I
estimate the free parameters of the A matrix of
the SVAR model, and put them in the vector called
"axbase"
*/
Nonlin(parmset=simszha) g12 g14 g21 g24 g31

compute nfree=5
dec frml[rect] afrml
frml afrml = ||1.0, g12, 0.0, g14|$
g21, 1.0, 0.0, g24|$
g31, 0.0, 1.0, 0.0|$
0.0, 0.0, 0.0, 1.0||

compute g12 = g14 = g21 = g24 = g31 = 0.0
cvmodel(parmset=simszha,pmethod=simplex,piters=10,iters=300,method=bfgs)
vmat afrml
compute delta=2.5
cvmodel(parmset=simszha,dfc=ncoef,pdf=delta,dmatrix=marginalized,method=bfgs)
vmat afrml
dec rect saxx
compute [vector] axbase=%parmspeek(simszha)
compute saxx=1.2*%decomp(%xx)
compute scladjust=%funcval
compute nu=12.0
****************************************************************************************************
declare vect weights(ndraws)
declare vect au(%nreg)
declare vect dbase d(nvar)
compute sumwt=0.0,sumwt2=0.0

/*Between here and the next line of *************,
a)Do a draw for the coefficients from a multivariate t
density and “poke” it back
into the parmset so that the AFRML can get the new
values
b)Compute (log kernels of)the true marginal posterior
density
and the importance function
*/

infobox(action=define,progress,lower=1,upper=ndraws)
'Monte Carlo Integration'
do draws = 1,ndraws
if %clock(draws,2)==1 {
compute grandom =nu/%rangamma(nu)
compute au =%ran(sqrt(grandom))
compute %parmspoke(simszha,axbase+saxx*au)
compute a =afrml(1)
compute dhat =a*vmat*tr(a)
compute ddiag =%xdiag(dhat)
compute
pdensity=.5*(%nobs-ncoef)*log(%det(a*tr(a)))-(.5*(%nobs-ncoef)+delta+1)*%sum(%log(ddiag))
compute
idensity=-((nu+nfree)/2.0)*log(nu+%dot(au,au))
******************************************************************************
*Compute the weight
compute wt
=exp(pdensity-scladjust-idensity-((nu+nfree)/2.0)*log(nu))

*Make a draw for D matrix
ewise d(i)
=(%nobs/2.0)*ddiag(i)/%rangamma(.5*(%nobs-ncoef)+delta+1)


*Generate the draw for a factor of sigma(the var-cov
matrix of the residuals) combining A and D
compute swish =inv(a)*%diag(%sqrt(d))


/*I think my problem is here-- I estimate 'betav' and
'beta', which respectively are the
covariance matrix of coefficients and the vector of
coefficients of the SUR model, out
side the IF-ELSE condition, and use these
outside-estimated coefficients to compute 'surdraw'.
Now If I want to compute "betav" in side the If-ELSE
condition,
replacing "betav' by "%decomp(%xx), then the resulting
matrix doesn't the covariance matrix of coefficients
of
the SUR model. Rather it becomes a 5 x5 matirx, which
comes from the number of free parameters
of the A matrix in the SVAR model. Similarly, I need
to be able to compute "beta", which will mean
the vector of coefficients of the SUR Model, in the
equation "compute surdraw=beta+betau".
*/

compute ransur=%ran(1.0)
compute betau=betav*ransur
compute surdraw=beta+betau
}
else {
compute surdraw=beta-betau
}
compute sumwt=sumwt+wt,sumwt2=sumwt2+wt**2
compute %matpoke(fullcoeff,slots,surdraw)

compute %modelsetcoeffs(surmodel,fullcoeff)

impulse(noprint,model=surmodel,decomp=swish,results=impulses)
nvar nstep
*
* Store the impulse responses
*
dim responses(draws)(nvar*nvar,nstep)
ewise
responses(draws)(i,j)=impulses((i-1)/nvar+1,%clock(i,nvar))(j)
compute weights(draws)=wt
infobox(current=draws)
end do draws
infobox(action=remove)

disp 'Effective sample size' sumwt**2/sumwt2
*
compute weights=weights/sumwt
*
dec vect[strings] xlabel(nvar) ylabel(nvar)
dec vect[integer] depvars
compute depvars=%modeldepvars(surmodel)
do i=1,nvar
compute ll=%l(depvars(i))
compute xlabel(i)=ll
compute ylabel(i)=ll
end do i

grparm(bold) hlabel 18 matrixlabels 14
grparm axislabel 24
spgraph(header='Impulse
responses',xpos=both,xlab=xlabel, $
ylab=ylabel,vlab='Responses
of',vfields=nvar,hfields=nvar)

dec vect[series] upper(nvar) lower(nvar) resp(nvar)
do i=1,nvar
compute minlower=maxupper=0.0
smpl 1 ndraws
do j=1,nvar
set lower(j) 1 nstep = 0.0
set upper(j) 1 nstep = 0.0
set resp(j) 1 nstep = 0.0
do k=1,nstep
set work 1 ndraws =
responses(t)((i-1)*nvar+j,k)
compute
frac=%wfractiles(work,weights,||.16,.84||)
compute lower(j)(k)=frac(1)
compute upper(j)(k)=frac(2)
compute resp(j)(k)=%avg(work)
end do k
compute
maxupper=%max(maxupper,%maxvalue(upper(j)))
compute
minlower=%min(minlower,%minvalue(lower(j)))
end do j
*
smpl 1 nstep
do j=1,nvar
graph(ticks,min=minlower,max=maxupper,number=0)
3 j i
# resp(j)
# upper(j) / 2
# lower(j) / 2
end do j
end do i
spgraph(done)