/* ARMA-GARCH */
new; cls;
phi={0.8,0.1}; theta={0.5,0.3,0.2,0.1};
a0=0.001;
a={0.2,0.1,0.1}; b={0.2,0.1,0.1,0.1};          
n=1000;
x=armagarch(phi,theta,a0,a,b,n);
library pgraph;
graphset;
title("ARMA-GARCH");
t=seqa(1,1,n);
xy(t,x);



proc armagarch(phi,theta,a0,a,b,n);
   local cutn,p,q,e,x,t;
   cutn=100;                             /* first some data to cut off */
   p=rows(phi);
   q=rows(theta);
   print "ARMA parameter:";
   print/lz "q=" rows(theta);;print/lz "p=" rows(phi);
   e=garch(a0,a,b,cutn+n);                            /* Now, e~GARCH. */
   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;

proc garch(a0,a,b,n);
   local cutn,e,s2,q,p,u,h,max,t;
   if a0<=0;
      errorlog "ERROR:Parameter a0 must be positive.";
      retp(".");
   endif;
   if not a>=0;
      errorlog "ERROR:Each element of a must be non negative.";
      retp(".");
   endif;
   if not b>=0;
      errorlog "ERROR:Each element of b must be non negative.";
      retp(".");
   endif;
   cutn=100;                  /* first some data to cut off */
   e=rndn(cutn+n,1);   
   s2=a0/(1-sumc(a)-sumc(b));
   q=rows(a); p=rows(b);
   print "GARCH parameter:";
   print/lz "q=" q;; print/lz "p=" p;
   u=zeros(cutn+n,1);
   h=zeros(cutn+n,1);
   max=maxc(q|p);
   u[1:max]=sqrt(s2)*e[1:max]; h[1:max]=s2*ones(max,1);
   t=max+1;
   do while t<=cutn+n;
      h[t]=a0+rev(a)'(u[t-q:t-1]^2)+rev(b)'h[t-p:t-1];
      u[t]=sqrt(h[t])*e[t];
      t=t+1;
   endo;
   u=u[cutn+1:cutn+n];
   retp(u);
endp;