/* KERNEL DENSITY LIBRARY ver 1.0
Author: Ruud H. Koning rhkoning@xs4all.nl.
Date: 13 July 1996
Provided without guarantee for public noncommercial use.
*/

/* kernel.src version 1.0 13 july 1996
Source file with functions of the kernel library:
    ukernel     procedure for univariate kernel estimation
    mkernel     procedure for multivariate kernel estimation
    bandw1      procedure for automatic bandwidth selection
    nw          procedure for nonparametric regression
    viewuknl    procedere for plotting the estimated density

*/

/* comment the following line out if you don't have the pgraph library
installed. You will not be able to use the procedure view_ukernel */
#include pgraph.ext

/* procedure to calculate a univariate kernel density estimate and its
derivative
usage:      {f, d, h} = ukernel(x, z, h, w, &kf);
input:      x:      T vector where density is to be estimated
            z:      n vector with observed data points
            h:      scalar bandwidth, if h<=0 bandwidth is determined by
                    procedure bandw1(z)
            w:      n vector with weights
            &kf:    pointer to weighting function
output:     f:      T vector with estimated density
            d:      T vector with estimated derivative
            h:      scaler bandwidth
*/

proc (3)=ukernel(x, z,h,w, &kf);
 local arg, i, n, f, d, k, kff, kfd, pkff,argl,argh,kffl,kffh,fl,fh,
  hl,hh;
 local kf:proc;

 /* error checking here */
 if (cols(x)>1);
  errorlog "ukernel.g: x has too many columns";
  retp(-1,-1,-1 );
 endif;
 if (cols(z)>1);
  errorlog "ukernel.g: z has too many columns";
  retp( -1,-1,-1 );
 endif;

 /* initialization */
 i = 1;
 n = rows(x);
 k = cols(x);

 /* determine bandwidth */
 if (h<=0);
  h=bandw1(z);
 endif;

 f=zeros(n,1);
 d=zeros(n,1);
 do while (i<=n);
  arg = (x[i]-z)/h;
  {kff, kfd} = kf(arg);
  f[i] = meanc(kff.*w)/h;
  d[i] = meanc(kfd.*w)/(h^2);
  i=i+1;
 endo;
 retp( f,d,h );
endp;


/* procedure to calculate a multivariate kernel density estimate and its
derivative
usage:      {f, d, h} = mkernel(x, z, h, w, &kf);
input:      x:      T x k matrix where density is to be estimated
            z:      n x k matrix with observed data points
            h:      scalar, bandwidth parameter (same bandwidth for all
                    components), if h<=0 it is determined by eq 4.14 of
                    Silverman
            w:      n vector with weights
            &kf:    pointer to weighting function
output:     f:      T vector with estimated density
            d:      T x k matrix with estimated derivatives
            h:      scalar, bandwidth used
*/

proc (3)=mkernel(x, z, h, w, &kf);
 local arg,i,n,f,d,k,kff,kfd,pkff,p,deneq0;
 local kf:proc;

 i = 1;
 n = rows(x);
 k = cols(x);
 if (h<=0);
   p=1/(k+4);
   h=(4/(k+2))^p*n^(-p);
 endif;
 f = zeros(n,1);
 d = zeros(n,k);
 do while (i<=n);
  arg = (x[i,.]-z)/h;
  {kff, kfd} = kf(arg);
  pkff = prodc(kff');
  f[i] = meanc(pkff.*w)/(h^k);
  deneq0=kff .==0;
  kff=kff+deneq0;
  d[i,.] = meanc((kfd.*pkff.*w)./kff)'/(h^(k+1));
  i = i+1;
 endo;
 retp( f,d,h );
endp;


/* computes Nadaraya-Watson kernel regression estimator

usage:  {g,d,h} = nw(x,z,y,h,&kf);
input:  x:          T-vector where regression function is evaluated
        z:          n-vector with data on independent variable
        y:          n-vector with data on dependent variable
        h:          scalar, bandwidth, if h<=0 the bandwidth is determined
                    automatically using procedure bandw1
        &kf:        pointer to univariate kernel density estimator
output: g:          T-vector with values of regression evaluated in x
        d:          T-vector with estimated derivative in x
        h:          scalar bandwidth used
*/

proc (3)=nw(x, z, y,h,&kf);
    local num, numder, denom, denomder, g, s, n, k, d,v;
    local kf: proc;

 n = rows(x);
 k = cols(x);
 g = zeros(n,1);
 d = zeros(n,1);

 if (k ne 1);
  print "error in nw.g: too many columns in x";
  retp( g, d );
 endif;
 if (cols(z) ne 1);
  print "error in nw.g: too many columns in z";
  retp( g, d );
 endif;

 {num,numder,h} = ukernel(x,z,h,y,&kf);
 {denom,denomder,h} = ukernel(x,z,h,1,&kf);
 s = (denom.==0);
 g = num./(denom+s);
 d = numder./(denom+s) - num.*denomder./(denom^2+s);
 retp( (1-s).*g, (1-s).*d,h );
endp;

/* bandw1
procedure to calculate the optimal bandwidth in kernel estimation of a density.
The optimal bandwidth is calculated according to eq. 3.31 of Silverman (1986)

usage:      h=bandw1(y);
input:      y:      n-vector whose density will be estimated;
output:     h:      scalar, optimal bandwidth choice;

*/

proc bandw1(y);
 local s, ys,n, a, iqr, qi1, qi3;

 if (cols(y)>1);
  errorlog "input error in bandw1.g: too many columns";
  retp( -1 );
 endif;
 s=sqrt(vcx(y));
 n=rows(y);
 ys=sortc(y,1);
 qi1=round(0.25*n);
 qi3=round(0.75*n);
 iqr=ys[qi3]-ys[qi1];
 retp( 0.9*minc(s|(iqr/1.34))/n^0.2 );
endp;

/* view_ukernel
procedure to plot kernel density estimate. Library pgraph must be activated.

usage: call view_ukernel(x,f,h);
input:      x:  n-vector with data points
            f:  n x k matrix (k=1 or k=3) with estimated densities for
                different bandwidths
            h:  k-vector (k=1 or k=3) with bandwidths
output:     none
globals:    all globals of the pgraph library

*/

proc (0)=view_ukernel(x,f,h);
 local data,k,xlow,xhigh,flow,fhigh,xlegend,ylegend;

 /* error checking */
 if (cols(f) ne rows(h));
  errorlog "error in viewkrnl.g: rows(f) unequal to cols(h)";
  retp();
 endif;
 /* set global variables pgraph */
 _pdate=0;

 k=cols(f);
 xlow=minc(x);
 xhigh=maxc(x);
 flow=minc(minc(f));
 fhigh=maxc(maxc(f));
 xlegend=xlow;
 ylegend=flow+0.8*(fhigh-flow);
 _plegctl=1|4|xlegend|ylegend;
 if (k==1);
  _plegstr="h=" $+ ftos(h,"%*.*lf",5,3);
  else;
  _plegstr="h=" $+ ftos(h[1],"%*.*lf",5,3) $+ "\0h="
   $+ftos(h[2],"%*.*lf",5,3)$+"\0h="$+ ftos(h[3],"%*.*lf",5,3);
 endif;
 data=sortc(x~f,1);
 xy(data[.,1],data[.,2:cols(data)]);
endp;



