Bivariate density estimation

[jonasdovern]

Dear RATS users,

Is there a way to estimate the bivariate density function for two series of data in RATS? (What the density-function is doing in the univariate case for one series of data.)

Best, Jonas

[TomDoan]

This includes the estimation of a bivariate density function using two outputs from a Gibbs sampling routine. Operations using xgrid, ygrid and fgrid are the ones that are needed for the estimation of the density function. The choice of grids and the "h" bandwidth will depend upon the application.

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")

Another example is in the simsmodel.prg (or simsmodel.rpf) program in the Sims and Zha, Econometrica 1999 replication, which is included with RATS and is also at:

http://www.estima.com/forum/viewtopic.php?f=4&t=335

[jonasdovern]

Tom, thanks a lot for the example.

[jonasdovern]

I made an easy to use procedure from the code that Tom shared. Might be helpful ...

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

The procedure estimates the joint density of the two series using a Gaussian kernel.