/* UNIVARIATE NON-PARAMETRIC REGRESSION PACKAGE                       */
/* BY Cameron Rookley Nov 1996                                        */
/* see bnonparm.g for BIVARIATE VERSION                               */
/*     mnonparm.g for MULTIVARIATE VERSION

   Note: The following files must all be included in your src directory
         unonparm.g, unonparm.dec, sgrid1.g, sgrid1.dec.

         See unonparm.txt for notes and hints on using this package.

   Purpose: to estimate univariate functions using non-parametric
            regression using either a Nadayara Watson weighting scheme.

   Format:  {mx,h,h1,h2,der1,der2}=unonparm(y,x,h0,h0d,h0d2)

   where:   y  is an NX1 vector, representing the dependent variable
               for the non-parametric regression.

            x  is an NX1 vector, representing the independent variable for
               the nonparametric regression.

            h0 is an Lx1 vector of bandwidths to use in estimating m(x),
               If this is set to zero a grid search will be carried
               out using cross-validation.

    h0d & h0d2 are Mx1 and PX1 vectors of bandwidths used to estimate
               1st & 2nd derivities of y with respect to x, if _der1 ne 0
               and/or _der2 ne 0 (see globals below).  If these are set
               to zero then bandwidths will be chosen using either cross
               validation or the Rice Criteria (Default).

    OUTPUTS:

          mx is an NxL matrix of fitted values for the dependent variable
             for each bandwidth supplied by h0, or is an Nx1 vector of
             fitted values corresponding to the optimal bandwidth if h0=0.

   h,h1 & h2 are scalers representing optimal bandwidths from
             cross validation or Rice criteria (for derivitives).

        der1 is a NxM matrix of estimated derivitives of y with respect
             to X. These derivities will only be estimated if_ders is
             non-zero.

        der2 is a NxP matrix of estimated second derivitives of y with
             respect to X for each given bandwidth. These derivities will
             only be estimated if_der2 is non-zero.

  Global Variables:

      _fast     Scaler specifying whether to use the fast-fourier
                transform.  This will massively decrease computation
                time if set to 1, but at the cost of some accuracy if
                independent variable is not an equidistant grid.
                Default=0 (do not use).
                Note if _fast=1, _kern (see below) will be
                reset to 1 for computational efficiency.  Currently
                derivitive estimation is not set up using FFT, so ensure
                _ders and _ders2 are set to 0 if data set is large.

      _resist   Scaler specifying whether to smooth while being resistant
                to outliers. (Default=0 no resistant smoothing).  If set
                to 1, then residuals exceeding +- _rescut standard deviatons
                of the 1st step residuals will be treated as outliers.  If
                _resist=2, then _rescut gives the cut off point in absolute
                terms for outliers.

      _rescut   Cutoff for residuals when using resistant smoothing.
                (Default=2 when _resist=1)

      _kern     scaler specifying which type of Kernel to use for m(x)
                (Default=1)
                0: Epan./Bartlett (1963) K(u)=0.75(1-u^2)I(|u|<=1)
                1: Gaussian kernel K(u)=exp(-u^2/2)/(sqrt(2*pi))
                3: 3'rd order kernel
                   K(u)=(15/32)*(3-3u-10u^2+10u^3+7u^4-7u^5)
                4: 4'th order kernel
                   K(u)=(3/sqrt(8*pi))*(1-(z^2)/3)*exp(-(z^2)/2)
                5: 5'th order kernel
                   K(u)=(25/256)*(15-105u^2+189u^4-99u^6)
                       -0.2099(-5u+35u^3-63u^5+33u^7)


      _gmin     scalers which determine the nature of the grid
      _gmax     search for bandwidth selection.  Using sgrid1.g
      _nints    (see description) we will search from _gmin to _gmax
      _niters   in _nints intervals to find optimal h. _niters finer
                Defaults: _gmin=0.00001, _gmax=2*stddevc(x), _nints=20,
                          _niters=20

      _bwcm1    scalers representing Bandwidth Choice Methods
      _bwcm2    for estimating 1st and 2nd derivities.
                If element=0 (default) then a Factor Method based
                on Rice criteria will be used.
                If element=1 a separate grid search will be
                carried out, using cross-validation and difference
                quotients.  Note: if non-unique elements of indep. var
                exist then these values will be reset to default.

      _ders     scaler specifying whether to calculate derivitives of
                y with respect to x. (Default=0)

      _ders2    same as _ders but for 2nd der's of y wrt x. (Default=0)

               _gmind1, _gmaxd1, _nintsd1, _niterd1
               _gmind2, _gmaxd2, _nintsd2, _niterd2
                Same as above but for finding bandwidths for 1st and 2nd
                derivitive estimation. (Same defaults)

      _kernd1   Same as _kern but used for 1st derivitive
                estimation. (Default=3)

      _kernd2   Same as _kernd1 but used for 2nd derivitive estimation
                (Default=4)

      _wtfnh    percentage of observations to include for weighting
                cross-validation functions (default=90%).
                This is used to prevent boundary behaviour from
                affecting band-width choice.

     _gridout   Specifies whether intermediate output from cross-validation
                grid searches is to be printed. (Default=1, yes print)

*/
proc(6)=unonparm(dep,ind,h0,h0d,h0d2);
local h,cval,gradcv,rccv,mhat,j,hnew;
local gridint;
local hd,gradcvd,rccvd,hdmat,_globhd,hd1,hd2;
local der1,der2,lab;
local cvflag,cd2flag,xys;
local bwcflag;
local c;
#include unonparm.dec;
_gridlb=1/(10^4);
/* _smflag=1; */
_ep=1/(10^250);
bwcflag=0;
_n=rows(ind);
if _fast==1;
   _kern=1;
