/**************     Local Polynomial Regression     ************************
Author: Cameron Rookley, 1999
        dh74673@goodnet.com
Submitted: FEb 1999
Code provided gratis without performance guarantees for public
non-commercial use.
***********************************************************/
/* the following program provides an example of
   local polynomial regression and various support
   routines.  For convenience all functions have been
   pasted onto the bottom of this test script.

   You may want to copy the individual functions into
   separate files with file names corresponding to
   the function names, with a .g extension on the end,
   and place the files in your \gauss\src directory for
   subsequent calling

   Cameron Rookley, 1999
   dh74673@goodnet.com  */

new;
cls;
library pgraph;
/* globals for lpoly function - see header of last
   function for details */

_nobmin=10; _kern=0; _progout=0; _der1=0; _der2=0; _heff=0;
ep=10^(-4);

/* x1 and x2 are the independent variables */
x1=rndn(1500,1);
x2=rndn(1500,1);

/* this is the true surface to be estimated */
truey=sqrt(abs(x1.*x2));

/* this is the dep. variable - true+noise */
/* you may want to play around with different levels of noise */

dep=truey+0.1.*rndn(1500,1);

/* this is the first step to setting up the mesh
   to estimate and graph the function at */
x1seq=linspace(minc(x1)-ep,maxc(x1)+ep,30);
x2seq=linspace(minc(x2)-ep,maxc(x2)+ep,30);

/* this creates a cartesian set of grid points in
   2 columns */
targmat=cartmatn(x1seq,x2seq);
trueytarg=sqrt(abs(targmat[.,1].*targmat[.,2]));

/* arbitrarily chosen bandwidths, play around with these
   to see the impact of over and under smoothing */
bw1=1; bw2=1;

p=0; /* order of the local polynomial to fit 1 typically works best
        play around with p=0 (observe how crappy it is at the borders)
        try p=2 with different bandwidth combinations */

der=0; /* function to be estimated der=0 estimates
          function, der=1 estimates 1st derivatives */


/* local polynomial regression call - smooths dep with respect to
   x1 and x2, estimating the surface at the points given in
   targmat.  Choosing the number of points in targmat to be
   much smaller than the data set greatly decreases computation time.
   bilinint is called afterwards to get full set of fitted values */

{yhat,varp2,delta}=lpoly(dep,x1~x2,bw1~bw2,p,targmat,der);

/* convert fitted values at target points into matrices for
   surface plotting and interpolating */

{x,y,z}=surfvec(yhat,targmat);
{x,y,z2}=surfvec(trueytarg,targmat);

/* plot the actual */
surface(x,y,z2);

/* plot the fit */
surface(x,y,z);

/* bilinear interpolate full dataset from target points */
{fullyhat}=bilinint(z,x,y,x1~x2);

/* compare fullyhat to true y, and dep */
corrx(fullyhat~truey~dep);



stop;


proc(1)=linspace(min,max,nobs);
retp(seqa(min,(max-min)/(nobs-1),nobs));
endp;

/* CARTMATN.G   Program to create a set of cartesian grid pairs from 2
unique
               vectors

   format: cmat=cartmatn(x1,x2) where x1 and x2 are vectors of unique
elements
                                    nx1    mx1
                               for the mesh. cmat is a (mxn)x2 matrix used
                               which contains every bivariate point in the
                               grid */
proc(1)=cartmatn(x1,x2);
local r1,r2,cmat,i;
r1=rows(x1); r2=rows(x2);
cmat=(x1[1].*ones(r2,1))~x2;
i=2;
do while i le r2;
   cmat=cmat|((x1[i].*ones(r2,1))~x2);
   i=i+1;
endo;
retp(cmat);
endp;

/* given data in the form of y   and    x which corresponds to
                           nx1         nx2
   a cartesian grid of target points, graph the surface for it

Format:  {x1,x2,z}=surfvec(y,x)

from this the command surface(x,y,z) can be issued as dimensions are
now properly set, where x1 corr's to unique elements of 1st col of x, and
x2 corr's to unique elements of 2nd col. of x */

proc(3)=surfvec(y,x);
local x1seq,x2seq,zhat,i,j,thisy;
x1seq=unique(x[.,1],1);
x2seq=unique(x[.,2],1);
if rows(x1seq)*rows(x2seq) ne rows(y);
   print "error with dimensions of inputs.  x may not be a stacked";
   print "cartesian grid";
endif;
zhat=zeros(rows(x2seq),rows(x1seq));
i=1;
do while i le rows(x2seq);
   j=1;
   do while j le rows(x1seq);
      thisy=selif(y,(x[.,1].==x1seq[j]).*(x[.,2].==x2seq[i]));
      zhat[i,j]=thisy;
      j=j+1;
   endo;
   i=i+1;
