/* 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 demonstrates 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).
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 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.
Companion proc RNDSSTA computes symmetric stable random variables
for 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;