/* Stable Distribution Parameter Estimation by Quantile Method
   (McCulloch 96)
   Written and submitted by J. Huston McCulloch for public, non-
   commercial use, August 16, 1996.
   Contact author at mcculloch.2@osu.edu; 614-292-0382.    */

proc(5) = stab(x,symm);
/* McCulloch 12/95 GAUSS version (converted from FORTRAN)
  TAKES A STAB AT STABLE DISTRIBUTION PARAMETERS USING J.H. MCCULLOCH,
  "SIMPLE CONSISTENT ESTIMATORS OF STABLE DISTRIBUTION PARAMETERS."
  (COMMUNICATIONS IN STATISTICS:  SIMULATION & COMPUTATION 15#4,
  1109-1136 (1986))
  INPUT VARIABLES FOR STAB:
    x is a column vector of iid stable random variables with characteristic
      exponent alpha in [.5, 2].
    symm = 0 if beta is to be estimated.
           1 if only the symmetric case (beta = 0) is to be considered.
  OUTPUT VARIABLES:
     alpha = characteristic exponent.
     beta = skewness parameter (= 0 if symmetrical)
     c = scale parameter
     delta = location parameter
     zeta = shifted location parameter (= delta if beta = 0)
          = delta + beta*c*tan(pi*alpha/2), for alpha NE 1,
          = delta, for alpha = 1;
  PARAMETERIZATION:
     The distribution of x has the log characteristic function
        ln Eexp(i*x*t) = i*delta*t + psi(c*t),
     where the standard (c = 1, delta = 0) log c.f. is
        psi(t) = -abs(t)^alpha*[1-i*b*sign(t)*tan(pi*alpha/2)],
                    for alpha NE 1,
               = -abs(t)*[1+i*beta*(2/pi)*sign(t)*ln(abs(t)),
                    for alpha = 1.
     (See McCulloch, "On the Parameterization of the Afocal Stable
     Distribution," _Bulletin of the London Mathematical Society_,
     forthcoming 1996.  (1.1) in McCulloch (1986) should have a c before
     the t inside the log in the last line for internal consistency.)
     With this parameterization, if X has parameters (alpha, beta,
     c, delta), then (X-delta)/c has parameters (alpha, beta, 1, 0),
     for all alpha.

   Alpha is constrained to lie between 0.5 and 2.0.
     If the input data is inconsistent with any alpha value in this
     range, alpha will be returned either as 0.5 or as 2.0.  The former
     should be interpreted as meaning that alpha is indeterminately less
     than 0.5.  The latter indicates that the
     data is either Gaussian or too platykurtic (that's the
     opposite of leptokurtic) to be consistent with any stable distribution.
   Beta is constrained to lie between -1.0 and +1.0.  Note that if alpha
     is 2.0, beta has no effect on the distribution.  In this case 0.0
     is returned for beta regardless of the data.
   C is constrained to be positive.  If the interquartile range is
     0, the apparent scale is 0 and (-1,-1,-1,-1,-1) is returned.
     This is an error condition.                  */

local ena, enb, enc, za, ah, bh, q05, q25, q50, q75, q95, xx;
local d, cn, cnu, an, bn, sign, alpha, beta, c, delta, zeta;
local i, j, ii, jj, aa, b, bb, t1, t2, s1, s2, dt, ds, s, t;

@  ENA IS THE INDEX NU SUB ALPHA DESCRIBED IN MCCULLOCH (1986).  @
let ena[16,5] = {
           2.4388     2.4388     2.4388     2.4388     2.4388,
           2.5120     2.5117     2.5125     2.5129     2.5148,
           2.6080     2.6093     2.6101     2.6131     2.6174,
           2.7369     2.7376     2.7387     2.7420     2.7464,
           2.9115     2.9090     2.9037     2.8998     2.9016,
           3.1480     3.1363     3.1119     3.0919     3.0888,
           3.4635     3.4361     3.3778     3.3306     3.3161,
           3.8824     3.8337     3.7199     3.6257     3.5997,
           4.4468     4.3651     4.1713     4.0052     3.9635,
           5.2172     5.0840     4.7778     4.5122     4.4506,
           6.3140     6.0978     5.6241     5.2195     5.1256,
           7.9098     7.5900     6.8606     6.2598     6.1239,
          10.4480     9.9336     8.7790     7.9005     7.6874,
          14.8378    13.9540    12.0419    10.7219    10.3704,
          23.4831    21.7682    18.3320    16.2163    15.5841,
          44.2813    40.1367    33.0018    29.1399    27.7822 };
