/* Piecewise Linear & Polynomial OLS Overfit */
/* 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];
order={5,1,1};                      /* You could try order=1 or order={1,1,1} instead. */
location={1 10,
          1  5,
          6 10};
call piecewise(y,x,order,location);




proc piecewise(y,x,order,location);
   local n,z,x1,i,j,beta,xx,xs,points,npiece,xp,yp;
   n=rows(x);
   if maxc(maxc(location))>n or minc(minc(location))<1;
      errorlog "ERROR: Location index is out of range.";
      retp(".");
   endif;
   if not location[.,1]<location[.,2];
      errorlog "ERROR: Some location index is reversed or the same.";
      retp(".");
   endif;
/* sort them out in terms of x */
   z=x~y;
   z=sortc(z,1);
   x=z[.,1];
   y=z[.,2];
/* graphic settings */
   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);
      _plctrl=-1;
      _pcolor=15;
      _psymsiz=1;
      xy(x,y);
/* main */
   npiece=rows(location);
   if rows(order)==1;
      order=order*ones(npiece,1);
   endif;
   j=1;
   do while j<=npiece;
      xp=x[location[j,1]:location[j,2]];
      yp=y[location[j,1]:location[j,2]];
   /* creating X(=x1) */
      x1=ones(rows(xp),1);
      i=1;
      do while i<=order[j];
        x1=x1~xp^i;
        i=i+1;
      endo;
   /* OLS    beta=inv(X'X)*X'y */
      beta=inv(x1'x1)*x1'yp;
   /* display the coefficients */
      print/rz "data from" location[j,1] "            to" location[j,2];
      print/lz "Polynomial order=" order[j];
      print "beta=" beta;
      print;
   /* graph */
      points=200;                                                                 
      xx=seqa(minc(xp),(maxc(xp)-minc(xp))/(points-1),points-1)|maxc(xp);
      xs=ones(points,1);
      i=1;
      do while i<=order[j];
         xs=xs~xx^i;
         i=i+1;
      endo;
      setwind(1);
         _pcolor=j+2;
         _plctrl=0;
         xy(xx,xs*beta);
      j=j+1;
   endo;
   endwind;
   retp(0);
endp;