// 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. // // Translated from GAUSS to SCILAB by jcosmao@monaco.edu - June 2007 function [retp] = stabfit(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. // Warning: The function may not work if x does not contain at least 100 datas. // This is due to the 'perctl' function. // ENA IS THE INDEX NU SUB ALPHA DESCRIBED IN MCCULLOCH (1986). ena = [ 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. enb = [ 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 enc = [ 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. za = [ 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'; pi=3.141592654; // CN, AN, AND BN ARE THE SAMPLE VALUES OF THE INDICES NU SUB C, // ALPHA, AND BETA, RESPECTIVELY. xx = sort(x); q = perctl(xx,[05 25 50 65 95]); q05 = q(1,1); q25 = q(2,1); q50 = q(3,1); q75 = q(4,1); q95 = q(5,1); d = q95 - q05; cn = q75 - q25; if cn == 0 warning('Warning - interquartile range is 0') retp = [-1,-1,-1,-1,-1]; return end an = d/cn; if an < 2.4388 retp= [2, 0, cn/1.9078, q50, q50]; return end if symm == 1 for i = 2:15 if an <= ena(1,i); break; end end t = (an - ena(1, i-1)) / (ena(1, i) - ena(1, i-1)); alpha = (22 - i - t) / 10; if alpha < .5; alpha = .5; end cnu = enc(1, i-1) * (1-t) + enc(1,i) * t; c = cn/cnu; retp=[alpha, 0, c, q50, q50]; return end 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; end if bn == 0; sign = 0; end 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); for j = 1:5 for i = 2:15 if an <= ena(j,i); break; end end t = (an - ena(j,i-1))/(ena(j,i)-ena(j,i-1)); ah(j) = 2 - (i-2+t)*.1; end bh = zeros(16,1); for i = 2:16 for j = 2:4 if bn <= enb(j,i); break; end end t = (bn - enb(j-1,i))/(enb(j,i) - enb(j-1,i)); bh(i) = (j - 2 + t) * .25; end bh(1) = 2 * bh(2) - bh(3); for j = 2:5 i = floor((2 - ah(j)) * 10 + 2); i = max([i 2]); i = min([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; end end 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; end beta = bb + t * .25; if beta > 1 beta = 1; end 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); end retp=[alpha, beta, c, delta, zeta]; endfunction