/*
Title: rndKMsta for GAUSS 3.6+, rev. 4/24/02.
Author: J. Huston McCulloch (mcculloch.2@osu.edu)
Date:  May 7, 2001, revised April 24, 2002.
This code is provided gratis without any guarantees or warrantees.  
Proprietary modifications of this code are not permitted.
Please cite this code should it be used extensively in a research 
project.  
*/

/* This file contains the stable random number generator subroutine 
    rndKMsta for GAUSS 3.6+. 
   It is similar to the author's RNDSTA and RNDSSTA, but is based on  
    the GAUSS 3.6+ "Kiss Monster" rndKMu random number generator.
   Revised 4/24/02 to prevent 0 returns from rndKMu from creating infinite
    results.      
   At the head of the file is a demonstration driver program that
    calls rndKMsta, lets you type in alpha and beta, and makes nice
    graphs and histograms of the results.  Values of alpha in [1.6, 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 rndKMsta below for arguments, 
    returns, and details. 
   See J. Huston McCulloch "Financial Applications of Stable Distributions",
    in _Handbook of Statistics_, vol. 14, 1996, pp. 393-425, for a survey of
    pertinent literature.
*/

/* Driver program to demonstrate rndKMsta */
/* Driver revised 4/24/02 for GAUSS 3.6+ graphics. */ 
new;
library pgraph;
graphset;
_pdate = "";
_pltype = 6;
_pgrid = {1,1};

r = 1000;
state = hsec;  @ seeds rndKMu with clock on first pass @
x = 1; 
do while x > 0; @ deliberate infinite loop -- exit with alpha = 0 below. @
print "";
print "Enter a value of alpha in [.1,2] ";
print "(Enter 0 to exit.)";
alpha = con(1,1);
if alpha == 0; 
  print "";
  print "Exiting rndKMsta demo program."; 
  end;
  endif;
print "Enter a value of beta in [-1,1] ";
beta = con(1,1);
x = seqa(1,1,r);
{y, state} = rndKMsta(r, 1, alpha, beta, state);
xlabel("");
ylabel("x");
xy(x,y);
print "";
print "Graph may now be viewed on TKF viewer window.";
print "Type any key to compute histogram."
wait;
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);
print "";
print "";
print "Histogram may now be viewed on TKF viewer window.";
print "Type any key to try new parameter values.";
wait;
endo; 
end;
   
proc(2) = rndKMsta(r, c, alpha, beta, state);
/* Returns rXc matrix x of iid standard stable pseudo-random
    numbers with characteristic exponent alpha in [.1,2], and skewness
    parameter beta in {-1, 1], using method of Chambers, Mallows and Stuck
    (JASA 1976, 340-4), in conjunction with the GAUSS 3.6+ rndKMu 
    uniform random number generator.  Encoded in GAUSS by J. Huston McCulloch, 
    Ohio State University Econ Dept. (mcculloch.2@osu.edu), 4/01, 
    directly from the CMS equation (4.1).
    Revised 4/26/02 to eliminate 0 returns from rndKMu that may 
    result in infinite returns from this routine.  
     
   Outputs:  {x, state}   

   state is the state variable used by rndKMu.  
    On input on the first call, this may be 
    an integer seed in [0, 2^32-1].  Alternatively, setting seed = -1 
    lets the clock determine the seed.  On subsequent calls, it is 
    the 500x1 state vector returned by rndKMsta itself as the second output.   
    Each call uses up 2*r*c uniform pseudo-random numbers from rndKMu, 
    so the period is 2^(3859-1) = 2^3858.
 
   In the simple symmetric case beta = 0, each r.v. in x 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 symmetric cases, the median is 0.
    If you are only interested in the symmetric case, you may read 
    no further. 
   
   In the general case with beta /= 0, the CMS method 
    was applied in such a way that each x will have the
    standard (delta =0, c = 1) characteristic function
    log Eexp(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.
     In the symmetric case beta = 0, the above issues do not arise, and the 
    program automatically uses a faster and simpler version of the 
    general formula.  Separate procs for the symmetric and asymmetric 
    cases are not required.  
     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,state} = rndKMu(r,c,state);
w = w + 2^-33; 
@ adding 2^-33 changes range of rndKMu from [0,1) to (0,1). @  
w = -ln(w);
{phi,state} = rndKMu(r,c,state);
phi = (phi + (2^-33 -.5))*pi; 

@ Symmetric cases: @
 @ Gaussian case: @
  if alpha == 2;             
    retp(2 * sqrt(w) .* sin(phi), state);
   endif;
 @ Cauchy case: @ 
  if alpha == 1 and beta == 0;  
    retp(tan(phi), state);
   endif;
 @ Symmetric stable case: @
  if beta == 0;        
    retp((cos((1-alpha) * phi) ./ w) ^ (1/alpha - 1) .* sin(alpha * phi)
             ./ cos(phi) ^ (1/alpha), state);
   endif;

@ Asymmetric cases: @
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, state);
 endif;
bphi = (pi/2) + beta * phi;
x = (2/pi) * (bphi .* tan(phi) - beta * ln((pi/2) * w .* cosphi ./ bphi));
if alpha == 1;
  retp(x, state);
 else;                 @ fudges case abs(alpha-1) <= 1.e-8 @ 
  zeta = beta * tan(pi * alpha/2);
  retp(x + zeta, state);
 endif;
endp;