/* Nonparametric Kernel Regression (its bandwidth is manually set) */
/* 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 */
call kernreg(y,x,h,"epanechnikov");   /* select one of kernel types from                            */
                                      /* {uniform,gauss,triangular,biweight,triweight,epanechnikov} */




proc kernreg(y,x,h,types);
   local z,n,invs,i,u,kv,w,points,xx,ghat;
   z=x~y;
   z=sortc(z,1);
   x=z[.,1];
   y=z[.,2];
   n=rows(x);
   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(x)-minc(x))+minc(x)|maxc(x)+0.05*(maxc(x)-minc(x)),-0.25*(maxc(y)-minc(y))+minc(y)|maxc(y)+0.25*(maxc(y)-minc(y)));
   setwind(1);
      title("Nonparametric Kernel");
      xy(xx,ghat);
   setwind(1);
      _plctrl=-1;
      _pcolor=15;
      _psymsiz=1;
      xy(x,y);
   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;