/* LOWESS */
/* Notice that this method needs enough points for small alpha. */
/* 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];
points=9;                                             /* # of points between x[k] and x[k+1] */
alpha=0.8;                                            /* 0<alpha<1 */
iter=3;                                               /* # of iterations */
call lowessinterp(y,x,points,alpha,iter);    


proc lowessinterp(y,x,points,alpha,iter);
   local z,n,xx,yy,x1,i;
/* sort them out in terms of x */
   z=x~y;
   z=sortc(z,1);
   x=z[.,1];
   y=z[.,2];
   n=rows(x);
/* graph */
   xx=zeros((points+1)*(n-1)+1,1);
   yy=zeros((points+1)*(n-1)+1,1);
   i=1;
   do while i<=n-1;
      x1=seqa(x[i],(x[i+1]-x[i])/(points+1),points+1);
      xx[(i-1)*(points+1)+1:i*(points+1)]=x1;
      i=i+1;
   endo;
   xx[rows(xx)]=x[n];
   i=1;
   do while i<=rows(yy);
      yy[i]=lowess(y,x,alpha,iter,xx[i]);
      i=i+1;
   endo;
   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("LOWESS");
      xy(xx,yy);
   setwind(1);
      _plctrl=-1;
      _pcolor=15;
      _psymsiz=1;
      xy(x,y);
   endwind;
   retp(xx~yy);
endp;

proc lowess(y,x,alpha,iter,x0);
   local z,n,k,xs,h,w,i,j,u,x1,P,beta,y0,yhat,e,r;
/* sort them out in terms of x */
   z=x~y;
   z=sortc(z,1);
   x=z[.,1];
   y=z[.,2];
   n=rows(x);
/* LOESS */
   k=floor(alpha*n);
   xs=sortc(abs(x-x0),1);
   h=xs[k];
   w=zeros(n,1);
   i=1;
   do while i<=n;
      u=abs(x[i]-x0)/h;
      if u>=0 and u<=1;
         w[i]=(1-u^3)^3;
      else;
         w[i]=0;
      endif;
      i=i+1;
   endo;
   x1=ones(n,1)~x;
   P=zeros(rows(x),rows(x)); 
   P=diagrv(P,w);
   beta=inv(x1'*P*x1)*(x1'*P*y);
   yhat=x1*beta;
/* LOWESS */
   j=1;
   do while j<=iter;
      e=y-yhat;
      r=zeros(n,1);
      i=1;
      do while i<=n;
         u=abs(e[i])/(6*median(abs(e)));
         if u>=0 and u<=1;
            r[i]=15/16*(1-u^2)^2;
         else;
            r[i]=0;
         endif;
         i=i+1;
      endo;
      P=diagrv(P,r.*w);
      beta=inv(x1'*P*x1)*(x1'*P*y);
      yhat=x1*beta;
      j=j+1;
   endo;
   y0=beta[1]+beta[2]*x0;
   retp(y0);
endp;