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



proc armaegarch_GED(phi,theta,a0,a,b,nu,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=egarch_GED(a0,a,b,gam,nu,cutn+n);           /* Now, e~EGARCH_GED. */
   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 egarch_GED(a0,a,b,gam,nu,n);
   local cutn,e,s2,q,p,u,lnsig2,max,t,sig;
   if a0<=0;
      errorlog "ERROR:Parameter a0 must be positive.";
      retp(".");
   endif;
   if rows(a)/=rows(gam);
      errorlog "ERROR: The number of element a and its gam must be the same.";
      retp(".");
   endif;
   if rows(a)<rows(b);
      errorlog "Warning: If p>q, the period p-q is ignored for z=u/sig calculation";
   endif;
   cutn=100;                  /* first some data to cut off */
   e=rndged0(nu,cutn+n);                                      /* Now, e~GED(nu). */  
   s2=a0/(1-sumc(a)-sumc(b));
   q=rows(a); p=rows(b);
   print "EGARCH parameter:";
   print/lz "q=" q;; print/lz "p=" p;
   u=zeros(cutn+n,1);
   lnsig2=zeros(cutn+n,1);
   max=maxc(q|p);
   u[1:max]=sqrt(s2)*e[1:max]; lnsig2[1:max]=ln(s2)*ones(max,1);
   t=max+1;
   do while t<=cutn+n;
      sig=sqrt(exp(lnsig2[t-q:t-1]));
      lnsig2[t]=a0+rev(a)'((abs(u[t-q:t-1])+rev(gam).*u[t-q:t-1])./sig)+rev(b)'lnsig2[t-p:t-1];
      u[t]=sqrt(exp(lnsig2[t]))*e[t];
      t=t+1;
   endo;
   u=u[cutn+1:cutn+n];
   retp(u);
endp;


proc rndged0(nu,n);
   local lowerbound,upperbound,max,x,i,y,U;
   lowerbound=-4; upperbound=4;            /* You could change, but normally O.K. */
   max=pdfged(0,nu);
/* acceptance-rejection method */
   x=zeros(n,1);
   i=1;
   do while i<=n;
      y=(upperbound-lowerbound)*rndu(1,1)+lowerbound;
      U=rndu(1,1);
      if U<=pdfged(y,nu)/max;
         x[i]=y;
         i=i+1;
      endif;
   endo;
   retp(x);
endp;

proc pdfged(x,nu);
   local lambda;
   if nu<=0;
      errorlog "ERROR: Parameter nu must be positive.";
      retp(".");
   endif;
   lambda=1/2^(1/nu).*(gamma(1/nu)./gamma(3/nu))^(1/2);
   retp( nu.*exp(-1/2*abs(x./lambda)^nu)./(lambda.*2^(1+1/nu).*gamma(1/nu)) );
endp;