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;