Simulations of time series with known ARMA structure are often needed as part of other simulations or as learning data sets for developing time series analysis skills. The following program generates a time series by using the IML functions NORMAL, ARMACOV, HANKEL, PRODUCT, RATIO, TOEPLITZ, and ROOT.

proc iml;
reset noname;
start armasim(y,n,phi,theta,seed);
/*-----------------------------------------------------------*/
/* IML Module: armasim                                       */
/* Purpose: Simulate n data points from ARMA process         */
/*          exact covariance method                          */
/* Arguments:                                                */
/*                                                           */
/* Input: n    : series length                               */
/*        phi  : AR coefficients                             */
/*        theta: MA coefficients                             */
/*        seed : integer seed for normal deviate generator   */
/* Output: y: realization of ARMA process                    */
/* ----------------------------------------------------------*/

   p=ncol(phi)-1;
   q=ncol(theta)-1;
   y=normal(j(1,n+q,seed));

      /* Pure MA or white noise */
   if p=0 then y=product(theta,y)[,(q+1):(n+q)];
   else do;                /* Pure AR or ARMA */

         /* Get the autocovariance function */
      call armacov(gamma,cov,ma,phi,theta,p);
      if gamma[1]<0 then
      do;
         print 'ARMA parameters not stable.';
         print 'Execution terminating.';
         stop;
      end;

         /* Form covariance matrix */
      gamma=toeplitz(gamma);

         /* Generate covariance between initial y and */
         /* initial innovations                       */
      if q>0 then
      do;
         psi=ratio(phi,theta,q);
         psi=hankel(psi[,-((-q):(-1))]);
         m=max(1,(q-p+1));
         psi=psi[-((-q):(-m)),];
         if p>q then psi=j(p-q,q,0)//psi;
         gamma=(gamma||psi)//(psi`||i(q));
      end;

         /* Use Cholesky root to get startup values */
      gamma=root(gamma);
      startup=y[,1:(p+q)]*gamma;
      e=y[,(p+q+1):(n+q)];

         /* Generate MA part */
      if q>0 then
      do;
         e=startup[,(p+1):(p+q)]||e;
         e=product(theta,e)[,(q+1):(n+q-p)];
      end;

      y=startup[,1:p];
      phi1=phi[,-(-(p+1):(-2))]`;

         /* Use difference equation to generate */
         /* remaining values                    */
      do ii=1 to n-p;
         y=y||(e[,ii]-y[,ii:(ii+p-1)]*phi1);
      end;
   end;
   y=y`;
finish armasim; /* ARMASIM */

run armasim(y,10,{1 -0.8},{1 0.5},1234321);
print ,'Simulated Series:', y;

The results are shown in Output 9.12.1.