@ file stabmed @
@ This file contains procs rndsta and rblob2 for web posting @
@ driver program to illustrate proc stabmed @
new;
mean = 0.001;
stdev = 0.008;
skew = -0.369;
kurt = 4.983;
nobs = 100;
{alpha, beta, c, delta} = stabmed(mean, stdev, skew, kurt, nobs);
format /rdn 10,6;
print "alpha = " alpha;
print "beta = " beta;
print "c = " c;
print "delta = " delta;
end;
proc(4) = stabmed(mean, stdev, skew, kurt, nobs);
/* Finds stable distribution whose median sample skewness, kurtosis, standard
deviation and mean match given values, by Monte Carlo simulation.
Written by J. Huston McCulloch, Econ. Dept., Ohio State Univ., 7/00,
mcculloch.2@osu.edu.
Inputs:
mean = sample mean.
stdev = sample standard deviation, based on unbiased variance.
skew = sample skewness, equal to third moment about mean divided
by 1.5 power of second moment about mean. Both moments are divided
by nobs, rather than nobs-1 as in computing stdev.
kurt = sample kurtosis, equal to fourth moment about mean divided
by square of second moment about mean, again dividing by nobs.
Do not use "excess kurtosis", from which 3 has been subtracted.
nobs = number of observations in sample.
Outputs:
alpha = characteristic exponent, theoretically restricted to (0, 2],
but here restricted to [.1, 2].
beta = skewness parameter, restricted to [-1, 1]. In the Gaussian
case alpha = 2, this parameter is unidentified, and the returned
value may be disregarded.
c = scale parameter (= std. dev./sqrt(2) in Gaussian case alpha = 2).
delta = location parameter (= population mean when alpha > 1).
If no solution is found, the vector {-888.88, -888.88, -888.88, -888.88}
is returned. This will typically occur for kurt < 3, or when
skewness is excessive relative to kurtosis.
Parameterization is as in McCulloch, "On the Parametrization of the
Afocal Stable Distributions," Bulletin of the London Mathematical
Society, vol. 28 (1996), pp. 651-55.
Delta and c are estimated using the medians of the standardized stable
distributions having the estimated alpha and beta. Consequently
delta will not exactly match the given mean, even when alpha > 1 so
that the population mean exists, and c*sqrt(2) will not exactly match
the given stdev even in the Gaussian case alpha = 2.
This program is a very noisy way to estimate stable distribution parameters,
and is offered only to give users a sense that sample moments are
not entirely without information about the parameters of a distribution
whose corresponding population moments do not exist. Estimated values
are quite sensitive to nobs.
For efficient estimation, use ML, eg with John Nolan's program STABLE.EXE,
at www.cas.american.edu/~jpnolan . Nolan's "S" parameterization
corresponds to that used here (with his "gamma" in place of my "c"),
except for the location parameter in the "afocal" cases alpha = 1,
beta =/ 0. (see above paper for details).
Monte Carlo estimates are based on 1001 replications of size nobs.
Change nreps in procs stabmed and lossfun if another value is desired.
The program
uses a Nelder-Mead polytope (aka simplex) algorithm to minimize the
sum of squared deviations between the given skew and kurt and the
median Monte Carlo value. Estimation may take several minutes with
nobs = 1000. Progress of the
Nelder-Mead algorithm is monitored on-screen, so that the program will
not appear to be brain dead.
To make the Monte Carlo loss function a smooth function of the stable
parameters, the same seed is used for all parameter values. A different
seed would give slightly different results.
*/
local theta0, dtheta, theta, loss, its;
local m, v, m1, x, xx, rep, nreps, n2, c, delta;
@ kurt, skew, and nobs are globalized for passing to lossfunc @
declare kurt_, skew_, nobs_;
kurt_ = kurt;
skew_ = skew;
nobs_ = nobs;
/* Theta is parameter vector {alpha, beta}.
Initial values of theta for rblob2 implied by theta0 and dtheta are
{1.5, 0}, {2, 0} and {1.5, .5}. */
theta0 = {1.5 0.0};
dtheta = {.5 .5};
/* Alpha (theta[1]) is "smartly" restricted to [.1, 2] by rblob2, using
its last 2 arguments.
Beta (theta[2]) is indirectly restricted to [-1, 1] by having lossfunc
return an impossibly bad loss value. */
{theta, loss, its} = rblob2(theta0, dtheta, &lossfunc, .0001, 1, .1, 2);
if loss < -.001;
print "No stable distribution match found by STABMED.";
retp(-888.88, -888.88, -888.88, -888.88);
endif;
/* Compute c and delta using estimated alpha and beta: */
nreps = 1001;
n2 = (nreps-1)/2 + 1;
rndseed 123456789.;
x = rndu(100, 1);
m = zeros(nreps,1);
v = zeros(nreps,1);
rep = 1;
do until rep > nreps;
x = rndsta(theta[1], theta[2], nobs_);
m1 = sumc(x)/nobs_;
xx = x - m1;
m[rep] = m1;
v[rep] = sumc(xx.^2)/(nobs_-1);
rep = rep+1;
endo;
m = sortc(m, 1);
v = sortc(v, 1);
c = stdev/sqrt(v[n2]);
delta = mean - m[n2]*c;
retp(theta[1], theta[2], c, delta);
endp;
proc lossfunc(theta);
local nreps, n2, kurt, skew, rep, x, xx, mean, m2, m3, m4;
if theta[1] > 2 or theta[1] < .1 or abs(theta[2]) > 1;
retp(-1.e100);
endif;
nreps = 1001;
n2 = (nreps-1)/2 + 1;
rndseed 123456789.;
x = rndu(100, 1);
kurt = zeros(nreps,1);
skew = zeros(nreps,1);
rep = 1;
do until rep > nreps;
x = rndsta(theta[1], theta[2], nobs_);
mean = sumc(x)/nobs_;
xx = x - mean;
m2 = sumc(xx.^2)/nobs_;
m3 = sumc(xx.^3)/nobs_;
m4 = sumc(xx.^4)/nobs_;
kurt[rep] = m4/(m2^2);
skew[rep] = m3/(m2^1.5);
rep = rep+1;
endo;
kurt = sortc(kurt,1);
skew = sortc(skew,1);
retp(- (kurt[n2]-kurt_)^2 - (skew[n2]-skew_)^2);
endp;
proc(3) = rblob2(x0, dx0, &func, tol, view, x1min, x1max);
@ Restricted uphill simplex search @
@ Finds 1xn xhi that maximizes func(x), by uphill simplex method.
Returns xhi, yhi, its.
View = 1 for on-screen progress, 0 for none, 2 for each x with pause.
Argument x0 is 1xn vector of initial guesses.
dx0 is 1xn vector of perturbations to x0.
func expects a row vector
converges when dx < tol * abs(dx0) @
@ x[1] is constrained to lie between x1min and x1max @
local func:proc;
local n, its, i1, y, i, j, ylo, yhi, centroid;
local x, dx, xnew, ynew, a,maxgrad;
format /rdn 10,3;
n = cols(x0);
if cols (dx0) /= n; print "blob: x0 " n " dx0 " cols(dx0); end; endif;
x = x0|(ones(n,1)*x0+diagrv(zeros(n,n),dx0));
its = 0;
i1 = 1;
y = zeros(n+1,1);
start:
i = i1;
do until i > n+1;
y[i] = func(x[i,.]);
its = its + 1;
i = i + 1;
endo;
rank:
y = rev(sortc(y~x,1));
x = y[.,2:n+1];
y = y[.,1];
ylo = y[n+1];
yhi = y[1];
if view == 2; print x~y; wait; endif;
@ test@;
@
maxgrad = inv(x[2:n+1,.] - x[1,.]) * (y[2:n+1] - y[1]);
if feq(x[1,1], x1min) or feq(x[1,1],x1max); maxgrad[1] = 0; endif;
maxgrad = maxc(abs(maxgrad)); @
if maxc((maxc(x) - minc(x))./abs(dx0')) < tol; @ AND maxgrad < tol; @
@ if maxc(maxc(x) - minc(x)) < tol; AND maxgrad < tol; @
retp(x[1,.], yhi, its);
endif;
@ begin reflection @
centroid = sumc(x[1:n,.])'/n;
dx = centroid - x[n+1,.];
@ if centroid is already on boundary, skip to halve @
if feq(centroid[1],x1min) or feq(centroid[1],x1max); goto halve; endif;
a = 1;
xnew = centroid + dx;
if xnew[1] > x1max;
a = (x1max - centroid[1]) / dx[1];
elseif xnew[1] < x1min;
a = (x1min - centroid[1]) / dx[1];
endif;
xnew = centroid + a * dx;
ynew = func(xnew);
if ynew <= ylo; goto halve; endif;
x[n+1,.] = xnew;
y[n+1] = ynew;
its = its + 1;
if view; print "reflect " its xnew ynew yhi; endif;
if ynew <= yhi; goto rank; endif;
@ try doubling @
if a < 1; goto rank; endif;
a = 2;
xnew = centroid + 2 * dx;
if xnew[1] > x1max;
a = (x1max - centroid[1]) / dx[1];
elseif xnew[1] < x1min;
a = (x1min - centroid[1]) / dx[1];
endif;
xnew = centroid + a * dx;
ynew = func(xnew);
if ynew <= y[n+1]; goto rank; endif;
x[n+1,.] = xnew;
y[n+1] = ynew;
its = its + 1;
if view; print "double " its xnew ynew yhi; endif;
goto rank;
halve: ;
xnew = centroid -.5 * dx;
ynew = func(xnew);
if ynew <= ylo; goto contract; endif;
x[n+1,.] = xnew;
y[n+1] = ynew;
its = its + 1;
if view; print "halve " its xnew ynew yhi; endif;
goto rank;
contract:
x = .5 * (x[1,.] + x);
i1 = 2;
if view; print "contract" its x[1,.] yhi yhi; endif;
goto start;
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;