/* B-Spline 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];
allpoints=100;                                  /* # of all points between min(x) and max(x) */
eliminate={1,2};
call bsplineinterp1(y,x,allpoints,eliminate);





proc bsplineinterp1(y,x,allpoints,eliminate);
   local z,xall,yall,eindex,i,xx,yy,t,yi,xi;
/* 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;
/* t step (-1<=t<=rows(x)) */
   t=seqa(-1,(rows(x)+1)/(allpoints-1),allpoints);
   t[rows(t)]=rows(x);
/* graph */
   xx=zeros(allpoints,1);
   yy=zeros(allpoints,1);
   i=1;
   do while i<=rows(yy);
      {yi,xi}=bspline(y,x,t[i]);
      yy[i]=yi;
      xx[i]=xi;
      i=i+1;
   endo;
/* graph */
   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);
      title("B-Spline Interpolation");
      xy(xx,yy);
   setwind(1);
      _plctrl=-1;
      _pcolor=15;
      _psymsiz=1;
      xy(xall,yall);
   endwind;
   retp(xx~yy);
endp;

proc(2)=bspline(y,x,t);
   local z,n,w,j,k,yi,xi;
   if t<-1 or t>rows(x);
      errorlog "ERROR: Out of range.";
   endif;
/* sort them out in terms of x */
   z=x~y;
   z=sortc(z,1);
   x=z[.,1];
   y=z[.,2];
/* B-spline interpolation */
   n=rows(x)-1;
   xi=0;
   yi=0;
   j=-2;
   do while j<=n+2;
      k=j;
      /* endpoint index is 0 or n */
      if j<0;
         k=0;
      endif;
      if j>n;
         k=n;
      endif;
      /* weight */
      if abs(t-j)<1;
         w=(3*abs(t-j)^3-6*abs(t-j)^2+4)/6;
      elseif abs(t-j)<2;
         w=-(abs(t-j)-2)^3/6;
      else;
         w=0;
      endif;
      /* iterated sum */
      xi=xi+w*x[k+1];
      yi=yi+w*y[k+1];
      j=j+1;
   endo;
   retp(yi,xi);
endp;