/* This program generates the ACF and PACF for an 
                MA(1)*SMA(1) process (s=12) */

data b;
array r{40};
array q(40,40);

 do k=1 to 40;
   theta1=-0.8; /*you can input the values of two parameters here*/ 
   Theta=-0.8;  /* theta1 = regular MA parameter, Theta = seasonal
                                                     MA parameter */
   if k=1 then r{k}=-theta1/(1+theta1**2);
   else if k=11 then r{k}=(theta1*Theta)/((1+theta1**2)*(1+Theta**2));
   else if k=12 then r{k}=-Theta/(1+Theta**2);
   else if k=13 then r{k}=(theta1*Theta)/((1+theta1**2)*(1+Theta**2));
   else r{k}=0;
   rho=r{k};
              do j= k to 1 by -1;
                   if k=1 then q(1,1)=r{k};
				   else;
				   do;
                        sum1=0;
	                    sum2=0;
	                    do i=1 to k-1;
	                        subsum1=q(k-1,i)*r(k-i);
		                    sum1=sum1+subsum1;
		                    subsum2=q(k-1,i)*r(i);
	                     	sum2=sum2+subsum2;
	                    end;

                        if k=j then 
                                 do;
                                 q(k,j)=(r(k)-sum1)/(1-sum2);
								 thi=q(k,j);
                                 end;
	                    else q(k,j)=q(k-1,j)-q(k,k)*q(k-1,k-j);
                    end;
             end;
	
	 output;
  
end;
run;

data b;
  set b;
  lag=k;
  phijj=thi;


symbol interpol=needle
       cv=blue
	   ci=black
	   value=dot
	   height=1
	   width=1;

proc gplot data=b;
  title1 'Theoretical ACF for MA(1)*SMA(1) Process';
  title2 'with MA1 = -0.8 and SMA1 = -0.8 (S=12)';
  plot rho*lag;   
  run;
  title1 'Theoretical PACF for MA(1)*SMA(1) Process';
  title2 'with MA1 = -0.8 and SMA1 = -0.8 (S=12)';
  plot phijj*k;
run;