endif;
_aseq=seqa(1,1,_n);
_d=cols(ind);
if _d>1;
   print "multivariate data encountered, please use nonparm function";
   stop;
endif;
if _ders ne 0;
    if _kernd1 ne 3 and _kernd1 ne 5;
      print "Invalid kernel choice for _kernd1";
      print "Please choose a kernel of order 3 or order 5";
      stop;
    endif;
    if _bwcm1 ne 0 and _bwcm1 ne 1;
        print "Invalid value for _bwcm1, please set all values to 0 or 1";
        stop;
    endif;
    if _gmaxd1==0;
      _gmaxd1=_gmax;
    endif;
    if _gmind1==0;
      _gmind1=_gmin;
    endif;
endif;
if _ders2 ne 0;
   if _kernd2 ne 4 and _kernd2 ne 6;
      print "Invalid kernel choice for der. element of _kernd2";
      print "Must be a kernel of order 4 or order 6 for 2nd der's";
      stop;
   endif;
   if _bwcm2 ne 0 and _bwcm2 ne 1;
      print "Invalid value for _bwcm2 please set all elements to 0 or 1";
      stop;
   endif;
   if _gmaxd2==0;
     _gmaxd2=_gmax;
   endif;
   if _gmind2==0;
      _gmind2=_gmin;
   endif;
endif;

_x=ind; clear ind;
_y=dep; clear dep;

xys=sortc(_x~_y,1);
_x=xys[.,1];
_y=xys[.,2];

/* find cutoffs for percentage of observations to include in
   cv calculations to limit boundary effects */
{_lcut,_ucut}=getcuts(_x);

if _gmax==0;
  _gmax=stdc(_x);
endif;
if _gmin==0;
  _gmin=0.1.*stdc(_x);
endif;
lab="h      "~"ch(h)   ";
/* choose optimal h for m(x) where necessary */
if h0==0; /* need to optimize */
    if _gridout==1;
     print "finding optimal bandwidth for m(x):";
     print;
     print $lab;
     print;
    endif;
    {hnew,cval}=sgrid1(_gmin,_gmax,_nints,_niters,&cv);
    if hnew==-99 and cval==-99; /* ie grid search fails */
       bwcflag=1; /* level fail flag */
    endif;
    h=hnew;
else;
  h=h0;
endif;

