/* Polynomial OLS Overfit (order=1 through 8 or 9) 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];
eliminate={1,2};
call polynols1(y,x,1,eliminate);
call polynols1(y,x,2,eliminate);
call polynols1(y,x,3,eliminate);
call polynols1(y,x,6,eliminate);



proc polynols1(y,x,order,elimination);
   local z,xall,yall,eindex,n,x1,i,beta,xx,xs,points,fmat;
/* 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);
/* creating X(=x1) */
   x1=ones(n,1);
   i=1;
   do while i<=order;
     x1=x1~x^i;
     i=i+1;
   endo;
/* OLS    beta=inv(X'X)*X'y */
   beta=inv(x1'x1)*x1'y;
/* display the coefficients */
   print/lz "Polynomial order=" order;
   print "beta=" beta;
   print;
/* graph */
   points=200;                                                                 
   xx=seqa(minc(x),(maxc(x)-minc(x))/(points-1),points-1)|maxc(x);
   xs=ones(points,1);
   i=1;
   do while i<=order;
     xs=xs~xx^i;
     i=i+1;
   endo;
   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);
      fmat="Polynomial (order= %*.*lG )";
      title(ftos(order,fmat,1,0));
      xy(xx,xs*beta);
   setwind(1);
      _plctrl=-1;
      _pcolor=15;
      _psymsiz=1;
      xy(xall,yall);
   endwind;
   retp(xx~xs*beta);
endp;