A Forum for RATS Software and Econometrics Discussions
https://www.estima.com/forum/
Code: Select all
*
* Replication file for Kadiyala and Karlsson, "Numerical Methods for
* Estimation and Inference in Bayesian-VAR models", Journal of Applied
* Econometrics, 1997, vol 12, no 2, pp 99-132.
*
* Simulation of Normal-diffuse prior from page 106
*
dec symm sigmahat(2,2)
compute sigmahat=1.8/28.0*%identity(2)
compute [vect] betahat=||1.0,1.0||
compute [vect] betabar=||8.0,9.0||
compute [symm] betasig=%identity(2)
compute [symm] zz=||.91||
*
* Bandwidth for bivariate Gaussian window for kernel density estimator
*
compute h=.5
*
* Set up grid in the plane
*
dec vect xgrid(45)
dec vect ygrid(45)
dec rect fgrid(45,45)
ewise xgrid(i)=.2*i
ewise ygrid(i)=.2*i
compute fgrid=%const(0.0)
*
compute ndraws=10000
compute nburn =100
*
* Compute draws using Gibbs sampler
*
dec vect u(2)
set b1 1 ndraws = 0.0
set b2 1 ndraws = 0.0
compute [symm] sigmabase = 28.0*sigmahat
compute [symm] sigmascale = sigmabase
*
do draws=-nburn+1,ndraws
*
* Draw sigma (capital psi in KK's notation)
*
compute sigma=%ranwisharti(%decomp(inv(sigmascale)),27)
*
* Draw beta (KK's gamma) given sigma
*
compute [symm] betaxxinv = %kroneker(inv(sigma),zz)+inv(betasig)
compute [vect] betamean = inv(betaxxinv)*(%kroneker(inv(sigma),zz)*betahat+inv(betasig)*betabar)
compute [vect] betadraw = betamean+%decomp(inv(betaxxinv))*(u=%ran(1.0))
if draws>0 {
compute b1(draws)=betadraw(1)
compute b2(draws)=betadraw(2)
*
* Update the density across the grid.
*
ewise fgrid(i,j)=fgrid(i,j)+exp(-.5/h**2*((betadraw(1)-xgrid(i))**2+(betadraw(2)-ygrid(j))**2))
}
*
* Compute the scale parameter for the next draw for sigma
*
compute sigmascale=sigmabase+.91*%innerxx(||betadraw(1)-betahat(1),betadraw(2)-betahat(2)||)
end do draws
*
* If we really wanted f to be a true density function (integrating to
* one), we would need to divide fgrid by ndraws * (2*pi) * h**2. Since
* we're just interested in the shape, we can skip that step.
*
gcontour(x=xgrid,y=ygrid,f=fgrid,header="Density of Draws from Gibbs Sampling")
*
* Compute the marginal density using equation (9).
*
compute ggrid=%zeros(45,45)
ewise ggrid(i,j)=exp(%logdensity(betasig,||xgrid(i)-betabar(1),ygrid(j)-betabar(2)||))*$
%det(28.0*sigmahat+.91*%innerxx(||xgrid(i)-betahat(1),ygrid(j)-betahat(2)||))**(-28.0/2.0)
gcontour(x=xgrid,y=ygrid,f=ggrid,header="Direct Calculation of Density")Code: Select all
procedure %%bvdensity xser yser start end
type series xser yser
type integer start end
option real bandwidth .5
option real xmin
option real xmax
option real ymin
option real ymax
option integer ygridsize 100
option integer xgridsize 100
option integer number 30
local integer i j k nobs lstart lend
local real lxmin lxmax lymin lxmin
local vec xgrid ygrid
local rec fgrid
*
Written by Jonas Dovern, Kiel Economics Research & Forecasting GmbH & Co. KG, 01/2012.
*
*** Setting up grids
inquire(regressorlist) lstart>>start lend>>end
# xser yser
disp lstart lend
comp nobs = lend-lstart+1
dim xgrid(xgridsize) ygrid(ygridsize)
comp fgrid = %zeros(xgridsize,ygridsize)
statistics(noprint,fract) xser
if %defined(xmin)
comp lxmin = xmin
else
comp lxmin = %minimum
if %defined(xmax)
comp lxmax = xmax
else
comp lxmax = %maximum
statistics(noprint,fract) yser
if %defined(ymin)
comp lymin = ymin
else
comp lymin = %minimum
if %defined(ymax)
comp lymax = ymax
else
comp lymax = %maximum
ewise xgrid(i) = lxmin+(i-1)*(lxmax-lxmin)/(xgridsize-1)
ewise ygrid(i) = lymin+(i-1)*(lymax-lymin)/(ygridsize-1)
*
* Calculating density
infobox(action=define,progress,lower=1,upper=nobs) "Estimating density"
do k=lstart,lend
infobox(action=modify,current=k)
ewise fgrid(i,j) = fgrid(i,j)+exp(-.5/bandwidth^2.*((xser(k)-xgrid(i))^2+(yser(k)-ygrid(j))^2))
end do k
ewise fgrid(i,j) = fgrid(i,j)/(nobs*2*%pi*bandwidth^2)
infobox(action=remove)
* Plotting density
gcontour(x=xgrid,y=ygrid,f=fgrid,header="Contour plot of density estimate",number=number)
end procedure