/* Polynomial OLS Overfit (order=1 through 8 or 9) passing the point (x0,y0) */
/* 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]; x0=2; y0=2;
call fixedols(y,x,y0,x0,1);
call fixedols(y,x,y0,x0,2);
call fixedols(y,x,y0,x0,3);
call fixedols(y,x,y0,x0,5);



proc fixedols(y,x,y0,x0,order);
   local n,x1,beta,xx,xs,points,fmat,b,j,i,nn;
/* OLS    beta=inv(X'X)*X'y going through (x0,y0) */
   n=rows(x);
   x1=(x-x0)^seqa(1,1,order)';
   beta=inv(x1'x1)*x1'(y-y0);
/* expansion of polynomial to get the coefficients */
   b=zeros(order+1,order);
   j=1;
   do while j<=order;
      i=0; nn=order-(j-1);
      do while i<=nn;
         b[i+j,j]=beta[order-j+1]*nn!/((nn-i)!*i!)*(-x0)^i;
         i=i+1;
      endo;
      j=j+1;
   endo;
   b=rev(sumc(b')); 
   b[1]=b[1]+y0;
/* display the coefficients */
   print/lz "Polynomial order=" order;
   print "beta=" b;
   print;
/* graph */
   points=200;                                                                 
   xx=seqa(minc(x),(maxc(x)-minc(x))/(points-1),points-1)|maxc(x);
   xs=(xx-x0)^seqa(1,1,order)';
   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);
      fmat="Polynomial (order= %*.*lG )";
      title(ftos(order,fmat,1,0));
      xy(xx,xs*beta+y0);
   setwind(1);
      _plctrl=-1;
      _pcolor=15;
      _psymsiz=1;
      xy(x,y);
   setwind(1);
      _pcolor=12;
      xy(x0|miss(zeros(n-1,1),0),y0|miss(zeros(n-1,1),0));
   endwind;
   retp(xx~xs*beta+y0);
endp;