Simulating a Univariate ARMA Process

/****************************************************************/
/*          S A S   S A M P L E   L I B R A R Y                 */
/*                                                              */
/*    NAME: ARMASIM                                             */
/*   TITLE: Simulating a Univariate ARMA Process                */
/* PRODUCT: IML                                                 */
/*  SYSTEM: ALL                                                 */
/*    KEYS: MATRIX  TIME    SUGI6                               */
/*   PROCS: IML                                                 */
/*    DATA:                                                     */
/*                                                              */
/* SUPPORT: TJW                         UPDATE: 29APR87         */
/*     REF:                                     SEP2013         */
/*                                                              */
/****************************************************************/


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;
call randseed(seed, 1);
y = randfun(1 || (n+q), "Normal");

/* 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"];