if _ders ne 0;
    if h0d==0;
     if _gridout==1 and h0d==0;
       print "Finding optimal bandwidth for der. of y wrt x";
       print;
       print $lab;
       print;
     endif;
     if rows(uniquem(_x,1))<_n;
        if _gridout==1;
         print "Non-unique elements encountered.  Switching to Rice Criteria";
        endif;
        _bwcm1=0;
     endif;
     if _bwcm1==1;
          {hnew,cval}=sgrid1(_gmind1,_gmaxd1,_nintsd1,_niterd1,&cvd);
          hd=hnew;
          if hnew==-99 and cval==-99; /* ie grid search fails */
             if _gridout==1;
               print "Grid search failed, switching to Rice criteria";
             endif;
            _bwcm1==0; /* use rice criteria */
          endif;
     endif;
     if _bwcm1==0; /* use rice criteria see Muller pg. 103 */
          if _kernd1==3;
            hd=0.7083*h;
          elseif _kernd1==5;
            hd=0.7396*h;
          endif;
     endif;
    else;
        hd=h0d;
    endif;
else;
  hd=0;
endif;

/* find optimal bandwidth for second der's */
if _ders2 ne 0;
   if h0d2==0;
     if _gridout==1;
       print "Finding optimal bandwidth for SECOND der. of y wrt x";
       print;
       print $lab;
       print;
     endif;
     if rows(uniquem(_x,1))<_n;
        _bwcm2=0;
     endif;
     if _bwcm2==1;
        {hnew,cval}=sgrid1(_gmind2,_gmaxd2,_nintsd2,_niterd2,&cvd2);
        hd2=hnew;
          if hnew==-99 and cval==-99; /* ie grid search fails */
            _bwcm2==0; /* use rice criteria */
          endif;
     endif;
     if _bwcm2==0; /* use rice criteria */
          if _kernd2==4;
            hd2=0.8919*h;
          elseif _kernd2==6;
            hd2=0.9507*h;
          else;
            print "Invalid kernel choice for der. element of _kernd2";
            print "Must be a kernel of order 4 or order 6";
            stop;
          endif;
     endif;
   else;
       hd2=h0d2;
   endif;
else;
  hd2=0;
endif;

/* estimate mhat,m'hat, and m''hat */
if bwcflag==0;
  mhat=zeros(rows(_y),rows(h));
  c=1;
  do while c le rows(h);
    mhat[.,c]=nadwat(_y,_x,h[c]);
    c=c+1;
  endo;
else;
  mhat=-99;
endif;
if _ders>0;
  der1=zeros(rows(_y),rows(hd));
  c=1;
  do while c le rows(hd);
    der1[.,c]=nadwatd(_y,_x,hd[c]);
    c=c+1;
  endo;
else;
  der1=zeros(_n,1);
endif;
if _ders2>0;
   der2=zeros(rows(_y),rows(hd2));
   c=1;
   do while c le rows(hd2);
     der2[.,c]=nadwatd2(_y,_x,hd2[c]);
     c=c+1;
   endo;
else;
   der2=zeros(_n,1);
endif;


retp(mhat,h,hd,hd2,der1,der2);
endp;
/* ******************* */
/* END OF MAIN ROUTINE */
/* ******************* */


/* ************************************************ */
/* **********PROCEDURES *************************** */
/* ************************************************ */

/* Procedure to find cut-off points of x, to correspond with
   _wtfnh */
proc(2)=getcuts(x);
local tx,lrow,urow,lb,ub;
lrow=maxc(1|(round(((1-_wtfnh)/2)*_n)));
urow=minc(_n|(_n-round(((1-_wtfnh)/2)*_n)));
tx=sortc(x,1);
lb=tx[lrow];
ub=tx[urow];
retp(lb,ub);
endp;

/* ***************************************************  */
/* Procedure to find Cross-Validation function for m(x) */
proc cv(h);
local thish,i,ki,m,cv,numv,denv;
local index,keepdum,xt;
local fftdata,fftkern;
numv=zeros(_n,1);
denv=zeros(_n,1);
i=1;
do while i le _n;
    if _kern==0 or _kern==3 or _kern==5 or _kern==6;
       index=indexcat((_x-_x[i])./h,-1.01|1);
    elseif _kern==1 or _kern==2 or _kern==4;
       index=indexcat((_x-_x[i])./h,-5.8|5.8);
    else;
       print "invalid value for _kern encountered";
       stop;
    endif;
   if missrv(index,-99) ne -99;
      ki=getki(_x[index],_x[i],h,_kern);
      numv[index]=numv[index]+(index.ne i).*ki.*_y[i];
      denv[index]=denv[index]+(index.ne i).*ki;
   endif;
   i=i+1;
 endo;
 m=(denv.>0).*(numv./(denv+_ep));
 cv=sumc(((_y-m).^2).*(_x.<_ucut).*(_x.>_lcut));
