/* Multiplicative Model(seasonal adjustment by dropping max and min values from SR) */
new; cls;
let data[91,2]=
198001 4958.9
198002 4857.8
198003 4850.3
198004 4936.6
198101 5032.5
198102 4997.3
198103 5056.8
198104 4997.1
198201 4914.3
198202 4935.5
198203 4912.1
198204 4915.6
198301 4972.4
198302 5089.8
198303 5180.4
198304 5286.8
198401 5402.3
198402 5493.8
198403 5541.3
198404 5583.1
198501 5629.7
198502 5673.8
198503 5758.6
198504 5806.0
198601 5858.9
198602 5883.3
198603 5937.9
198604 5969.5
198701 6013.3
198702 6077.2
198703 6128.1
198704 6234.4
198801 6275.9
198802 6349.8
198803 6382.3
198804 6465.2
198901 6543.8
198902 6579.4
198903 6610.6
198904 6633.5
199001 6716.3
199002 6731.7
199003 6719.4
199004 6664.2
199101 6631.4
199102 6668.5
199103 6684.9
199104 6720.9
199201 6783.3
199202 6846.8
199203 6899.7
199204 6990.6
199301 6988.7
199302 7031.2
199303 7062.0
199304 7168.7
199401 7229.4
199402 7330.2
199403 7370.2
199404 7461.1
199501 7488.7
199502 7503.3
199503 7561.4
199504 7621.9
199601 7676.4
199602 7802.9
199603 7841.9
199604 7931.3
199701 8016.4
199702 8131.9
199703 8216.6
199704 8272.9
199801 8396.3
199802 8442.9
199803 8528.5
199804 8667.9
199901 8733.2
199902 8775.5
199903 8886.9
199904 9040.1
200001 9097.4
200002 9205.7
200003 9218.7
200004 9243.8
200101 9229.9
200102 9193.1
200103 9186.4
200104 9248.8
200201 9363.2
200202 9392.4
200203 9485.6
;
y=data[.,2];              /* y=y[2:rows(y),.]-y[1:rows(y)-1,.]; */
t0=1995; step=1/4; N=4;
call multmodel(y,t0,step,N);






proc(4)=multmodel(x,t0,step,N);
   local nx,MA,i,j,SR,SS,Sx,x1,b,count,S,y,T,C,R,label;
   nx=rows(x);
/* MA calculation */
   MA=miss(zeros(nx,1),0);
   if (N-1)%2;            /* When N is odd. */
      i=N/2+1;
      do while i<=nx-N/2;
         MA[i]=(sumc(x[i-N/2:i+N/2-1])/N+sumc(x[i-N/2+1:i+N/2])/N)/2;
         i=i+1;
      endo;
   else;                  /* When N is even. */
      i=(N-1)/2+1;
      do while i<=nx-(N-1)/2;
         MA[i]=sumc(x[i-(N-1)/2:i+(N-1)/2])/N;
         i=i+1;
      endo;
   endif;
/* S calculation */
   SR=x./MA;
   SS=zeros(N,1);
   if (N-1)%2;
      j=1;
      do while j<=N;
         Sx=miss(0,0);
         i=N/2+1;
         do while i+(j-1)<=nx-N/2;
            Sx=Sx|SR[i+(j-1)];
            i=i+N;
         endo;
         Sx=Sx[2:rows(Sx)];
         Sx=delif(Sx,Sx.==maxc(Sx));
         Sx=delif(Sx,Sx.==minc(Sx));
         SS[j]=meanc(Sx);
         j=j+1;
      endo;
      S=miss(zeros(nx,1),0);
      j=1;
      do while j<=N;
         i=N/2+1;
         do while i+(j-1)<=nx-N/2;
            S[i+(j-1)]=SS[j];
            i=i+N;
         endo;
         j=j+1;
      endo;
   else;
      j=1;
      do while j<=N;
         Sx=miss(0,0);
         i=(N-1)/2+1;
         do while i+(j-1)<=nx-(N-1)/2;
            Sx=Sx|SR[i+(j-1)];
            i=i+N;
         endo;
         Sx=Sx[2:rows(Sx)];
         Sx=delif(Sx,Sx.==maxc(Sx));
         Sx=delif(Sx,Sx.==minc(Sx));
         SS[j]=meanc(Sx);
         j=j+1;
      endo;
      S=miss(zeros(nx,1),0);
      j=1;
      do while j<=N;
         i=(N-1)/2+1;
         do while i+(j-1)<=nx-(N-1)/2;
            S[i+(j-1)]=SS[j];
            i=i+N;
         endo;
         j=j+1;
      endo;   
   endif;
/* T calculation */  
   y=delif(MA,MA.==miss(0,0));
   x1=ones(rows(y),1)~seqa(1,1,rows(y));
   b=inv(x1'x1)*x1'y;
   T=x1*b;
   if (N-1)%2;
      T=miss(zeros(N/2,1),0)|T|miss(zeros(N/2,1),0);
   else;
      T=miss(zeros((N-1)/2,1),0)|T|miss(zeros((N-1)/2,1),0);
   endif;
/* C & R calculation */
   C=MA./T;
   R=x./(S.*T.*C);
/* display the result */
   print/lz "N=" N;
   label={"x" "S" "T" "C" "R"};
   print $label;;
   print x~S~T~C~R;
/* graph */
   library pgraph;
   graphset;
   begwind;
   window(4,1,1);
   setwind(1);
      _plegctl=1;
      _plegstr="data\000MA\000T";
      _pltype=6;
      ylabel("data, MA, T");
      xy(seqa(t0,step,nx),x~MA~T);
   setwind(2);
      graphset;
      ylabel("C");
      xy(seqa(t0,step,nx),C~ones(nx,1));
   setwind(3);
      graphset;
      ylabel("S");
      xy(seqa(t0,step,nx),S);
   setwind(4);
      graphset;
      ylabel("R");
      xy(seqa(t0,step,nx),R~ones(nx,1));
   endwind;
   retp(S,T,C,R);
endp;