/* Newton Interpolation 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];
points=9;                                             /* # of points between x[k] and x[k+1] */
eliminate={1,2};
call newtoninterp1(y,x,points,eliminate);





proc newtoninterp1(y,x,points,eliminate);
   local z,xall,yall,eindex,n,i,b,x1,xx,yy;
/* 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);
/* 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]=newtonint(y,x,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(yy)-minc(yy))+minc(yy)|maxc(yy)+0.25*(maxc(yy)-minc(yy)));
   setwind(1);
      title("Newton Interpolation");
      xy(xx,yy);
   setwind(1);
      _plctrl=-1;
      _pcolor=15;
      _psymsiz=1;
      xy(xall,yall);
   endwind;
   retp(xx~yy);
endp;

proc newtonint(y,x,x0);
   local z,n,i,j,k,num,denom,s,c,y0;
/* sort them out in terms of x */
   z=x~y;
   z=sortc(z,1);
   x=z[.,1];
   y=z[.,2];
   n=rows(x);
/* calculation of c */
   c=zeros(n,1);
   c[1]=y[1];
   j=2;
   do while j<=n;
      denom=1;
      i=1;
      do while i<=j-1;
         denom=denom*(x[j]-x[i]);
         i=i+1;
      endo;
      if j==2;
         num=y[2]-c[1];
      else;
         num=y[j]-c[1];
         i=3;
         do while i<=n;
            s=c[i-1]; 
            k=1;
            do while k<=i-2;
               s=s*(x[j]-x[k]);
               k=k+1;
            endo;
            num=num-s;
            i=i+1;
         endo;
      endif;
      c[j]=num/denom;
      j=j+1;
   endo;
/* polynomial f(x)=c1+c2*(x-x1)+c3*(x-x1)*(x-x2)+c4*(x-x1)*(x-x2)*(x-x3)+... evaluated at x=x0 */
   y0=c[1];
   i=2;
   do while i<=n;
      s=1;
      k=1;
      do while k<=i-1;
         s=s*(x0-x[k]);
         k=k+1;
      endo;
      y0=y0+c[i]*s;
      i=i+1;
   endo;
   retp(y0);
endp;