retp(-cv);  /* maximize the negative of cv */
endp;


/* ****************************************************  */
/* Procedure to find Cross-Validation function for m'(x) */
proc cvd(h);
local kdy,sk,kly,skd,keepdum,i,m,cv;
local dyvec,mdyvec,keepdum2,index;
local kl,kd;

kdy=zeros(_n,1);
sk=zeros(_n,1);
kly=zeros(_n,1);
skd=zeros(_n,1);

i=1;
do while i le _n;
    if _kernd1==0 or _kernd1==3 or _kernd1==5;
       index=indexcat((_x-_x[i])./h,-1.01|1);
    elseif _kernd1==1 or _kernd1==2 or _kernd1==4;
       index=indexcat((_x-_x[i])./h,-5.8|5.8);
    else;
       print "Invalid value for _kernd1 encountered.";
       stop;
    endif;
   if missrv(index,-99) ne -99;
      keepdum2=(index.ne i).*(index.ne(i-1));
      kl=getki(_x[index],_x[i],h,_kernd1);
      kd=(((1/h).*(index.ne i))+(index.==i)).*
                                getkd(_x[index],_x[i],h,_kernd1);
      skd[index]=skd[index]+keepdum2.*kd;
      sk[index]=sk[index]+keepdum2.*kl;
      kdy[index]=kdy[index]+keepdum2.*kd.*_y[i];
      kly[index]=kly[index]+keepdum2.*kl.*_y[i];
   endif;
   i=i+1;
endo;
m=(sk.>0).*((kdy./(sk+_ep))-(skd.*kly)./((sk.^2)+_ep));
dyvec=(_y[2:_n]-_y[1:_n-1])./(_x[2:_n]-_x[1:_n-1]);
mdyvec=submat(dyvec~m[1:_n-1]~_x[1:_n-1],seqa(1,2,floor(_n/2)),0);
cv=(1/(floor(_n/2))).*sumc( ((mdyvec[.,1]-mdyvec[.,2]).^2).*
                   (mdyvec[.,3].<_ucut).*(mdyvec[.,3].>_lcut));
retp(-cv);  /* maximize the negative of cv */
endp;

/* ****************************************************  */
/* Procedure to find Cross-Validation function for m''(x) */
proc cvd2(h);
local keepdum,i,ki,m,cv;
local dyvec,mdyvec,keepdum2,index,xt;
local dy2vec;
local kl,kd,kd2,kdy,skl,kly,skd,kd2y,skd2;

kdy=zeros(_n,1);
skl=zeros(_n,1);
kly=zeros(_n,1);
skd=zeros(_n,1);
kd2y=zeros(_n,1);
skd2=zeros(_n,1);
i=1;
do while i le _n;
   if _kernd2==4;
       index=indexcat((_x-_x[i])./h,-5.8|5.8);
   elseif _kernd2==6;
       index=indexcat((_x-_x[i])./h,-1|1);
   endif;
   if missrv(index,-99) ne -99;
     keepdum2=(index.ne i).*(index.ne(i-1)).*(index.ne(i+1));
     kl=getki(_x[index],_x[i],h,_kernd2);
     kd=(((1/h).*(index.ne i))+(index.==i)).*getkd(_x[index],_x[i],h,_kernd2);
     kd2=(((1/(h^2)).*(index.ne i))+(index.==i)).*getkd2(_x[index],_x[i],h,_kernd2);
     kd2y[index]=kd2y[index]+keepdum2.*kd2.*_y[i];
     skd2[index]=skd2[index]+keepdum2.*kd2;
     kdy[index]=kdy[index]+keepdum2.*kd.*_y[i];
     skl[index]=skl[index]+keepdum2.*kl;
     kly[index]=kly[index]+keepdum2.*kl.*_y[i];
     skd[index]=skd[index]+keepdum2.*kd;
   endif;
   i=i+1;
