/* Brown's one-Prameter quadratric */
/* (Forecast Smoothing) Criteria: MAPE */
new; cls;
let data[12,1]=
23
25
36
31
26
28
48
36
31
42
53
43
;
nahead=10;
t=seqa(1,1,rows(data)+nahead);

library pgraph;
graphset;
_pltype=6;
_plegctl=1;
_plegstr="actual\000predicted\000n ahead";
title("QUAD");
prec=0.01;
xy(t,(data|miss(zeros(nahead,1),0))~quadMAPE(data,nahead,prec));




proc quadMAPE(x,nahead,prec);
   local n,m1,m2,m3,m,mp,e,i,j,b,c,t,alpha,alphastar,count,MAPE,MAPE_1;
   n=rows(x);
   MAPE_1=1e+8;
   count=1;
      j=1;
      do while j<=(1/prec);
         alpha=prec*(j-1);
         m1=zeros(n,1);
         m1[1]=x[1];
         i=2;
         do while i<=n;
            m1[i]=alpha*x[i-1]+(1-alpha)*m1[i-1]; 
            i=i+1;
         endo;
         m2=zeros(n,1);
         m2[1]=m1[1];
         i=2;
         do while i<=n;
            m2[i]=alpha*m1[i-1]+(1-alpha)*m2[i-1]; 
            i=i+1;
         endo;
         m3=zeros(n,1);
         m3[1]=m2[1];
         i=2;
         do while i<=n;
            m3[i]=alpha*m2[i-1]+(1-alpha)*m3[i-1]; 
            i=i+1;
         endo;
         m=3*m1-3*m2+m3;
         e=x[2:n]-m[2:n];
         MAPE=sumc(abs(e)./x[2:n]*100)/rows(e);
         if MAPE<MAPE_1;
            alphastar=alpha;
            MAPE_1=MAPE;
         endif;
         call sleep(0.0001);
         print/lz "iteration" count "out of" (1/prec);;
         print/rz "alpha=" alpha;;
         print "         MAPE=" MAPE;
         count=count+1;
         j=j+1;
      endo;
   print;
   print/lz "alpha*=" alphastar;
/* main */
   m1=zeros(n,1);
   m1[1]=x[1];
   i=2;
   do while i<=n;
      m1[i]=alphastar*x[i-1]+(1-alphastar)*m1[i-1]; 
      i=i+1;
   endo;
   m2=zeros(n,1);
   m2[1]=m1[1];
   i=2;
   do while i<=n;
      m2[i]=alphastar*m1[i-1]+(1-alphastar)*m2[i-1]; 
      i=i+1;
   endo;
   m3=zeros(n,1);
   m3[1]=m2[1];
   i=2;
   do while i<=n;
      m3[i]=alphastar*m2[i-1]+(1-alphastar)*m3[i-1]; 
      i=i+1;
   endo;
   m=3*m1-3*m2+m3;
   b=alphastar^2/(2*(1-alphastar)^2)*((6-5*alphastar)*m1[n]-(10-8*alphastar)*m2[n]+(4-3*alphastar)*m3[n]);
   c=alphastar^2/((1-alphastar)^2)*(m1[n]-2*m2[n]+m3[n]);
   t=seqa(0,1,nahead+1);
   mp=miss(zeros(n-1,1),0)|m[n]+0.5*b*t+c*t^2;
   m=m|miss(zeros(nahead,1),0);
   e=x[2:n]-m[2:n];
   MAPE=sumc(abs(e)./x[2:n]*100)/rows(e);
   print/lz "MAPE=" MAPE;
   retp(m~mp);
endp;