/* ARCH(q) */
new; cls;
a0=0.001;
a={0.3,0.2,0.1,0.1,0.1,0.1};               
n=1000;
u=arch(a0,a,n);
library pgraph;
graphset;
title("ARCH(q)");
t=seqa(1,1,n);
xy(t,u);



proc arch(a0,a,n);
   local cutn,e,s2,q,u,h,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;
   cutn=100;                  /* first some data to cut off */
   e=rndn(cutn+n,1);          
   s2=a0/(1-sumc(a));
   q=rows(a);
   print/lz "q=" rows(a);
   u=zeros(cutn+n,1);
   h=zeros(cutn+n,1);
   u[1:q]=sqrt(s2)*e[1:q];
   t=q+1;
   do while t<=cutn+n;
      h[t]=a0+rev(a)'(u[t-q:t-1]^2);
      u[t]=sqrt(h[t])*e[t];
      t=t+1;
   endo;
   u=u[cutn+1:cutn+n];
   retp(u);
endp;