endo;
m=(skl.>0).*((kd2y./(skl+_ep))-(kdy.*skd)./((skl.^2)+_ep)-
  (kly.*skd2)./((skl.^2)+_ep)-(kdy.*skd)./((skl.^2)+_ep)+
   2.*(kly.*(skd.^2))./((skl.^3)+_ep));
dyvec=(_y[2:_n]-_y[1:_n-1])./(_x[2:_n]-_x[1:_n-1]);
dy2vec=(dyvec[2:_n-1]-dyvec[1:_n-2])./(0.5.*(_x[3:_n]-_x[1:_n-2]));
mdyvec=submat(dy2vec~m[2:_n-1]~_x[2:_n-1],seqa(1,3,floor(_n/3)),0);
cv=(1/(floor(_n/3))).*sumc(((mdyvec[.,1]-mdyvec[.,2]).^2).*
               (mdyvec[.,3].<_ucut[_j]).*(mdyvec[.,3].>_lcut[_j]));
retp(-cv); /* maximize the negative of cv */
endp;

/* ***** Procedure to get kernel of x-x[i,.] ***** */
proc getki(x,xi,h,kern);
local ki;
if kern==0;
     ki=(bart((x-xi)./h));
elseif kern==1 or kern==2;
     ki=(gaussk((x-xi)./h));
elseif kern==3;
     ki=(threeok((x-xi)./h));
elseif kern==4;
     ki=(fourok((x-xi)./h));
elseif kern==5;
     ki=(fiveok((x-xi)./h));
elseif kern==6;
     ki=(sixok((x-xi)./h));
else;
     print "Unknown kernel code encountered. Check kern.";
     stop;
endif;
retp(ki);
endp;

/* ***** Procedure to get Derivitive of kernel of x-x[i,.] ***** */
proc getkd(x,xi,h,kern);
local ki;
if kern==3;
     ki=(threeok1((x-xi)./h));
elseif kern==4;
     ki=(fourok1((x-xi)./h));
elseif kern==5;
     ki=(fiveok1((x-xi)./h));
elseif kern==6;
     ki=(sixok1((x-xi)./h));
else;
     print "Unknown or too low, kernel code encountered. Check _kernd1.";
     stop;
endif;
retp(ki);
endp;

/* ***** Procedure to get 2nd Derivitive of kernel of x-x[i,.] ***** */
proc getkd2(xloc,xi,h,kern);
local ki;
    if kern==4;
       ki=(fourok2((xloc-xi)./h));
    elseif kern==6;
       ki=(sixok2((xloc-xi)./h));
    else;
       print "Invalid kernel code encountered for second derivitive.";
       print "Check _kernd2 to ensure 4th or 6th order kernel chosen.";
       stop;
    endif;
retp(ki);
endp;

/* procedure to calculate m according to
   Nadaraya Watson weighting scheme */
proc nadwat(y,x,h);
local numv,denv,m,i,ki,j;
local keepdum,index;
local optseq,n,xrange,b,a,c,fftdata;
local res,rescut,fftres,psinum,pdarg,psiden;
if _fast ne 1;
 numv=zeros(_n,1);
 denv=zeros(_n,1);
 i=1;
 do while i le _n;
   if _kern==0 or _kern==6;
       index=indexcat((_x-_x[i])./h,-1.01|1);
   elseif _kern==1 or _kern==2 or _kern==4;
       index=indexcat((_x-_x[i])./h,-6|6);
   elseif _kern==3 or _kern==5;
       print "must choose an even ordered kernal for est. of m(x)";
       stop;
   endif;
   if missrv(index,-99) ne -99;
      ki=getki(_x[index],_x[i],h,_kern);
      numv[index]=numv[index]+ki.*_y[i];
      denv[index]=denv[index]+ki;
   endif;
   i=i+1;
 endo;
 m=(denv.>0).*(numv./(denv+_ep));
