## 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;
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"];