/* Henderson's 13-term Moving Average */
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\000Henderson 13";
title("Henderson 13");
xy(t,data~henderson13(data));







proc henderson13(x);
   local w,t,k,i,m,w1,w2,w3,w4,w5,w6;
/* Symmetric part in the middle */
   w={0.24006,0.21434,0.14736,0.06549,0,-0.02786,-0.01935};
   t=rows(x);
   k=rows(w)-1;
   w=rev(w[2:k+1])|w;
   m=zeros(rows(x),cols(x));
   i=k+1;
   do while i<=t-k;
      m[i,.]=w'x[i-k:i+k,.];
      i=i+1;
   endo;
/* Asymmetric part at end points */
w1={-0.09186,-0.05811,0.01202,0.11977,0.24390,0.35315,0.42113}; 
w2={-0.04271,-0.03863,0.00182,0.07990,0.17436,0.25392,0.29223,0.27910};
w3={-0.01603,-0.02487,0.00267,0.06784,0.14939,0.21605,0.24144,0.21540,0.14810};
w4={-0.00813,-0.02019,0.00413,0.06608,0.14441,0.20784,0.23002,0.20076,0.13024,0.04483};
w5={-0.01099,-0.02204,0.00330,0.06626,0.14559,0.21004,0.23324,0.20498,0.13547,0.05108,-0.01694};
w6={-0.01643,-0.02577,0.00127,0.06594,0.14698,0.21314,0.23803,0.21149,0.14368,0.06099,-0.00532,-0.03401};
   m[1,.]=rev(w1)'x[1:k+1,.]; m[t,.]=w1'x[t-k:t,.];
   m[2,.]=rev(w2)'x[1:k+2,.]; m[t-1,.]=w2'x[t-k-1:t,.];
   m[3,.]=rev(w3)'x[1:k+3,.]; m[t-2,.]=w3'x[t-k-2:t,.];
   m[4,.]=rev(w4)'x[1:k+4,.]; m[t-3,.]=w4'x[t-k-3:t,.];
   m[5,.]=rev(w5)'x[1:k+5,.]; m[t-4,.]=w5'x[t-k-4:t,.];
   m[6,.]=rev(w6)'x[1:k+6,.]; m[t-5,.]=w6'x[t-k-5:t,.];
   retp(m);
endp;