else;
 y=_y|zeros(8/h,1); n=rows(y);
 optseq=2.^(seqa(1,1,25));
 i=1;
 do while i le 25;
   if n<optseq[i];
      y=y|zeros(optseq[i]-n,1);
      i=26;
   else;
      i=i+1;
   endif;
 endo;
 n=rows(y);
 xrange=maxc(x)-minc(x);
 xrange=xrange*(n/_n);
 b=-2*((pi*h)/(xrange))^2;
 a=1;
 c=exp(b.*(seqa(1,1,n/2).^2));
 fftdata=fft(y.*a);
 numv=real(ffti(fftdata.*(1|c[1:(n/2)-1]|rev(c))));
 denv=real(ffti(fft(y.ne 0).*(1|c[1:(n/2)-1]|rev(c))));
 m=numv./(real(denv)+_ep);
 m=real(m[1:_n]);
endif;
if _resist ne 0;
  res=_y-m;
  if _resist==1;
    rescut=sqrt(vcx(res)).*_rescut;
  elseif _resist==2;
    rescut=_rescut;
  else;
    print "resist must be an integer from 0 to 2";
    stop;
  endif;
  res=maxc(res'|-rescut.*ones(1,_n));
  res=minc(res'|rescut.*ones(1,_n));
  if _fast ne 1;
    psinum=zeros(_n,1);
    psiden=zeros(_n,1);
    i=1;
    do while i le _n;
       if _kern==0 or _kern==6;
         keepdum=(((abs((_x-_x[i])./h')).le 1));
       elseif _kern==1 or _kern==2 or _kern==4;
         keepdum=(((abs((_x-_x[i])./h')).le 5.8));
       elseif _kern==3 or _kern==5;
         print "must choose an even ordered kernal for est. of m(x)";
         stop;
       endif;
       xt=selif(_x~_aseq,keepdum);
       if cols(xt)>1;
        index=xt[.,2];
        xt=xt[.,1];
        ki=getki(xt,_x[i],h,_kern);
        psinum[index]=psinum[index]+ki.*res[i];
        psiden[index]=psiden[index]+ki.*(res[i] > -rescut).*(res[i] < rescut);
       endif;
       i=i+1;
    endo;
    m=m+(psiden.>0).*(psinum./(psiden+_ep));
  else;
    res=res|zeros(n-_n,1);
    fftres=fft(res.*a);
    psinum=real(ffti(fftres.*(1|c[1:(n/2)-1]|rev(c))));
    pdarg=(fft((a.*(res[1:_n].< rescut).*(res[1:_n].> -rescut))|zeros(n-_n,1))).*
                     (1|c[1:(n/2)-1]|rev(c));

    psiden=real(ffti(pdarg));
    m=m+(psiden[1:_n].>0).*(psinum[1:_n]./(psiden[1:_n]+_ep));
  endif;
endif;
retp(m);
endp;

/* Procedure to find m'(x) using Nadayra Watson weights */
proc nadwatd(y,x,h);
local kdy,sk,kly,skd,keepdum,i,md;
local index;
local kl,kd;
kdy=zeros(_n,1);
sk=zeros(_n,1);
kly=zeros(_n,1);
skd=zeros(_n,1);
i=1;
do while i le _n;
    if _kernd1==3 or _kernd1==5;
       index=indexcat((_x-_x[i])./h,-1.01|1);
    else;
       print "Must use an odd-ordered kernal for _kernd1.";
       stop;
   endif;
   if missrv(index,-99) ne -99;
     kl=getki(_x[index],_x[i],h,_kernd1);
     kd=(((1/h).*(index.ne i))+(index.==i)).*
                     getkd(_x[index],_x[i],h,_kernd1);
     skd[index]=skd[index]+kd;
     sk[index]=sk[index]+kl;
     kdy[index]=kdy[index]+kd.*_y[i];
     kly[index]=kly[index]+kl.*_y[i];
   endif;
   i=i+1;
endo;
md=(sk.>0).*((kdy./(sk+_ep))-(skd.*kly)./((sk.^2)+_ep));
retp(md);
endp;

/* Procedure to find m''(x) using Nadayra Watson weights */
proc nadwatd2(y,x,h);
local keepdum,i,ki,md2;
local index,xt;
local kl,kd,kd2,kdy,skl,kly,skd,kd2y,skd2;
kdy=zeros(_n,1);
skl=zeros(_n,1);
kly=zeros(_n,1);
skd=zeros(_n,1);
kd2y=zeros(_n,1);
skd2=zeros(_n,1);
i=1;
do while i le _n;
   if _kernd2==4;
     index=indexcat((_x-_x[i])./h,-5.8|5.8);
   elseif _kernd2==6;
     index=indexcat((_x-_x[i])./h,-1.01|1);
   else;
     print "Must choose either a 4th or 6th order kernel for m''(x)";
     stop;
   endif;
   if missrv(index,-99) ne -99;
    kl=getki(_x[index],_x[i],h,_kernd2);
    kd=(((1/h).*(index.ne i))+(index.==i)).*
                         getkd(_x[index],_x[i],h,_kernd2);
    kd2=(((1/(h^2)).*(index.ne i))+(index.==i)).*
                         getkd2(_x[index],_x[i],h,_kernd2);
    kd2y[index]=kd2y[index]+kd2.*_y[i];
    skd2[index]=skd2[index]+kd2;
    kdy[index]=kdy[index]+kd.*_y[i];
    skl[index]=skl[index]+kl;
    kly[index]=kly[index]+kl.*_y[i];
    skd[index]=skd[index]+kd;
   endif;
   i=i+1;
endo;
md2=(skl.>0).*((kd2y./(skl+_ep))-(kdy.*skd)./((skl.^2)+_ep)-
  (kly.*skd2)./((skl.^2)+_ep)-(kdy.*skd)./((skl.^2)+_ep)+
   2.*(kly.*(skd.^2))./((skl.^3)+_ep));
retp(md2);
endp;

/* **** Bartlett kernel procedure */
proc bart(u);
retp((abs(u).le1).*(0.75.*(1-u.^2)));
endp;

/* ****  Gaussian kernel procedure */
proc gaussk(u);
retp(exp((-(u.^2))./2)./sqrt(2*pi));
endp;

/* **** Third Order kernel         */
proc threeok(u);
  retp((15/32).*(3-3.*u-10.*(u.^2)+10.*(u.^3)+7.*(u.^4)-7.*(u.^5)));
endp;

/* ****  Fourth-order kernel       */
proc fourok(u);
retp((3/sqrt(8*pi)).*(1-(u.^2)/3).*exp(-(u.^2)./2) );
endp;

/* *** Fifth order kernel          */
proc fiveok(u);
retp( (35/256).*(15-105.*(u.^2)+189.*(u.^4)-99.*(u.^6))-
      0.2099.*(-5.*u+35.*(u.^3)-63.*(u.^5)+33.*(u.^7)));
endp;

/* *** Sixth order Kernal          */
proc sixok(u);
retp( (315/64).*(77.*(u.^6)-135.*(u.^4)+63.*(u.^2)-5));
endp;

proc sixok1(u);
retp( (315/64).*(462.*(u.^5)-540.*(u.^3)+126.*u));
endp;

proc sixok2(u);
retp( (315/64).*(2310.*(u.^4)-1620.*(u.^2)+126));
endp;


/* *** First derivitive of Third Order kernel */
proc threeok1(u);
retp( (15/32).*(-3-20.*u+30.*(u.^2)+28.*(u.^3)-35.*(u.^4)));
endp;

/* *** First derivitive of Fifth Order kernel */
proc fiveok1(u);
retp( (35/256).*(-210.*u+756.*(u.^3)-594.*(u.^5))-
      0.2099.*(-5+105.*(u.^2)-315.*(u.^4)+231.*(u.^6)));
endp;


/* 1st derivitive of 4'th Order kernel */
proc fourok1(u);
retp(-u.*fourok(u)-((2.*u)./(sqrt(8*pi))).*exp((-u.^2)./2));
endp;

/* 2nd derivitive of 4'th order kernel */
proc fourok2(u);
retp( (-fourok(u)-(2/sqrt(8*pi)).*exp(-(u.^2)./2))-
       u.*(fourok1(u)- ((2.*u)./sqrt(8*pi)).*exp(-(u.^2)./2)));
endp;

proc fgaussk(u);
 retp(exp(-(u.^2)/.2));
endp;