endo;
retp(x1seq',x2seq,zhat);
endp;

/* LPOLY.G  Procedure for univariate or multivariate nonparametric
            regression via local polynomial fitting.  See EBBS
            for data-driven bandwidth selection, with lpoly nested
            within.


Format: {yhat,varp2,delta}=lpoly(y,x,h,p,targmat,der)

  Outputs:   yhat nx1 vector of smoothed values for function of interest
                  given by der.  If der=0 yhat estimated, If der=1 an
                  estimate of the 1st der's is provided.  If der=2
                  an estimate of the second der's is provided.
                  This will either correspond to all
                  elements of x if targmat=0, or an mx1 vector
                  of smoothed values corresponding to non-zero targmat.

            varp2 Second component of estimated variance function.
                  (See note below)

Note: var(mhat'^k)=(sigma^2)*(k!)^2*e[k]'(inv(x'w*x)*(x'w.^2*x)*inv(x'w*x)
          varp2 is the component to the right of sigma^2.  Sigma^2 can
          be estimated from residuals or by smoothing squared resid's
          with correction. See ebbs.

            delta  diagonal of smoothing matrix suggested by Wand(1995)
used
                   for correcting bias, when estimating variance

   Inputs:   y is the dependent variable

             x is a nxk matrix of independent variables

             h is a common scaler bandwidth for each variable, or a kx1
               vector of bandwidths for each variable, or a nxk vector
               of local bandwidths

             p is a scaler indicating the degree of the polynomial to
               fit locally.  This must be less than or equal to 3.
               p = 1 or 2 recommended for estimating function

       targmat is an mxk matrix of target points to get smoothed values at
               or if targmat=0, then y will be smoothed at all data points.
               For large datasets, the use of targmat is highly
recommended,
               if feasible to drastically reduce computation time.

           der Indicates function of interest to estimate.  If der=0, then
               actual function estimated.  If der=1 then matrix of 1st
               derivitives estimates.  If der=2 then matrix of 2nd der's
               estimated, corr. to column 1, column2 etc. and then unique
               cross derivitives in the order
x1x2,x1x3..x1xn,x2x3,..x2xn...

        note:  Make sure targmat and x are of the same scale.


Globals: _kern     Kernel to use for weighting.  _kern=0 Gaussian
                                                 _kern=1 Epan.

         _nobmin   Specifies minimum number of observations required to
                   estimate local polynomial at each point.  If in a region
                   with sparse data, having this set to a certain number
will
                   cause the bandwidth to be locally cranked up to capture
                   more observations (Default=1+(k*p)+#of cross terms)

         _progout  If non-zero, the progress of the smooth in terms of
%'age
                   observations smoothed will be reported.
                   Use only if data set is large.
*/
proc(3)=lpoly(y,x,h,p,targmat,der);
local m,k,b1,b2,kind,yhat,y1,y2,n,i,cont,curh,kernarg,mkernarg;
local keepind,rk,xdiff,ki,w,ind,dep,krow,beta;
local ipi,invipi,varhat,delta,pc,sc;
local thisrow,varmat,wmat,smmat;
local kerncut;
/* #include lpoly.dec; */
n=rows(y);
if targmat==0;
  m=n;
  targmat=x;
else;
  m=rows(targmat);
endif;
k=cols(x);
if _nobmin==-99;
   if p==0; _nobmin=2; elseif p==1; _nobmin=k+1; elseif p==2;
      _nobmin=1+0.5*(k^2)+1.5*k; elseif p==3;
      _nobmin=1+1.5*k+1.5*(k^2);
   endif;
endif;
if rows(h)==1 and cols(h)==1;
     h=h.*ones(m,k);
elseif rows(h)==1 and cols(h)==k;
     h=h.*ones(m,k);
elseif rows(h)==k and cols(h)==1;
     h=h'.*ones(m,1);
elseif rows(h)==m and cols(h)==cols(targmat);
     /* we're okay */
elseif rows(h)==m and cols(h)==1;
     h=h.*ones(1,k);
else;
     print "rows of bandwidth vector inconsistent with cols of x";
     stop;
endif;
if der==0;
   yhat=zeros(m,1);
   if p ge 1;
      _der1=zeros(m,k);
   endif;
   if p ge 2;
      b2=seqa(k+2,1,0.5*(k^2+k));
      _der2=zeros(m,rows(b2));
   endif;
elseif der==1;
   if p < 1;
     print "resetting polynomial to 1, to estimate 1st derivitive";
     p=1;
   endif;
   b1=seqa(2,1,k);
   yhat=zeros(m,k);
elseif der==2;
   if p < 2;
      print "resetting polynomial to 2, to estimate 2nd derivitives";
      p=2;
   endif;
   b2=seqa(k+2,1,0.5*(k^2+k));
   yhat=zeros(m,rows(b2));
endif;
varhat=yhat;/* ie. zeros with same dimensions as yhat */
_heff=zeros(m,k);
delta=zeros(m,1);
if _kern==0;
   kerncut=4;
else;
   kerncut=1;
endif;
/* START OF MAIN LOOP */
/* USE CONDITIONAL BRANCHING AROUND MAIN LOOP FOR EFFICIENCY */
/* A LOT OF REPETITIVE CODE BELOW */
if _kern==0; /* copy this down */
if p==0 and der==0;
i=1;
do while i le m;
  curh=h[i,.]';
  cont=1;
  do while cont==1;
   kernarg=(x-targmat[i,.])./curh';
   mkernarg=maxc(abs(kernarg'));
   keepind=indexcat(mkernarg,-kerncut|kerncut);
   thisrow=maxindc(keepind.==i);
   rk=rows(keepind);
   if rk ge _nobmin;
      xdiff=x[keepind,.]-targmat[i,.];
      ki=prodc(pdfn((xdiff./curh')'));
      w=sqrt(ki);
      ind=ones(rk,1);
      ind=ind.*w;
      ipi=ind'ind;
      if abs(det(ipi)) > (10^(-20));
        invipi=inv(ipi);
        smmat=invipi*(ind.*w)';
        delta[i]=sumc(smmat[1,.]'.^2)-2.*smmat[1,thisrow];
        _heff[i,.]=curh';
        dep=y[keepind];
        beta=smmat*dep;
        varmat=invipi*(ind.*w)'(ind.*w)*invipi;
        yhat[i]=beta[1];
        varhat[i]=varmat[1,1];
        cont=0;
      else;
        curh=curh*1.1;
        cont=1;
      endif;
   else;
     /* pump up band-width to capture more observations */
     curh=curh*1.1;
     cont=1;
   endif;
  endo;
  i=i+1;
endo;
retp(yhat,varhat,delta);
elseif p==1 and der==0;
i=1;
do while i le m;
  /* i; */
  curh=h[i,.]';
  cont=1;
  do while cont==1;
   kernarg=(x-targmat[i,.])./curh';
   mkernarg=maxc(abs(kernarg'));
   keepind=indexcat(mkernarg,-kerncut|kerncut);
   thisrow=maxindc(keepind.==i);
   rk=rows(keepind);
   if rk ge _nobmin;
      xdiff=x[keepind,.]-targmat[i,.];
      ki=prodc(pdfn((xdiff./curh')'));
      w=sqrt(ki);
      ind=ones(rk,1);
      ind=ind~xdiff;
      ind=ind.*w;
      ipi=ind'ind;
      if abs(det(ipi)) > (10^(-20));
        invipi=inv(ipi);
        smmat=invipi*(ind.*w)';
        delta[i]=sumc(smmat[1,.]'.^2)-2.*smmat[1,thisrow];
        _heff[i,.]=curh';
        dep=y[keepind];
        beta=smmat*dep;
        varmat=invipi*(ind.*w)'(ind.*w)*invipi;
        yhat[i]=beta[1];
        varhat[i]=varmat[1,1];
        _der1[i,.]=beta[2:k+1]';
        cont=0;
      else;
        curh=curh*1.1;
        cont=1;
      endif;
   else;
     /* pump up band-width to capture more observations */
     curh=curh*1.1;
     cont=1;
   endif;
  endo;
  i=i+1;
endo;
retp(yhat,varhat,delta);
elseif p==2 and der==0;
i=1;
do while i le m;
  curh=h[i,.]';
  cont=1;
  do while cont==1;
   kernarg=(x-targmat[i,.])./curh';
   mkernarg=maxc(abs(kernarg'));
   keepind=indexcat(mkernarg,-kerncut|kerncut);
   thisrow=maxindc(keepind.==i);
   rk=rows(keepind);
   if rk ge _nobmin;
      xdiff=x[keepind,.]-targmat[i,.];
      ki=prodc(pdfn((xdiff./curh')'));
      w=sqrt(ki);
      ind=ones(rk,1);
      ind=ind~xdiff~(xdiff.^2);
      pc=1; /* prim. column */
      do while pc le k;
            sc=pc+1; /* sec. column */
            do while sc le k;
               ind=ind~xdiff[.,pc].*xdiff[.,sc];
               sc=sc+1;
            endo;
            pc=pc+1;
      endo;
      ind=ind.*w;
      ipi=ind'ind;
      if abs(det(ipi)) > (10^(-20));
        invipi=inv(ipi);
        smmat=invipi*(ind.*w)';
        delta[i]=sumc(smmat[1,.]'.^2)-2.*smmat[1,thisrow];
        _heff[i,.]=curh';
        dep=y[keepind];
        beta=smmat*dep;
        varmat=invipi*(ind.*w)'(ind.*w)*invipi;
        yhat[i]=beta[1];
        varhat[i]=varmat[1,1];
        _der1[i,.]=beta[2:k+1]';
        _der2[i,.]=beta[b2]';
        cont=0;
      else;
        curh=curh*1.1;
        cont=1;
      endif;
   else;
     /* pump up band-width to capture more observations */
     curh=curh*1.1;
     cont=1;
   endif;
  endo;
  i=i+1;
endo;
retp(yhat,varhat,delta);
elseif p==1 and der==1;
i=1;
do while i le m;
  curh=h[i,.]';
  cont=1;
  do while cont==1;
   kernarg=(x-targmat[i,.])./curh';
   mkernarg=maxc(abs(kernarg'));
   keepind=indexcat(mkernarg,-kerncut|kerncut);
   thisrow=maxindc(keepind.==i);
   rk=rows(keepind);
   if rk ge _nobmin;
      xdiff=x[keepind,.]-targmat[i,.];
      ki=prodc(pdfn((xdiff./curh')'));
      w=sqrt(ki);
      ind=ones(rk,1);
      ind=ind~xdiff;
      ind=ind.*w;
      ipi=ind'ind;
      if abs(det(ipi)) > (10^(-20));
        invipi=inv(ipi);
        smmat=invipi*(ind.*w)';
        delta[i]=sumc(smmat[1,.]'.^2)-2.*smmat[1,thisrow];
        _heff[i,.]=curh';
        dep=y[keepind];
        beta=smmat*dep;
        varmat=invipi*(ind.*w)'(ind.*w)*invipi;
        yhat[i,.]=beta[b1]';
        varhat[i,.]=(diag(varmat[b1,b1]))';
        cont=0;
      else;
        curh=curh*1.1;
        cont=1;
      endif;
   else;
     /* pump up band-width to capture more observations */
     curh=curh*1.1;
     cont=1;
   endif;
  endo;
  i=i+1;
endo;
retp(yhat,varhat,delta);
elseif p==2 and der==1;
i=1;
do while i le m;
  curh=h[i,.]';
  cont=1;
  do while cont==1;
   kernarg=(x-targmat[i,.])./curh';
   mkernarg=maxc(abs(kernarg'));
   keepind=indexcat(mkernarg,-kerncut|kerncut);
   thisrow=maxindc(keepind.==i);
   rk=rows(keepind);
   if rk ge _nobmin;
      xdiff=x[keepind,.]-targmat[i,.];
      ki=prodc(pdfn((xdiff./curh')'));
      w=sqrt(ki);
      ind=ones(rk,1);
      ind=ind~xdiff~(xdiff.^2);
      pc=1; /* prim. column */
      do while pc le k;
            sc=pc+1; /* sec. column */
            do while sc le k;
               ind=ind~xdiff[.,pc].*xdiff[.,sc];
               sc=sc+1;
            endo;
            pc=pc+1;
      endo;
      ind=ind.*w;
      ipi=ind'ind;
      if abs(det(ipi)) > (10^(-20));
        invipi=inv(ipi);
        smmat=invipi*(ind.*w)';
        delta[i]=sumc(smmat[1,.]'.^2)-2.*smmat[1,thisrow];
        _heff[i,.]=curh';
        dep=y[keepind];
        beta=smmat*dep;
        varmat=invipi*(ind.*w)'(ind.*w)*invipi;
        yhat[i,.]=beta[b1]';
        varhat[i,.]=(diag(varmat[b1,b1]))';
        cont=0;
      else;
        curh=curh*1.1;
        cont=1;
      endif;
   else;
     /* pump up band-width to capture more observations */
     curh=curh*1.1;
     cont=1;
   endif;
  endo;
  i=i+1;
endo;
retp(yhat,varhat,delta);
elseif p==2 and der==2;
i=1;
do while i le m;
  curh=h[i,.]';
  cont=1;
  do while cont==1;
   kernarg=(x-targmat[i,.])./curh';
   mkernarg=maxc(abs(kernarg'));
   keepind=indexcat(mkernarg,-kerncut|kerncut);
   thisrow=maxindc(keepind.==i);
   rk=rows(keepind);
   if rk ge _nobmin;
      xdiff=x[keepind,.]-targmat[i,.];
      ki=prodc(pdfn((xdiff./curh')'));
      w=sqrt(ki);
      ind=ones(rk,1);
      ind=ind~xdiff~(xdiff.^2);
      pc=1; /* prim. column */
      do while pc le k;
            sc=pc+1; /* sec. column */
            do while sc le k;
               ind=ind~xdiff[.,pc].*xdiff[.,sc];
               sc=sc+1;
            endo;
            pc=pc+1;
      endo;
      ind=ind.*w;
      ipi=ind'ind;
      if abs(det(ipi)) > (10^(-20));
        invipi=inv(ipi);
        smmat=invipi*(ind.*w)';
        delta[i]=sumc(smmat[1,.]'.^2)-2.*smmat[1,thisrow];
        _heff[i,.]=curh';
        dep=y[keepind];
        beta=smmat*dep;
        varmat=invipi*(ind.*w)'(ind.*w)*invipi;
        yhat[i,1:k]=2.*beta[b2[1:k]]';
        varhat[i,1:k]=4.*(diag(varmat[b2[1:k],b2[1:k]]))';
        if k>1;
           yhat[i,k+1:cols(yhat)]=beta[b2[k+1:rows(b2)]]';
           varhat[i,k+1:cols(yhat)]=
              (diag(varmat[b2[k+1:rows(b2)],b2[k+1:rows(b2)]]))';
        endif;
        cont=0;
      else;
        curh=curh*1.1;
        cont=1;
      endif;
   else;
     /* pump up band-width to capture more observations */
     curh=curh*1.1;
     cont=1;
   endif;
  endo;
  i=i+1;
endo;
retp(yhat,varhat,delta);
else;
/* old general loop */
i=1;
do while i le m;
  if _progout==1;
    print i/m;
  endif;
  if cols(h)>1;
     curh=h[i,.]';
  elseif cols(h)==1 and rows(h)==m;
     curh=h[i];
  else;
     curh=h;
  endif;
  cont=1;
  do while cont==1;
   kernarg=(x-targmat[i,.])./curh';
   mkernarg=maxc(abs(kernarg'));
   if _kern==0;
     keepind=indexcat(mkernarg,-4|4);
   else;
     keepind=indexcat(mkernarg,-1|1);
   endif;
   thisrow=maxindc(keepind.==i);
   rk=rows(keepind);
   if rk ge _nobmin;
      xdiff=x[keepind,.]-targmat[i,.];
      if _kern==0;
         ki=prodc(pdfn((xdiff./curh')'));
      else;
         ki=eppkern(xdiff./curh');
      endif;
      w=sqrt(ki);
      ind=ones(rk,1);
      if p ge 1;
         ind=ind~xdiff;
      endif;
      if p ge 2;
         ind=ind~xdiff.^2;
         pc=1; /* prim. column */
         do while pc le k;
            sc=pc+1; /* sec. column */
            do while sc le k;
               ind=ind~xdiff[.,pc].*xdiff[.,sc];
               sc=sc+1;
            endo;
            pc=pc+1;
         endo;
      endif;
      if p ge 3;
         ind=ind~(xdiff.^3);
         pc=1; /* prim. column */
         do while pc le k;
            sc=1; /* sec. column */
            do while sc le k;
               if sc ne pc;
                  ind=ind~(xdiff[.,pc].^2).*xdiff[.,sc];
               endif;
               sc=sc+1;
            endo;
            pc=pc+1;
          endo;
      endif;
      ind=ind.*w;
      ipi=ind'ind;
      if abs(det(ipi)) > (10^(-20));
        invipi=inv(ipi);
        smmat=invipi*(ind.*w)';
        delta[i]=sumc(smmat[1,.]'.^2)-2.*smmat[1,thisrow];
        _heff[i,.]=curh';
        dep=y[keepind];
        beta=smmat*dep;
        varmat=invipi*(ind.*w)'(ind.*w)*invipi;
        if der==0;
          yhat[i]=beta[1];
          varhat[i]=varmat[1,1];
          if p ge 1;
             _der1[i,.]=beta[2:k+1]';
          endif;
          if p ge 2;
             _der2[i,.]=beta[b2]';
          endif;
        elseif der==1;
          yhat[i,.]=beta[b1]';
          varhat[i,.]=(diag(varmat[b1,b1]))';
        elseif der==2;
          yhat[i,1:k]=2.*beta[b2[1:k]]';
          varhat[i,1:k]=4.*(diag(varmat[b2[1:k],b2[1:k]]))';
          if k>1;
           yhat[i,k+1:cols(yhat)]=beta[b2[k+1:rows(b2)]]';
           varhat[i,k+1:cols(yhat)]=
              (diag(varmat[b2[k+1:rows(b2)],b2[k+1:rows(b2)]]))';
          endif;
        endif;
        cont=0;
      else;
        curh=curh*1.1;
        cont=1;
      endif;
   else;
     /* pump up band-width to capture more observations */
     curh=curh*1.1;
     cont=1;
   endif;
  endo;
  i=i+1;
endo;
retp(yhat,varhat,delta);
endif;
else;
if p==0 and der==0;
i=1;
do while i le m;
  curh=h[i,.]';
  cont=1;
  do while cont==1;
   kernarg=(x-targmat[i,.])./curh';
   mkernarg=maxc(abs(kernarg'));
   keepind=indexcat(mkernarg,-kerncut|kerncut);
   thisrow=maxindc(keepind.==i);
   rk=rows(keepind);
   if rk ge _nobmin;
      xdiff=x[keepind,.]-targmat[i,.];
      ki=eppkern(xdiff./curh');
      w=sqrt(ki);
      ind=ones(rk,1);
      ind=ind.*w;
      ipi=ind'ind;
      if abs(det(ipi)) > (10^(-20));
        invipi=inv(ipi);
        smmat=invipi*(ind.*w)';
        delta[i]=sumc(smmat[1,.]'.^2)-2.*smmat[1,thisrow];
        _heff[i,.]=curh';
        dep=y[keepind];
        beta=smmat*dep;
        varmat=invipi*(ind.*w)'(ind.*w)*invipi;
        yhat[i]=beta[1];
        varhat[i]=varmat[1,1];
        cont=0;
      else;
        curh=curh*1.1;
        cont=1;
      endif;
   else;
     /* pump up band-width to capture more observations */
     curh=curh*1.1;
     cont=1;
   endif;
  endo;
  i=i+1;
endo;
retp(yhat,varhat,delta);
elseif p==1 and der==0;
i=1;
do while i le m;
  curh=h[i,.]';
  cont=1;
  do while cont==1;
   kernarg=(x-targmat[i,.])./curh';
   mkernarg=maxc(abs(kernarg'));
   keepind=indexcat(mkernarg,-kerncut|kerncut);
   thisrow=maxindc(keepind.==i);
   rk=rows(keepind);
   if rk ge _nobmin;
      xdiff=x[keepind,.]-targmat[i,.];
      ki=eppkern(xdiff./curh');
      w=sqrt(ki);
      ind=ones(rk,1);
      ind=ind~xdiff;
      ind=ind.*w;
      ipi=ind'ind;
      if abs(det(ipi)) > (10^(-20));
        invipi=inv(ipi);
        smmat=invipi*(ind.*w)';
        delta[i]=sumc(smmat[1,.]'.^2)-2.*smmat[1,thisrow];
        _heff[i,.]=curh';
        dep=y[keepind];
        beta=smmat*dep;
        varmat=invipi*(ind.*w)'(ind.*w)*invipi;
        yhat[i]=beta[1];
        varhat[i]=varmat[1,1];
        _der1[i,.]=beta[2:k+1]';
        cont=0;
      else;
        curh=curh*1.1;
        cont=1;
      endif;
   else;
     /* pump up band-width to capture more observations */
     curh=curh*1.1;
     cont=1;
   endif;
  endo;
  i=i+1;
endo;
retp(yhat,varhat,delta);
elseif p==2 and der==0;
i=1;
do while i le m;
  curh=h[i,.]';
  cont=1;
  do while cont==1;
   kernarg=(x-targmat[i,.])./curh';
   mkernarg=maxc(abs(kernarg'));
   keepind=indexcat(mkernarg,-kerncut|kerncut);
   thisrow=maxindc(keepind.==i);
   rk=rows(keepind);
   if rk ge _nobmin;
      xdiff=x[keepind,.]-targmat[i,.];
      ki=eppkern(xdiff./curh');
      w=sqrt(ki);
      ind=ones(rk,1);
      ind=ind~xdiff~(xdiff.^2);
      pc=1; /* prim. column */
      do while pc le k;
            sc=pc+1; /* sec. column */
            do while sc le k;
               ind=ind~xdiff[.,pc].*xdiff[.,sc];
               sc=sc+1;
            endo;
            pc=pc+1;
      endo;
      ind=ind.*w;
      ipi=ind'ind;
      if abs(det(ipi)) > (10^(-20));
        invipi=inv(ipi);
        smmat=invipi*(ind.*w)';
        delta[i]=sumc(smmat[1,.]'.^2)-2.*smmat[1,thisrow];
        _heff[i,.]=curh';
        dep=y[keepind];
        beta=smmat*dep;
        varmat=invipi*(ind.*w)'(ind.*w)*invipi;
        yhat[i]=beta[1];
        varhat[i]=varmat[1,1];
        _der1[i,.]=beta[2:k+1]';
        _der2[i,.]=beta[b2]';
        cont=0;
      else;
        curh=curh*1.1;
        cont=1;
      endif;
   else;
     /* pump up band-width to capture more observations */
     curh=curh*1.1;
     cont=1;
   endif;
  endo;
  i=i+1;
endo;
retp(yhat,varhat,delta);
elseif p==1 and der==1;
i=1;
do while i le m;
  curh=h[i,.]';
  cont=1;
  do while cont==1;
   kernarg=(x-targmat[i,.])./curh';
   mkernarg=maxc(abs(kernarg'));
   keepind=indexcat(mkernarg,-kerncut|kerncut);
   thisrow=maxindc(keepind.==i);
   rk=rows(keepind);
   if rk ge _nobmin;
      xdiff=x[keepind,.]-targmat[i,.];
      ki=eppkern(xdiff./curh');
      w=sqrt(ki);
      ind=ones(rk,1);
      ind=ind~xdiff;
      ind=ind.*w;
      ipi=ind'ind;
      if abs(det(ipi)) > (10^(-20));
        invipi=inv(ipi);
        smmat=invipi*(ind.*w)';
        delta[i]=sumc(smmat[1,.]'.^2)-2.*smmat[1,thisrow];
        _heff[i,.]=curh';
        dep=y[keepind];
        beta=smmat*dep;
        varmat=invipi*(ind.*w)'(ind.*w)*invipi;
        yhat[i,.]=beta[b1]';
        varhat[i,.]=(diag(varmat[b1,b1]))';
        cont=0;
      else;
        curh=curh*1.1;
        cont=1;
      endif;
   else;
     /* pump up band-width to capture more observations */
     curh=curh*1.1;
     cont=1;
   endif;
  endo;
  i=i+1;
endo;
retp(yhat,varhat,delta);
elseif p==2 and der==1;
i=1;
do while i le m;
  curh=h[i,.]';
  cont=1;
  do while cont==1;
   kernarg=(x-targmat[i,.])./curh';
   mkernarg=maxc(abs(kernarg'));
   keepind=indexcat(mkernarg,-kerncut|kerncut);
   thisrow=maxindc(keepind.==i);
   rk=rows(keepind);
   if rk ge _nobmin;
      xdiff=x[keepind,.]-targmat[i,.];
      ki=eppkern(xdiff./curh');
      w=sqrt(ki);
      ind=ones(rk,1);
      ind=ind~xdiff~(xdiff.^2);
      pc=1; /* prim. column */
      do while pc le k;
            sc=pc+1; /* sec. column */
            do while sc le k;
               ind=ind~xdiff[.,pc].*xdiff[.,sc];
               sc=sc+1;
            endo;
            pc=pc+1;
      endo;
      ind=ind.*w;
      ipi=ind'ind;
      if abs(det(ipi)) > (10^(-20));
        invipi=inv(ipi);
        smmat=invipi*(ind.*w)';
        delta[i]=sumc(smmat[1,.]'.^2)-2.*smmat[1,thisrow];
        _heff[i,.]=curh';
        dep=y[keepind];
        beta=smmat*dep;
        varmat=invipi*(ind.*w)'(ind.*w)*invipi;
        yhat[i,.]=beta[b1]';
        varhat[i,.]=(diag(varmat[b1,b1]))';
        cont=0;
      else;
        curh=curh*1.1;
        cont=1;
      endif;
   else;
     /* pump up band-width to capture more observations */
     curh=curh*1.1;
     cont=1;
   endif;
  endo;
  i=i+1;
endo;
retp(yhat,varhat,delta);
elseif p==2 and der==2;
i=1;
do while i le m;
  curh=h[i,.]';
  cont=1;
  do while cont==1;
   kernarg=(x-targmat[i,.])./curh';
   mkernarg=maxc(abs(kernarg'));
   keepind=indexcat(mkernarg,-kerncut|kerncut);
   thisrow=maxindc(keepind.==i);
   rk=rows(keepind);
   if rk ge _nobmin;
      xdiff=x[keepind,.]-targmat[i,.];
      ki=eppkern(xdiff./curh');
      w=sqrt(ki);
      ind=ones(rk,1);
      ind=ind~xdiff~(xdiff.^2);
      pc=1; /* prim. column */
      do while pc le k;
            sc=pc+1; /* sec. column */
            do while sc le k;
               ind=ind~xdiff[.,pc].*xdiff[.,sc];
               sc=sc+1;
            endo;
            pc=pc+1;
      endo;
      ind=ind.*w;
      ipi=ind'ind;
      if abs(det(ipi)) > (10^(-20));
        invipi=inv(ipi);
        smmat=invipi*(ind.*w)';
        delta[i]=sumc(smmat[1,.]'.^2)-2.*smmat[1,thisrow];
        _heff[i,.]=curh';
        dep=y[keepind];
        beta=smmat*dep;
        varmat=invipi*(ind.*w)'(ind.*w)*invipi;
        yhat[i,1:k]=2.*beta[b2[1:k]]';
        varhat[i,1:k]=4.*(diag(varmat[b2[1:k],b2[1:k]]))';
        if k>1;
           yhat[i,k+1:cols(yhat)]=beta[b2[k+1:rows(b2)]]';
           varhat[i,k+1:cols(yhat)]=
              (diag(varmat[b2[k+1:rows(b2)],b2[k+1:rows(b2)]]))';
        endif;
        cont=0;
      else;
        curh=curh*1.1;
        cont=1;
      endif;
   else;
     /* pump up band-width to capture more observations */
     curh=curh*1.1;
     cont=1;
   endif;
  endo;
  i=i+1;
endo;
retp(yhat,varhat,delta);
else;
/* old general loop */
i=1;
do while i le m;
  if _progout==1;
    print i/m;
  endif;
  if cols(h)>1;
     curh=h[i,.]';
  elseif cols(h)==1 and rows(h)==m;
     curh=h[i];
  else;
     curh=h;
  endif;
  cont=1;
  do while cont==1;
   kernarg=(x-targmat[i,.])./curh';
   mkernarg=maxc(abs(kernarg'));
   if _kern==0;
     keepind=indexcat(mkernarg,-4|4);
   else;
     keepind=indexcat(mkernarg,-1|1);
   endif;
   thisrow=maxindc(keepind.==i);
   rk=rows(keepind);
   if rk ge _nobmin;
      xdiff=x[keepind,.]-targmat[i,.];
      if _kern==0;
         ki=prodc(pdfn((xdiff./curh')'));
      else;
         ki=eppkern(xdiff./curh');
      endif;
      w=sqrt(ki);
      ind=ones(rk,1);
      if p ge 1;
         ind=ind~xdiff;
      endif;
      if p ge 2;
         ind=ind~xdiff.^2;
         pc=1; /* prim. column */
         do while pc le k;
            sc=pc+1; /* sec. column */
            do while sc le k;
               ind=ind~xdiff[.,pc].*xdiff[.,sc];
               sc=sc+1;
            endo;
            pc=pc+1;
         endo;
      endif;
      if p ge 3;
         ind=ind~(xdiff.^3);
         pc=1; /* prim. column */
         do while pc le k;
            sc=1; /* sec. column */
            do while sc le k;
               if sc ne pc;
                  ind=ind~(xdiff[.,pc].^2).*xdiff[.,sc];
               endif;
               sc=sc+1;
            endo;
            pc=pc+1;
          endo;
      endif;
      ind=ind.*w;
      ipi=ind'ind;
      if abs(det(ipi)) > (10^(-20));
        invipi=inv(ipi);
        smmat=invipi*(ind.*w)';
        delta[i]=sumc(smmat[1,.]'.^2)-2.*smmat[1,thisrow];
        _heff[i,.]=curh';
        dep=y[keepind];
        beta=smmat*dep;
        varmat=invipi*(ind.*w)'(ind.*w)*invipi;
        if der==0;
          yhat[i]=beta[1];
          varhat[i]=varmat[1,1];
          if p ge 1;
             _der1[i,.]=beta[2:k+1]';
          endif;
          if p ge 2;
             _der2[i,.]=beta[b2]';
          endif;
        elseif der==1;
          yhat[i,.]=beta[b1]';
          varhat[i,.]=(diag(varmat[b1,b1]))';
        elseif der==2;
          yhat[i,1:k]=2.*beta[b2[1:k]]';
          varhat[i,1:k]=4.*(diag(varmat[b2[1:k],b2[1:k]]))';
          if k>1;
           yhat[i,k+1:cols(yhat)]=beta[b2[k+1:rows(b2)]]';
           varhat[i,k+1:cols(yhat)]=
              (diag(varmat[b2[k+1:rows(b2)],b2[k+1:rows(b2)]]))';
          endif;
        endif;
        cont=0;
      else;
        curh=curh*1.1;
        cont=1;
      endif;
   else;
     /* pump up band-width to capture more observations */
     curh=curh*1.1;
     cont=1;
   endif;
  endo;
  i=i+1;
endo;
retp(yhat,varhat,delta);
endif;
endif;
endp;

proc(1)=eppkern(x);
  retp( (prodc(abs(x').le 1)).*(0.75.^(cols(x))).*prodc(1-(x'.^2)));
endp;


/* bilinint.g

   Bivariate linear interpolation

   Format:

   {yhat}=bilinint( y ,  x1seq , x2seq,  x)

    (Lx1)          MXN    1XN     MX1   Lx2


  Outputs: yhat  is an estimate of the function at the target points x


  Inputs:

      y        Matrix of function values given by the cartesian grid
               corresponding to the column values x1seq, and row values
               x2seq (use surfvec if you need to set up y and sequence
               matrices in this format, starting with a vector for the
               dep var. and matrix pair for the indep. vars


        x      Matrix of target points at which to evaluate the function

*/
proc(1)=bilinint(y,x1seq,x2seq,x);
local x1lrow,x1l,x1r,x2brow,x2t,x2b,wl,wr,wt,wb,wtmat,n,ymat,i,ybox;
local totarea,yvec;
if rows(x2seq)==1 and cols(x2seq)>1;
   x2seq=x2seq';
endif;
if rows(x1seq)>1 and cols(x1seq)==1;
   x1seq=x1seq';
endif;
if cols(y) ne cols(x1seq);
      print "error with columns of y";
      stop;
endif;
if rows(y) ne rows(x2seq);
   print "error with rows of y";
   stop;
endif;
if minc(x[.,1])<minc(x1seq) or maxc(x[.,1])>maxc(x1seq);
   print "1st col. of x obs. falls outside of grid";
   stop;
endif;
if minc(x[.,2])<minc(x2seq') or maxc(x[.,2])>maxc(x2seq');
   print "1st col. of x obs. falls outside of grid";
   stop;
endif;
if maxc(x[.,1])==maxc(x1seq');
   x1seq=x1seq~(x1seq[cols(x1seq)]+1);
   y=y~y[.,cols(y)];
endif;
if maxc(x[.,2])==maxc(x2seq);
   x2seq=x2seq|(x2seq[rows(x2seq)]+1);
   y=y|y[rows(y),.];
endif;
x1lrow=minindc(((x[.,1]-x1seq)+(10^200).*(x[.,1].<x1seq))');
x1l=x1seq[x1lrow]';
x1r=x1seq[x1lrow+1]';
x2brow=minindc(((x[.,2]-x2seq')+(10^200).*(x[.,2].<x2seq'))');
x2b=x2seq[x2brow];
x2t=x2seq[x2brow+1];
wl=x[.,1]-x1l;
wr=x1r-x[.,1];
wt=x2t-x[.,2];
wb=x[.,2]-x2b;
/* tl tr bl br */
wtmat=(wr.*wb)~(wl.*wb)~(wr.*wt)~(wl.*wt);
n=rows(x);
ymat=zeros(n,4);
i=1;
do while i le n;
   ybox=y[(x2brow[i]+1)|x2brow[i],
                    x1lrow[i]~(x1lrow[i]+1)];
   ymat[i,.]=(ybox[1,.]~ybox[2,.]);
   i=i+1;
endo;
totarea=(x1r-x1l).*(x2t-x2b);
yvec=(totarea.>0).*sumc((ymat.*wtmat)')./(totarea);
retp(yvec);
endp;



