/* Smoothing Natural Spline Interpolation */
/* Notice: Some parameter relations may lead to no solution. Direct matrix inverse applied here. */
/* 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] */
k=5;                                                  /* k=2(linear),3,4,5,... */
w=ones(rows(x),1);         /* (positive)weight vector e.g. w={1,1,...,1} */
g=100;                     /* some constant */
call snaturalinterp(y,x,points,k,w,g);    



proc snaturalinterp(y,x,points,k,w,g);
   local z,n,xx,yy,x1,i,fmat,om;
/* sort them out in terms of x */
   z=x~y~w;
   z=sortc(z,1);
   x=z[.,1];
   y=z[.,2];
   w=z[.,3];
   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]=snatural(y,x,k,w,g,xx[i]);
      i=i+1;
   endo;
   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="Smoothing Natural Spline Interpolation k= %*.*lg";
      om=ftos(k,fmat,1,2);
      title(om);
      xy(xx,yy);
   setwind(1);
      _plctrl=-1;
      _pcolor=15;
      _psymsiz=1;
      xy(x,y);
   endwind;
   retp(xx~yy);
endp;

proc snatural(y,x,k,w,g,x0);
   local z,n,A,i,j,b,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 matrix */
   A=zeros(n+k,n+k);
   A[1:n,1:k]=x^seqa(0,1,k)';
   A[n+1:n+k,k+1:n+k]=(x^seqa(0,1,k)')';
   i=1;
   do while i<=n;
      j=1;
      do while j<=i;
         A[i,k+j]=(x[i]-x[j])^(2*k-1);
         j=j+1;
      endo;
      i=i+1;
   endo;
   i=1;
   do while i<=n;
      A[i,k+i]=(-1)^k*g*(2*k-1)!/w[i];
      i=i+1;
   endo;
   b=y|zeros(k,1);
   c=inv(A)*b;
/* evaluated at x=x0 */
   y0=c[1];
   i=2;
   do while i<=k;
      y0=y0+c[i]*x0^(i-1);
      i=i+1;
   endo;
   i=k+1;
   do while i<=n+k;
      y0=y0+c[i]*maxc(((x0-x[i-k])^(2*k-1))|0);
      i=i+1;
   endo;
   retp(y0);
endp;