/*
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;