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