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