@=====================================================================================@
@                        GAUSS Program for Density Estimation                         @
@=====================================================================================@

new;
library pgraph;
load data[462,3]=tbsm.dat;              @ read in data @
y=data[.,1];
ret=y[2:rows(y)]-y[1:rows(y)-1];        @ interest rate returns @
nn=rows(ret);                           @ number of observations @
yst=(ret-meanc(ret))/stdc(ret);         @ standardized returns @
@xy(0,yst); @                             @ plotting the data @

{px1,py1}=edf(yst);                     @ empirical distribution function @
@xy(px1,py1);@

@ Density Estimates @
optband=1.06*(nn^(-1/5));               @ optimal bandwidth for iid data @
{px2,py2}=dens(yst,3,optband);          @ density estimation with Uniform kernel @
@xy(px2,py2);@
{px3,py3}=dens(yst,1,optband);          @ density estimation with Gaussian kernel @
{px4,py4}=dens(yst,2,optband);          @ density estimation with Epanechnikov kernel @ 
@xy(px3,py3~py4~py2);@
py5=pdfn(px3);                          @ standard normal density @
xy(px3,py2~py3~py4~py5);  


@ Supplemental Procesures @

@ Empirical Distribution Function @
proc(2)=edf(y);
local strt,endd,kern,pts,px,py,i;
 strt=minc(y); endd=maxc(y);            @ min and max of the series @
 pts=100;                               @ number of gridpoints @
 px=seqa(strt,(endd-strt)/(pts-1),pts); @ grid of equally spaced points @
 py=zeros(pts,1);
 i=1; do while i<=pts;
  py[i]=meanc(y.<=px[i]);              @ empirical distribution function at point i @
 i=i+1; endo;
retp(px,py);
endp;

@ Kernel Density Estimate @
proc(2)=dens(y,krn,h);
local strt,endd,kern,pts,px,py,std,i,t;
 strt=minc(y); endd=maxc(y);            @ min and max of the series @
 pts=100;                               @ number of gridpoints @
 px=seqa(strt,(endd-strt)/(pts-1),pts); @ grid of equally spaced points @
 py=zeros(pts,1);
 i=1; do while i<=pts;
  if krn==1; t=nkernel(px[i,1],y,h)./h;
   elseif krn==2; t=ekernel(px[i,1],y,h)./h;
   elseif krn==3; t=ukernel(px[i,1],y,h)./h;
  endif;
  py[i,1]=sumc(t)/rows(y);              @ density kernel estimate at point i @
 i=i+1; endo;
retp(px,py);
endp;

@ Gaussian kernel @
proc nkernel(x,xi,h);
local u,kern;
 u=(x-xi)./h;
 kern=(1./sqrt(2*pi))*exp(-(1/2).*(u^2));
retp(kern);
endp;

@ Epanechnikov kernel @
proc ekernel(x,xi,h);
local u,kern;
 u=(x-xi)./h;
 kern=(0.75*(1-u^2)).*(abs(u).<=1);
retp(kern);
endp;

@ Uniform kernel @
proc ukernel(x,xi,h);
local u,kern;
 u=(x-xi)./h;
 kern=0.5*(abs(u).<=1);
retp(kern);
endp;