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.