proc ged x n
type vect *x
type integer n
*
*  Draws for a generalized error distribution by acceptance-
*  rejection method.
*
*  Written March, 2000 by Tom Doan, Estima.
*
*  Syntax: @GED(nu=shape parameter) x n
*    x = vector of reals (output)
*    n = desired dimension of n
*    nu = shape parameter. Should be >=1.
*
option real nu 2.0
local integer i
local real e r
local real scalef
if nu<1.0 {
   disp '**** NU = ' nu 'is illegal. It should be >=1 ****'
   return
}
*
*  For convenience, we're using a simple one-sided exponential
*  variate as the comparison. The kernel of the target distribution
*  is exp(-x**nu). This needs to be symmetrized and scaled to
*  a variance of 1.0 at the end. The following computes the
*  needed scale factor.
*
compute scalef=exp(.5*(%lngamma(1.0/nu)-%lngamma(3.0/nu)))
dim x(n)
do i=1,n
   loop
      compute e = -log(%uniform(0.0,1.0))
      compute r = exp(-(e**nu))/exp(-e)
      if r>%uniform(0.0,1.0)
         break
   end loop
   compute x(i)=%if(%uniform(0.0,1.0)>.5,e,-e)*scalef
end do i
end

