/* MVKS.G: multivariate kernel smoother 
           John Rust, University of Maryland, October, 2006 

	   y=f(x)+e

	   given n observations in vector ydata and xdata,

	   the multivariate Nadaraya-Watson kernel smoother is

	   ^
	   f(x) = sum


	   
	  
	   inputs:
	   
	   x    vector of points to compute nonparametric regression at
	   data vector of data points used to compute kernel density
	  
	   implicit global inputs

	   kerneltyp: type of kernel density to use in the density estimation

	              epanechnikov
		      gaussian
	   
	   */

   proc (1)=mvks(x,ydata,xdata);  
  
       local kti,kerneltyp,mdkernel,bwrule,rx,rd,nvars,npe,i,j,h,s,tmp,tmp1;

       kti=missrv(typecv("kerneltyp"),0);

       if (kti == 0);
          kerneltyp="epanechnikov";
	  kti=varput(kerneltyp,"kerneltyp");
       else;
          kerneltyp=varget("kerneltyp");
       endif;

       kti=missrv(typecv("mdkernel"),0);

       if (kti == 0);
          mdkernel="unidist";
	  kti=varput(mdkernel,"unidist");
       else;
          mdkernel=varget("mdkernel");
       endif;

       kti=missrv(typecv("bwrule"),0);

       if (kti == 0);
          bwrule="optimal";
	  kti=varput(bwrule,"optimal");
       else;
          bwrule=varget("bwrule");
       endif;

       rx=rows(x);
       rd=rows(xdata);
       nvars=cols(xdata);

       if (nvars /= cols(x));
         "Error: colums of 1st and 3rd arguments of mvks must be the same";
	 retp(miss(0,0));
       endif;

       npe=zeros(rx,1);
       
       h=1.06/((rd)^(.2));  /* Default Bandwidth parameter, Silverman rule
                               for 1 dimensional data */
       s=stdc(xdata);

       if (kerneltyp $== "gaussian");

          if (nvars > 1);
             if (bwrule $== "optimal");
	       h=((4/(nvars+2))^(1/(nvars+4)))/(rd^(1/(nvars+4)));
	       h=h*ones(nvars,1);
h=h*2.5;
	     endif;
	  endif;

          i=1; do until i > rx;
             tmp=zeros(rd,1);
             if (bwrule $== "adaptive");
	       h=adaptbw(x[i,.],xdata);
             endif;
	     j=1; do until j > nvars;
              tmp=tmp+(-(x[i,j]-xdata[.,j])^2)/(2*h[j]*h[j]*s[j]*s[j]);
             j=j+1; endo;
	     npe[i]=sumc(ydata.*exp(tmp))/sumc(exp(tmp));
	  i=i+1; endo;
       
       elseif (kerneltyp $== "tricube");
	  
          if (nvars > 1);
             if (bwrule $== "optimal");
	       h=3.12/(rd^(1/(nvars+4)));
	       h=h*ones(nvars,1);
	     endif;
	  endif;
	  
          i=1; do until i > rx;
           tmp=ones(rd,1);
           if (bwrule $== "adaptive");
	       h=adaptbw(x[i,.],xdata);
           endif;
	   if (mdkernel $== "product");
             j=1; do until j > nvars;
                tmp=tmp.*(abs((x[i,j]-xdata[.,j])/(s[j]*h[j])).<=1).*((ones(rd,1)-abs((x[i,j]-xdata[.,j])/(s[j]*h[j]))^3)^3);
	     j=j+1; endo;
	     tmp1=sumc(tmp);
           else;
	     j=1; do until j > nvars;
                tmp=tmp+((x[i,j]-xdata[.,j])/(s[j]*h[j]))^2;
	     j=j+1; endo;
             tmp=((ones(rd,1)-tmp)^3).*(tmp.<1);
	     tmp1=sumc(tmp);
	   endif;
	   
	   if (tmp1 <= 0);
	      npe[i]=0;
           else;
              npe[i]=sumc(ydata.*tmp)/tmp1;
           endif;
         
	  i=i+1; endo;
    
        else;
	  
          if (nvars > 1);
             if (bwrule $== "optimal");
	       h=2.4/(rd^(1/(nvars+4)));
	       h=h*ones(nvars,1);
	     endif;
	  endif;

          i=1; do until i > rx;
           tmp=ones(rd,1);
           if (bwrule $== "adaptive");
	       h=adaptbw(x[i,.],xdata);
           endif;
	   if (mdkernel $== "product");
             j=1; do until j > nvars;
                tmp=tmp.*((((x[i,j]-xdata[.,j])/(s[j]*h[j]))^2).<=5).*(ones(rd,1)-.2*(((x[i,j]-xdata[.,j])/(s[j]*h[j]))^2));
	     j=j+1; endo;
	     tmp1=sumc(tmp);
           else;
	     j=1; do until j > nvars;
                tmp=tmp+((x[i,j]-xdata[.,j])/(s[j]*h[j]))^2;
	     j=j+1; endo;
             tmp=(ones(rd,1)-tmp).*(tmp.<1);
	     tmp1=sumc(tmp);
	   endif;
	   
	   if (tmp1 <= 0);
	        npe[i]=0;
           else;
                npe[i]=sumc(ydata.*tmp)/tmp1;
           endif;
	  i=i+1; endo;

       endif;

       retp(npe);


   endp;