/* Holt-Winters' Multiplicative Method where parameter alpha, beta, gam are automatically chosen. */
/* (Forecast Smoothing) Criteria: MSE */
new; cls;
let data[12,1]=
23
25
36
31
26
28
48
36
31
42
53
43
;
t=seqa(1,1,rows(data));

library pgraph;
graphset;
_pltype=6;
_plegctl=1;
_plegstr="actual\000predicted";
title("Holt-Winters' Multiplicative Method");
p=4;                                     /* number of periods in S */
prec=0.1;
xy(t,data~hwmulMSE(data,p,prec));




proc hwmulMSE(x,p,prec);
   local n,L,T,S,m,e,i,j,k,q,alpha,beta,gam,alphastar,betastar,gamstar,count,MSE,MSE_1;
   n=rows(x);
   MSE_1=1e+8;
   count=1;
   q=1;
   do while q<=(1/prec);
      k=1;
      do while k<=(1/prec);
         j=1;
         do while j<=(1/prec);
            alpha=prec*(q-1);
            beta=prec*(k-1);
            gam=prec*(j-1);
            L=zeros(n,1); T=zeros(n,1); S=zeros(n,1); m=zeros(n,1);
            S[1:p]=x[1:p]/meanc(x[1:p]);
            L[p]=x[p]/S[p];
            T[p]=0;
            i=p+1;
            do while i<=n;
               L[i]=alpha*(x[i]/S[i-p])+(1-alpha)*(L[i-1]+T[i-1]);
               T[i]=beta*(L[i]-L[i-1])+(1-beta)*T[i-1];
               S[i]=gam*(x[i]/L[i])+(1-gam)*S[i-p];
               m[i]=(L[i-1]+T[i-1])*S[i-p]; 
               i=i+1;
            endo;
            L[1:p-1]=miss(zeros(p-1,1),0);
            T[1:p-1]=miss(zeros(p-1,1),0);
            m[1:p]=miss(zeros(p,1),0);
            e=x[p+1:n]-m[p+1:n];
            MSE=sumc(e^2)/rows(e);
            if MSE<MSE_1;
               alphastar=alpha;
               betastar=beta;
               gamstar=gam;
               MSE_1=MSE;
            endif;
            call sleep(0.0001);
            format/ros 8,8;
            print/lz "iteration" count "out of" (1/prec)^3;;
            print/rz "(alpha,beta,gam)=" alpha beta gam;;
            print "         MSE=" MSE;
            count=count+1;
            j=j+1;
         endo;
         k=k+1;
      endo;
      q=q+1;
   endo;
   print;
   print/lz "alpha*=" alphastar;
   print/lz "beta* =" betastar;
   print/lz "gam*  =" gamstar;
/* main */
   L=zeros(n,1); T=zeros(n,1); S=zeros(n,1); m=zeros(n,1);
   S[1:p]=x[1:p]/meanc(x[1:p]);
   L[p]=x[p]/S[p];
   T[p]=0;
   i=p+1;
   do while i<=n;
      L[i]=alphastar*(x[i]/S[i-p])+(1-alphastar)*(L[i-1]+T[i-1]);
      T[i]=betastar*(L[i]-L[i-1])+(1-betastar)*T[i-1];
      S[i]=gamstar*(x[i]/L[i])+(1-gamstar)*S[i-p];
      m[i]=(L[i-1]+T[i-1])*S[i-p]; 
      i=i+1;
   endo;
   L[1:p-1]=miss(zeros(p-1,1),0);
   T[1:p-1]=miss(zeros(p-1,1),0);
   m[1:p]=miss(zeros(p,1),0);
   e=x[p+1:n]-m[p+1:n];
   MSE=sumc(e^2)/rows(e);
   print/lz "MSE=" MSE;
   format /mb1 /ros 16,8;
   retp(m);
endp;