/* This program plots the ACF and PACF for an ARMA(1,1) process.  
   Phi1 is the autoregressive coefficient.  Theta1 is the 
   moving average coefficient.   For stationarity -1 < phi1 < 1.
   For invertibility -1 < theta1 < 1.  Also to avoid the 
   common factor problem phi1 and theta1 cannot equal each other.
   The user can, within these constraints, choose his or her own
   values of the parameters to see how the coefficient settings
   affect the appearances of the ACF and PACF graphs.  Also by
   setting theta1 = 0 you can examine the ACF and PACF of the
   AR(1) model while by setting phi1 = 0 you can examine the
   ACF and PACF of the MA(1) model.  */

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

   phi1=-0.5;/* You can input the values of the two coefficients here. */
   theta1=0.5;

  if phi1 = 0.0 then
    do lag = 1 to 30;
    if lag = 1 then r{lag}=-theta1/(1+theta1**2);
    else r{lag}=0;
    rho=r{lag};

              do j= lag to 1 by -1;
                   if lag=1 then q(1,1)=r{lag};
				   else;
				   do;
                        sum1=0;
	                    sum2=0;
	                    do i=1 to lag-1;
	                        subsum1=q(lag-1,i)*r(lag-i);
		                    sum1=sum1+subsum1;
		                    subsum2=q(lag-1,i)*r(i);
	                     	sum2=sum2+subsum2;
	                    end;

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

   else;

 do lag = 1 to 30;
    
   r{lag}=((1-theta1*phi1)*(phi1-theta1)*phi1**(lag-1))/(1-2*theta1*phi1+theta1**2);
   rho=r{lag};
              do j= lag to 1 by -1;
                   if lag=1 then q(1,1)=r{lag};
				   else;
				   do;
                        sum1=0;
	                    sum2=0;
	                    do i=1 to lag-1;
	                        subsum1=q(lag-1,i)*r(lag-i);
		                    sum1=sum1+subsum1;
		                    subsum2=q(lag-1,i)*r(i);
	                     	sum2=sum2+subsum2;
	                    end;

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

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

proc gplot data=b;
  title1 'Theoretical ACF for ARMA(1,1) Process';
  title2 'with phi1 = -0.5 and theta1 = 0.5'; 
  plot rho*lag;   
  run;
  title1 'Theoretical PACF for ARMA(1,1) Process';
  title2 'with phi1 = -0.5 and theta1 = 0.5';
  plot phijj*lag;
run;