@ 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;