/* Stable Random Number Generators (McCulloch, Aug96)

   Written and submitted by J. Huston McCulloch for public,
   non-commercial use.
   Contact author at mcculloch.2@osu.edu; 614-292-0382.  */

/* This file contains two stable random number generator subroutines
     RNDSSTA -- generates vector of standard symmetric stable random
               numbers with alpha in [.1,2] and beta = 0.
     RNDSTA -- generates vector of standard stable random numbers with
               alpha in [.1,2] and arbitrary beta in [-1,1].
   At the head of the file is a demonstration driver program that
    calls RNDSTA, lets you type in alpha and beta, and makes nice
    graphs and histograms of the results.  Values of alpha in [1.5, 1.9],
    with beta in [-.3, .3], simulate many financial asset returns.
    Values of alpha in [1, 1.2] with beta = -1 often simulate icicles.
   See the comments at the beginning of RNDSTA and RNDSSTA for details
    of these procs.
   See McCulloch "Financial Applications of Stable Distributions",
    forthcoming 1996 in Volume 14 of the _Handbook of Statistics_ for
    a survey of pertinent literature.
*/

/* Driver program to demonstrate rndsta */
new;
library pgraph;
graphset;
_pdate = "";
_pltype = 6;
_pgrid = {1,1};
r = 1000;
newdraw:
print "Enter a value of alpha in [.1,2] ";
print "(Type Ctrl + c to exit.)";
alpha = con(1,1);
print "Enter a value of beta in [-1,1] ";
beta = con(1,1);
x = seqa(1,1,r);
y = rndsta(alpha,beta,r);
ylabel("x");
xy(x,y);
yy = sortc(y,1);
ymedian = yy[r/2];   @ (roughly) @
dy = 9.5;
yy = maxc(yy'|(ymedian-dy)*ones(1,rows(y)));
yy = minc(yy'|(ymedian+dy)*ones(1,rows(y)));
xlabel("Outermost bins include tails, empty bins appear slightly nonempty");
ylabel("frequency");
{breaks,midpoint,count} = hist(yy,55);
xlabel("");
goto newdraw;
end;


proc rndssta(alpha, r);
/* Returns rX1 vector of iid standard symmetric stable pseudo-random
   numbers with characteristic exponent alpha in [.1,2], and skewness
   parameter beta = 0, using method of Chambers, Mallows and Stuck
   (JASA 1976, 340-4).
   Encoded in GAUSS by J. Huston McCulloch, Ohio State University Econ Dept.
   (mcculloch.2@osu.edu), 12/95, directly from the CMS equation (4.1),
   as simplified when beta = 0.
   Each r.v. has log characteristic function
        log E exp(ixt) = -abs(t)^alpha
   When alpha = 2, the distribution is Gaussian, with variance 2.
   When alpha = 1, the distribution is standard Cauchy.
   If alpha > 1, Ex = 0.  In all cases, the median is 0.
   This proc uses 2r uniform pseudo-random numbers from RNDU.  */

local phi, w;
if (alpha < .1) or (alpha > 2);
  format /rdn 10,4;
  print "Alpha must be in [.1,2] for proc RNDSSTA." ;
  print "Actual value is " alpha;
  if alpha < .1 and alpha > 0;
    print "If alpha < .1, probability of overflow becomes non-negligible.";
    endif;
  stop;
  endif;
w = - ln(rndu(r,1));
phi = (rndu(r,1) - .5) * pi;
if feq(alpha,2);
  retp(2 * sqrt(w) .* sin(phi));
elseif feq(alpha,1);
  retp(tan(phi));
endif;
retp((cos((1-alpha) * phi) ./ w) ^ (1/alpha - 1) .* sin(alpha * phi)
             ./ cos(phi) ^ (1/alpha));
endp;


proc rndsta(alpha, beta, r);
/* Returns rX1 vector x of iid standard stable pseudo-random
   numbers with characteristic exponent alpha in [.1,2], and
   skewness parameter beta in [-1,1].  Based on the method of Chambers,
   Mallows and Stuck (JASA 1976, 340-4).  Encoded in GAUSS by
   J. Huston McCulloch, Ohio State Univ. Econ. Dept.
   (mcculloch.2@osu.edu), 12/95.
     The CMS method was applied in such a way that x will have the
   standard (delta =0, c = 1) characteristic function
    log E exp(ixt) = -abs(t)^alpha*(1-i*beta*sign(t)*tan(pi*alpha/2)),
                               for alpha /= 1,
                   = -abs(t)*(1+i*beta*(2/pi)*sign(t)*log(abs(t))),
                               for alpha = 1.
   With this parameterization, Ex = delta = 0 when alpha > 1, so that
   if the distribution is positively skewed (beta > 0),
   med(x) must become very negative as alpha approaches 1 from above.
   For alpha < 1, delta is the "focus of stability", so that if beta > 0,
   med(x) must become very positive as alpha approaches 1 from below.
   If beta /= 0, there is therefore a discontinuity in the distribution
   as a function of alpha as alpha passes 1.
     CMS remove this discontinuity by subtracting
                zeta = beta*tan(pi*alpha/2)
   (equivalent to their -tan(alpha*phi0)) from x for alpha /= 1
   in their program RSTAB (also known as GGSTA in IMSL).
   The result is a random number whose
   distribution is a continuous function of alpha, but whose location
   parameter has no known useful interpretation other than computational
   convenience.  The present program restores the standard parameterization
   by using the CMS (4.1), but with zeta added back in (ie with their
   initial tan(alpha*phi0) deleted).
     When alpha is precisely unity, x is computed from the CMS (2.4).
   Rather than using the CMS D2 and exp2 functions to compensate for ill-
   conditioning of (4.1) when alpha is very near 1, the present program
   merely fudges these cases by computing x from (2.4) and adjusting for
   zeta when alpha is within 1.e-8 of 1.  This should make no difference
   for simulation results with samples of size less than approximately
   10^8, and then only when the desired alpha is within 1.e-8 of 1 but
   not equal to 1.
     Companion proc RNDSSTA computes symmetric stable random variables
   with beta = 0. It is faster, and the above issues do not arise in
   these cases.
     When alpha = 2, the distribution is Gaussian, with variance 2, regardless
   of beta.  When alpha = 1 and beta = 0, the distribution is standard Cauchy.
   When alpha = .5 and beta = 1, the distribution is Levy (reciprocal
   Chi-squared(1)). */

local phi, w, zeta, aphi, a1phi, x, cosphi, bphi;
if alpha < .1 or alpha > 2;
  format /rdn 10,4;
  print "Alpha must be in [.1,2] for proc RNDSTA.";
  print "Actual value is " alpha;
  if alpha < .1 and alpha > 0;
    print "If alpha <.1, the probability of overflow output is non-negligible.";
    endif;
  stop;
  endif;
if abs(beta) > 1;
  format /rdn 10,4;
  print "Beta must be in [-1,1] for proc RNDSTA.";
  print "Actual value is " beta;
  stop;
  endif;
w = -ln(rndu(r,1));
phi = (rndu(r,1)-.5)*pi;
cosphi = cos(phi);
if abs(alpha-1) > 1.e-8;
  zeta = beta * tan(pi*alpha/2);
  aphi = alpha * phi;
  a1phi = (1 - alpha) * phi;
  x = ((sin(aphi) + zeta * cos(aphi)) ./ cosphi) .*
    ((cos(a1phi) + zeta * sin(a1phi)) ./ (w .* cosphi))^((1-alpha)/alpha);
  retp(x);
  endif;
bphi = (pi/2) + beta * phi;
x = (2/pi) * (bphi .* tan(phi) - beta * ln((pi/2) * w .* cosphi ./ bphi));
if alpha == 1;
  retp(x);
  else;
  zeta = beta * tan(pi * alpha/2);
  retp(x + zeta);
  endif;
endp;

