/* Polynomial Time Series Moving Average PTSMA(k,order) */
new; cls;
let data[50,1]=
      60
      61
      60
      62
      60
      63
      64
      65
      69
      70
      69
      66
      64
      64
      62
      64
      65
      55
      53
      57
      60
      58
      63
      63
      63
      68
      69
      69
      68
      69
      70
      67
      66
      67
      68
      69
      69
      77
      80
      80
      78
      76
      76
      72
      71
      70
      73
      71
      78
      80
;
data=rev(data);               /* original data in reverse order */
t=seqa(1,1,rows(data));

library pgraph;
graphset;
_pltype=6;
_plegctl=1;
_plegstr="actual\000Cubic TSMA(10)";
title("Polynomial Time Series Moving Average");
xy(t,data~ptsma(data,10,3));




proc ptsma(x,k,order);
   local m,i,j,x1,y1,b,yhat;
   x1=ones(k,1);
   i=1;
   do while i<=order;
     x1=x1~seqa(1,1,k)^i;
     i=i+1;
   endo;
   m=zeros(rows(x),cols(x));
   j=1;
   do while j<=cols(x);
      i=k;
      do while i<=rows(x);
         /* OLS by every k data */
         y1=x[i-k+1:i,j];
         b=inv(x1'x1)*x1'y1;
         yhat=x1*b;
         /* Get each current estimate */
         m[i,j]=yhat[k];
         i=i+1;
      endo;
      j=j+1;
   endo;
   m[1:k-1,.]=miss(m[1:k-1,.],0);
   retp(m);
endp;