Simulations of time series with known autoregressive moving average (ARMA) structure are often needed as part of other simulations or as sample data sets for developing skills in time series analysis. You can use the ARMASIM function to simulate a univariate series from an ARMA model. The following module shows some of the computations that are required to simulate data from an ARMA model. The module uses many SAS/IML functions, including the ARMACOV, HANKEL, PRODUCT, RATIO, TOEPLITZ, and ROOT functions. A short simulated ARMA(1,1) series is shown in Output 9.12.1.

proc iml;
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.',
              '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;

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

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

end;
y = y`;
finish armasim; /* ARMASIM */

run armasim(y,10,{1 -0.8},{1 0.5}, 1234321);
print y[label="Simulated Series"];