/* LLR.GPR: program to do multidimensional interpolation/extrapolation
            via local linear regression
            John Rust, University of Maryland, October, 2006

   proc (1)=llr(pv,v,x);

   Inputs:
  
   pv: point in R^k at which a function v is to be approximated
   v:  nx1 vector of values of a function v evaluated at x
   x:  nxk matrix of points at which v is evaluated

   Output:

   v(x):  interpolated value at point x

   Note: n should be substantially larger than k   

   Parameters:

   wr:   set to 1 for weighted regression, where
         observations further from p are downweighted
         set to 0 for unweighted regression.

   h:    bandwidth  */

   proc (1)=llr(pv,v,x);

      local kti,kerneltyp,mdkernel,bwrule,result,h,i,j,p,rp,n,d,dp,ny,nx,s,wts,bv;

      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;

      n=rows(x);
      d=cols(x);
      s=stdc(x);

      rp=rows(pv);
      h=1.06/(rows(v)^(.2)); 
      result=zeros(rp,1);
      
      if (kerneltyp $== "gaussian");

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

          i=1; do until i > rp;
             dp=pv[i,.]-x;
	     wts=zeros(n,1);
             if (bwrule $== "adaptive");
	       h=adaptbw(pv[i,.],x);
             endif;
	     j=1; do until j > d;
              wts=wts+(-(pv[i,j]-x[.,j])^2)/(2*h[j]*h[j]*s[j]*s[j]);
             j=j+1; endo;
	     
             wts=exp(-(wts^2)/2);
             nx=(ones(n,1)~dp).*sqrt(wts); 
             ny=v.*sqrt(wts);
             bv=ny/nx;
             result[i]=bv[1]; 
	     
	  i=i+1; endo;
       
      elseif (kerneltyp $== "tricube");
	  
          if (d > 1);
             if (bwrule $== "optimal");
	       h=3.12/(n^(1/(d+4)));
	       h=h*ones(d,1);
	     endif;
	  endif;
	  
          i=1; do until i > rp;
	  
           wts=ones(n,1);
           dp=pv[i,.]-x;
           if (bwrule $== "adaptive");
	       h=adaptbw(pv[i,.],x);
           endif;
	   if (mdkernel $== "product");
             j=1; do until j > d;
                wts=wts.*(abs((pv[i,j]-x[.,j])/(s[j]*h[j])).<=1).*((ones(n,1)-abs((pv[i,j]-x[.,j])/(s[j]*h[j]))^3)^3);
	     j=j+1; endo;
           else;
	     j=1; do until j > d;
                wts=wts+((pv[i,j]-x[.,j])/(s[j]*h[j]))^2;
	     j=j+1; endo;
             wts=((ones(n,1)-wts)^3).*(wts.<1);
	   endif;

           nx=(ones(n,1)~dp).*sqrt(wts); 
           ny=v.*sqrt(wts);
           bv=ny/nx;
           result[i]=bv[1]; 
         
	  i=i+1; endo;
    
      else;
	  
          if (d > 1);
             if (bwrule $== "optimal");
	       h=2.4/(n^(1/(d+4)));
	       h=h*ones(d,1);
	     endif;
	  endif;

          i=1; do until i > rp;
           wts=ones(n,1);
           dp=pv[i,.]-x;
           if (bwrule $== "adaptive");
	       h=adaptbw(pv[i,.],x);
           endif;
	   if (mdkernel $== "product");
             j=1; do until j > d;
                wts=wts.*((((pv[i,j]-x[.,j])/(s[j]*h[j]))^2).<=5).*(ones(n,1)-.2*(((pv[i,j]-x[.,j])/(s[j]*h[j]))^2));
	     j=j+1; endo;
           else;
	     j=1; do until j > d;
                wts=wts+((pv[i,j]-x[.,j])/(s[j]*h[j]))^2;
	     j=j+1; endo;
             wts=(ones(n,1)-wts).*(wts.<1);
	   endif;
	 
           nx=(ones(n,1)~dp).*sqrt(wts); 
           ny=v.*sqrt(wts);
           bv=ny/nx;
           result[i]=bv[1]; 
	   
	  i=i+1; endo;

      endif;

      retp(result);

   endp;