@  ENB IS THE INDEX NU SUB BETA. @
let enb[16,5] = {
           0.0000     0.0000     0.0000     0.0000     0.0000,
           0.0000     0.0179     0.0357     0.0533     0.0710,
           0.0000     0.0389     0.0765     0.1133     0.1480,
           0.0000     0.0626     0.1226     0.1784     0.2281,
           0.0000     0.0895     0.1736     0.2478     0.3090,
           0.0000     0.1183     0.2282     0.3199     0.3895,
           0.0000     0.1478     0.2849     0.3942     0.4686,
           0.0000     0.1769     0.3422     0.4703     0.5458,
           0.0000     0.2062     0.3993     0.5473     0.6210,
           0.0000     0.2362     0.4561     0.6240     0.6934,
           0.0000     0.2681     0.5134     0.6993     0.7616,
           0.0000     0.3026     0.5726     0.7700     0.8248,
           0.0000     0.3415     0.6343     0.8339     0.8805,
           0.0000     0.3865     0.6994     0.8900     0.9269,
           0.0000     0.4408     0.7678     0.9362     0.9620,
           0.0000     0.5095     0.8381     0.9700     0.9847 };
@  ENC IS THE INDEX NU SUB C  @
let enc[16,5] = {
           1.9078     1.9078     1.9078     1.9078     1.9078,
           1.9140     1.9150     1.9160     1.9185     1.9210,
           1.9210     1.9220     1.9275     1.9360     1.9470,
           1.9270     1.9305     1.9425     1.9610     1.9870,
           1.9330     1.9405     1.9620     1.9970     2.0430,
           1.9390     1.9520     1.9885     2.0450     2.1160,
           1.9460     1.9665     2.0220     2.1065     2.2110,
           1.9550     1.9845     2.0670     2.1880     2.3330,
           1.9650     2.0075     2.1255     2.2945     2.4910,
           1.9800     2.0405     2.2050     2.4345     2.6965,
           2.0000     2.0850     2.3115     2.6240     2.9735,
           2.0400     2.1490     2.4610     2.8865     3.3565,
           2.0980     2.2445     2.6765     3.2650     3.9125,
           2.1890     2.3920     3.0040     3.8440     4.7755,
           2.3370     2.6355     3.5425     4.8085     6.2465,
           2.5880     3.0735     4.5340     6.6365     9.1440 };
@  ZA IS THE INDEX NU SUB ZETA.
   DELTA IS ESTIMATED FROM ZETA, SO NU SUB DELTA IS NOT USED. @
let za[16,5] = {
           0.0000     0.0000     0.0000     0.0000     0.0000,
           0.0000    -0.0166    -0.0322    -0.0488    -0.0644,
           0.0000    -0.0302    -0.0615    -0.0917    -0.1229,
           0.0000    -0.0434    -0.0878    -0.1321    -0.1785,
           0.0000    -0.0556    -0.1113    -0.1699    -0.2315,
           0.0000    -0.0660    -0.1340    -0.2060    -0.2830,
           0.0000    -0.0751    -0.1542    -0.2413    -0.3354,
           0.0000    -0.0837    -0.1733    -0.2760    -0.3896,
           0.0000    -0.0904    -0.1919    -0.3103    -0.4467,
           0.0000    -0.0955    -0.2080    -0.3465    -0.5080,
           0.0000    -0.0980    -0.2230    -0.3830    -0.5760,
           0.0000    -0.0986    -0.2372    -0.4239    -0.6525,
           0.0000    -0.0956    -0.2502    -0.4688    -0.7424,
           0.0000    -0.0894    -0.2617    -0.5201    -0.8534,
           0.0000    -0.0779    -0.2718    -0.5807    -0.9966,
           0.0000    -0.0610    -0.2790    -0.6590    -1.1980 };
@ The above are used in transposed form ([16,5]) @
ena = ena';
enb = enb';
enc = enc';
za = za';

@  CN, AN, AND BN ARE THE SAMPLE VALUES OF THE INDICES NU SUB C,
   ALPHA, AND BETA, RESPECTIVELY.      @
