/* Akima Polynomial 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 akimainterp(y,x,points);



proc akimainterp(y,x,points);
   local z,n,xx,yy,x1,i;
/* sort them out in terms of x */
   z=x~y;
   z=sortc(z,1);
   x=z[.,1];
   y=z[.,2];
   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]=akima(y,x,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);
      title("Akima Interpolation");
      xy(xx,yy);
   setwind(1);
      _plctrl=-1;
      _pcolor=15;
      _psymsiz=1;
      xy(x,y);
   endwind;
   retp(xx~yy);
endp;

proc akima(y,x,x0);
   local z,n,k,x1,x2,x_1,x_2,y1,y2,y_1,y_2,M1,M2,a;
   local mj,mj1,mj_1,mj_2,tj,tj1,a0,a1,a2,a3,y0;
/* sort them out in terms of x */
   z=x~y;
   z=sortc(z,1);
   x=z[.,1];
   y=z[.,2];
   n=rows(x);
   if x0<x[1] or x0>x[n];
      errorlog "ERROR: Out of range.";
      retp(".");
   endif;
   if x0==x[1];
      retp(y[1]);
   elseif x0==x[n];
      retp(y[n]);
   endif;
/* location indexing */
   k=sumc(x0.>x)+2;
/* extend endpoints */
   x1=x[n]+(x[n-1]-x[n-2]);
   x2=x[n]+(x[n]-x[n-2]);
   x_1=x[1]-(x[3]-x[2]);
   x_2=x[1]-(x[3]-x[1]);
   M1=ones(3,1)~x[n-2:n]~x[n-2:n]^2;
   M2=ones(3,1)~x[1:3]~x[1:3]^2;
   a=inv(M1)*y[n-2:n];
   y1=a[1]+a[2]*x1+a[3]*x1^2;
   y2=a[1]+a[2]*x2+a[3]*x2^2;
   a=inv(M2)*y[1:3];
   y_1=a[1]+a[2]*x_1+a[3]*x_1^2;
   y_2=a[1]+a[2]*x_2+a[3]*x_2^2;
   x=x_2|x_1|x|x1|x2;
   y=y_2|y_1|y|y1|y2;
/* coefficient calculation */
   k=k+1;
   mj=(y[k+1]-y[k])/(x[k+1]-x[k]);
   mj1=(y[k+2]-y[k+1])/(x[k+2]-x[k+1]);
   mj_1=(y[k]-y[k-1])/(x[k]-x[k-1]);
   mj_2=(y[k-1]-y[k-2])/(x[k-1]-x[k-2]);
   if abs(mj1-mj)+abs(mj_1-mj_2)==0;
      tj1=1/2*(mj_1+mj);
   else;
      tj1=(abs(mj1-mj)*mj_1+abs(mj_1-mj_2)*mj)/(abs(mj1-mj)+abs(mj_1-mj_2));
   endif;
   k=k-1;
   mj=(y[k+1]-y[k])/(x[k+1]-x[k]);
   mj1=(y[k+2]-y[k+1])/(x[k+2]-x[k+1]);
   mj_1=(y[k]-y[k-1])/(x[k]-x[k-1]);
   mj_2=(y[k-1]-y[k-2])/(x[k-1]-x[k-2]);
   if abs(mj1-mj)+abs(mj_1-mj_2)==0;
      tj=1/2*(mj_1+mj);
   else;
      tj=(abs(mj1-mj)*mj_1+abs(mj_1-mj_2)*mj)/(abs(mj1-mj)+abs(mj_1-mj_2));
   endif;
/* piecewise cubic interpolation between x[k] and x[k+1] */
   a0=y[k];
   a1=tj;
   a2=(3*mj-2*tj-tj1)/(x[k+1]-x[k]);
   a3=(tj+tj1-2*mj)/((x[k+1]-x[k])^2);
   y0=a0+a1*(x0-x[k])+a2*(x0-x[k])^2+a3*(x0-x[k])^3;
   retp(y0);
endp;