/* Nonparametric Kernel Regression (its bandwidth is manually set) eliminating some points */
/* Data from Holton(2003)"Value-at-Risk" */
new; cls;
let data[10,2]=
1.1  2.14
1.4  2.60
2.5  1.15
2.7  1.19
3.2  1.88
3.6  1.55
4.1  2.65
4.3  3.80
4.5  4.46
4.9  6.35
;
x=data[.,1]; y=data[.,2];
h=0.5;                                          /* manually set */
eliminate={1,2};
call kernreg1(y,x,h,"epanechnikov",eliminate);  /* select one of kernel types from                            */
                                                /* {uniform,gauss,triangular,biweight,triweight,epanechnikov} */




proc kernreg1(y,x,h,types,eliminate);
   local z,xall,yall,eindex,n,invs,i,u,kv,w,points,xx,ghat;
/* sort them out in terms of x */
   z=x~y;
   z=sortc(z,1);
   x=z[.,1];
   y=z[.,2];
/* elimination */
   if maxc(eliminate)>rows(x);
      errorlog "ERROR: Out of range.";
      retp(".");
   elseif eliminate==0 or eliminate==zeros(rows(eliminate),1);
      errorlog "No elimination:";
      xall=x;
      yall=y;
   else;
      xall=x;
      yall=y;
      eindex=zeros(rows(x),1);
      eindex[eliminate]=ones(rows(eliminate),1);
      x=delif(x,eindex);
      y=delif(y,eindex);
   endif;
   n=rows(x);
/* main */
   points=1000;
   xx=seqa(minc(x),(maxc(x)-minc(x))/(points-1),points-1)|maxc(x);
   invs=inv(chol(vcx(x)));
   ghat=zeros(points,1);
   i=1; 
   do while i<=points;
      u=(xx[i]-x)./h;
      kv=kernels(u*invs,types);
      w=kv./sumc(kv);
      ghat[i]=w'y;
      i=i+1;
   endo;
/* graph */
   library pgraph;
   graphset;
   pqgwin auto;
   begwind;
   window(1,1,0);
   scale(-0.05*(maxc(xall)-minc(xall))+minc(xall)|maxc(xall)+0.05*(maxc(xall)-minc(xall)),-0.25*(maxc(yall)-minc(yall))+minc(yall)|maxc(yall)+0.25*(maxc(yall)-minc(yall)));
   setwind(1);
      title("Nonparametric Kernel");
      xy(xx,ghat);
   setwind(1);
      _plctrl=-1;
      _pcolor=15;
      _psymsiz=1;
      xy(xall,yall);
   endwind;
   retp(xx~ghat);
endp;

proc kernels(u,types);
   local kv;
   if types$=="uniform";
      kv=0.5*(abs(u).<=1); 
   elseif types$=="gauss";              
      kv=(1/sqrt(2*pi))*exp(-0.5*u^2);
   elseif types$=="triangular";
      kv=(1-abs(u)).*(abs(u).<=1); 
   elseif types$=="biweight";         
      kv=(15/16)*((1-u^2)^2).*(abs(u).<=1);
   elseif types$=="triweight";
      kv=(35/32)*((1-u^2)^3).*(abs(u).<=1);
   elseif types$=="epanechnikov";
      kv=0.75*(1-u^2).*(abs(u).<=1);
   else;
      errorlog "ERROR:  Type must be either of the following:";
      errorlog "uniform/gauss/triangular/biweight/triweight/epanechnikov";
   endif;
   retp( kv );
endp;