/* Cubic Spline Interpolation */
/* 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] */
call cubspline(y,x,points);





proc cubspline(y,x,points);
   local z,n,p,M,i,b,x1,xs,xx,yy;
/* sort them out in terms of x */
   z=x~y;
   z=sortc(z,1);
   x=z[.,1];
   y=z[.,2];
   n=rows(x);
/* p column vector where M*b=p*/
   p=zeros((n-1)*4,1);
   p[1]=y[1];
   p[(n-1)*2]=y[n];
   i=2;
   do while i<=n-1;
      p[(i-1)*2:(i-1)*2+1]=y[i]|y[i];
      i=i+1;
   endo;
/* M */
   M=zeros((n-1)*4,(n-1)*4);
   /* (i) each polynomial passing through its respective end points */
   i=1;
   do while i<=n-1;
      M[(i-1)*2+1,(i-1)*4+1:(i-1)*4+4]=x[i]^3~x[i]^2~x[i]~1;
      M[(i-1)*2+2,(i-1)*4+1:(i-1)*4+4]=x[i+1]^3~x[i+1]^2~x[i+1]~1;
      i=i+1;
   endo;
   /* (ii) 1st derivatives matching at interior points */
   i=2;
   do while i<=n-1;
      M[(i-1)+2*(n-1),(i-2)*4+1:(i-2)*4+8]=3*x[i]^2~2*x[i]~1~0~-3*x[i]^2~-2*x[i]~-1~0;
      i=i+1;
   endo;
   /* (iii) 2nd derivatives matching at interior points */
   i=2;
   do while i<=n-1;
      M[(i-1)+2*(n-1)+(n-2),(i-2)*4+1:(i-2)*4+8]=6*x[i]~2~0~0~-6*x[i]~-2~0~0;
      i=i+1;
   endo;
   /* (iv) 2nd derivatives vanishing at the end points */
      M[(n-1)*4-1,1:4]=6*x[1]~2~0~0;
      M[(n-1)*4,(n-1)*4-3:(n-1)*4]=6*x[n]~2~0~0;
/* M*b=p => b=inv(M)*p */
   b=inv(M)*p;
/* display the coefficients */
   i=1;
   do while i<=n-1;
      print/lz "p" i "(x)      =";;
      print b[(i-1)*4+1] "x^3     +" b[(i-1)*4+2] "x^2     +" b[(i-1)*4+3] "x     +" b[(i-1)*4+4];
      i=i+1;
   endo;
/* 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;
      xs=x1^3~x1^2~x1~ones(points+1,1);
      yy[(i-1)*(points+1)+1:i*(points+1)]=xs*b[(i-1)*4+1:i*4];
      i=i+1;
   endo;
   xx[rows(xx)]=x[n];
   yy[rows(yy)]=y[n];
   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);
      title("Cubic Spline");
      xy(xx,yy);
   setwind(1);
      _plctrl=-1;
      _pcolor=15;
      _psymsiz=1;
      xy(x,y);
   endwind;
   retp(xx~yy);
endp;