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