xx = sortc(x,1);
q05 = frctl(xx,.05);
q25 = frctl(xx,.25);
q50 = frctl(xx,.50);
q75 = frctl(xx,.75);
q95 = frctl(xx,.95);
d = q95 - q05;
cn = q75 - q25;
if cn == 0;
  print "Warning -- interquartile range is 0";
  print "Press any key to continue"; wait;
  retp(-1,-1,-1,-1,-1);
endif;
an = d/cn;
if an < 2.4388;
  retp(2, 0, cn/1.9078, q50, q50);
endif;
if symm == 1; goto sstab; endif;

bn = (q05 + q95 - 2 * q50)/d;
@ The sign of beta is saved and restored at the end. @
sign = 1.0;
if bn < 0; sign = -1; endif;
if bn == 0; sign = 0; endif;
bn = abs(bn);
@   SUBSCRIPT J (OR JJ) RUNS FROM 1 TO 5 AS BETA RUNS FROM 0.0,
    0.25, ... 1.0.
    SUBSCRIPT I RUNS FROM 1 TO 16 AS ALPHA RUNS FROM 2.0, 1.9,
    ... 0.5.           @
ah = zeros(5,1);
j = 1; do until j > 5;
i = 2; do until i == 16;
if an <= ena[j,i]; break; endif;
i = i + 1; endo;
t = (an - ena[j,i-1])/(ena[j,i]-ena[j,i-1]);
ah[j] = 2 - (i-2+t)*.1;
j = j + 1; endo;

bh = zeros(16,1);
i = 2; do until i > 16;
j = 2; do until j == 5;
if bn <= enb[j,i]; break; endif;
j = j + 1; endo;
t = (bn - enb[j-1,i])/(enb[j,i] - enb[j-1,i]);
bh[i] = (j - 2 + t) * .25;
i = i + 1; endo;
bh[1] = 2 * bh[2] - bh[3];

j = 2; do until j > 5;
jj = j;
i = trunc((2 - ah[j]) * 10 + 2);
i = maxc(i|2);
i = minc(i|16);
aa = 2 - .1*(i-2);
s2 = -(ah[j] - aa) * 10;
b = (1 - s2) * bh[i-1] + s2 * bh[i];
if b < (j-1)/4; break; endif;
j = j + 1; endo;
j = jj;

bb = .25 * (j-2);
t1 = (bh[i-1] - bb) / .25;
t2 = (bh[i] - bb) / .25;
s1 = - (ah[j-1] - aa) * 10;
dt = t2 - t1;
ds = s2 - s1;
s = (s1 + t1 * ds) / (1 - ds * dt);
t = t1 + s * dt;
alpha = aa - s * .1;
if alpha < .5; alpha = .5; endif;
beta = bb + t * .25;
if beta > 1;
  beta = 1;
endif;
beta = beta * sign;
c = enc[j-1,i-1] * (1-s) * (1-t) + enc[j-1,i] * s * (1-t)
    + enc[j,i-1] * t * (1-s) + enc[j,i] * t * s;
c = cn/c;
zeta = za[j-1,i-1] * (1-s) * (1-t) + za[j-1,i] * s * (1-t)
       + za[j,i-1] * t * (1-s) + za[j,i] * t * s;
zeta = q50 + c * sign * zeta;
if alpha == 1;
  delta = zeta;
else;
  delta = zeta - beta * c * tan(pi*alpha / 2);
endif;
retp(alpha, beta, c, delta, zeta);

sstab:
i = 2; do until i == 16;
ii = i;
if an <= ena[1,i]; break; endif;
i = i + 1; endo;
i = ii;
t = (an - ena[1, i-1]) / (ena[1, i] - ena[1, i-1]);
alpha = (22 - i - t) / 10;
if alpha < .5; alpha = .5; endif;
cnu = enc[1, i-1] * (1-t) + enc[1,i] * t;
c = cn/cnu;
retp(alpha, 0, c, q50, q50);
endp;

proc frctl(x,p);
@ This proc obtains selected quantiles off the ordered vector x @
local zi, i, th;
zi = p*rows(x) + .5;
i = trunc(zi);
i = maxc(i|1);
i = minc(i|rows(x)-1);
th = zi - i;
retp((1-th)*x[i] + th * x[i+1]);
endp;
