/* ARMA(p,q) */
new; cls;
phi={0.7,0.1,0.1}; theta={-0.5,0.3,-0.2,0.1}; s2=0.1;         
n=1000;
x=arma(phi,theta,s2,n);
library pgraph;
graphset;
title("ARMA(p,q)");
t=seqa(1,1,n);
xy(t,x);



proc arma(phi,theta,s2,n);
   local cutn,p,q,e,x,t;
   cutn=100;                             /* first some data to cut off */
   p=rows(phi);
   q=rows(theta);
   print/lz "q=" rows(theta);;print/lz "p=" rows(phi);
   e=rndn(cutn+n,1)*sqrt(s2);
   x=zeros(cutn+n,1);
   t=maxc(p|q)+1;                              
   do while t<=cutn+n;
      x[t]=rev(phi)'x[t-p:t-1]+e[t]+rev(-theta)'e[t-q:t-1];
      t=t+1;
   endo;
   x=x[cutn+1:cutn+n];
   retp(x);
endp;