proc(2) = symstb(x,alpha);
@ June 1993: GAUSS version @
@ If you use this program, please cite @
@  J. Huston McCulloch, "Numerical Approximation of the Symmetric @
@  Stable Distribution and Density," OSU Econ Dept., Oct. 1994.  @
@ x may be a column vector; alpha is a scalar @
@ returns {sden,sdenp} @
@ sden = symmetric stable density @
@ sdenp = derivative of same wrt x @
@ Program symstb3 gives distribution fn as well @
@ Program symstb1 gives only density @
declare s_, sqpi_, a2_, cpxp0_, gpxp0_, cpxpp0_, gpxpp0_, cppp_, gppp_, c5i_;
declare a_, alf2_, alf1_, pd_, p_, pdd_;
declare first_ = 1;
declare oldalf_ = -999.99;
local r, q, znot, zn4, zn5;
local i, i1, j, j1;
local pialf, sp0, sppp0, xp0, xpp0, xppp0, spzp1;
local rp0, rpp0, rppp0, rp1, b, c;
local xx, xa1, xa1a, z, x1, x2, zp, x1p, x2p, rpz, sden;
local sign, zpp, x11, x1pp, x2a1, x2pp, cd, cdp, gd, gdp, rppz, cgr, sdenp;
fn afun(alpha) = 2^(1/alpha) - 1;
fn cden(x) = 1/(pi*(1+x.*x));
fn gden(x) = exp(maxc(-700*ones(1,rows(x)) | -(x.*x)'/4)) / (2*sqpi_);
fn zfun(x, alpha, a) = 1 - (1 + a*x) ^ (-alpha);
fn zpfun(x, alpha, a) = alpha * a * (1+a*x) ^ (-alpha -1);
fn xfun(z, alpha, a) = ((1-z)^(-1/alpha) - 1) / a;

if(not first_); goto starta; endif;
@ this s is the gij of the writeup @
s_ =    {1.8514190959D2,   -4.6769332663D2,
         4.8424720302D2,   -1.7639153404D2,
        -3.0236552164D2,    7.6351931975D2,
        -7.8560342101D2,    2.8426313374D2,
         4.4078923600D2,   -1.1181138121D3,
         1.1548311335D3,   -4.1969666223D2,
        -5.2448142165D2,    1.3224487717D3,
        -1.3555648053D3,    4.8834079950D2,
         5.3530435018D2,   -1.3374570340D3,
         1.3660140118D3,   -4.9286099583D2,
        -4.8988957866D2,    1.2091418165D3,
        -1.2285872257D3,    4.4063174114D2,
         3.2905528742D2,   -7.3211767697D2,
         6.8183641829D2,   -2.2824291084D2,
        -2.1495402244D2,    3.9694906604D2,
        -3.3695710692D2,    1.0905855709D2,
         2.1112581866D2,   -2.7921107017D2,
         1.1717966020D2,    3.4394664342D0,
        -2.6486798043D2,    1.1999093707D2,
         2.1044841328D2,   -1.5110881541D2,
         9.4105784123D2,   -1.7221988478D3,
         1.4087544698D3,   -4.2472511892D2,
        -2.1990475933D3,    4.2637720422D3,
        -3.4723981786D3,    1.0174373627D3,
         3.1047490290D3,   -5.4204210990D3,
         4.2221052925D3,   -1.2345971177D3,
        -5.1408260668D3,    1.1090264364D4,
        -1.0270337246D4,    3.4243449595D3,
         1.1215157876D4,   -2.4243529825D4,
         2.1536057267D4,   -6.8490996103D3,
        -1.8120631586D4,    3.1430132257D4,
        -2.4164285641D4,    6.9126862826D3,
         1.7388413126D4,   -2.2108397686D4,
         1.3397999271D4,   -3.1246611987D3,
        -7.2435775303D3,    4.3545399418D3,
         2.3616155949D2,   -7.6571653073D2,
        -8.7376725439D3,    1.5510852129D4,
        -1.3789764138D4,    4.6387417712D3};
s_ = reshape (s_,19, 4);
sqpi_ = sqrt(pi);
a2_ = afun(2);
cpxp0_ = 1/pi;
gpxp0_ = 1/(4*sqpi_*a2_);
cpxpp0_ = 2/pi;
gpxpp0_ = gpxp0_ * 1.5;
cppp_ = 4/pi;
gppp_ = gpxpp0_ * 2.5 - 1/(32*sqpi_*a2_^3);
c5i_ = {1, 5, 10, 10, 5, 1};
first_ = 0;

starta:
if alpha > 2.000001;
  format /rdn 20,15;
  print "alpha too big: " alpha;
  end;
  endif;
if alpha >= 2;
  sden = gden(x);
  sdenp = -sden .* x / 2;
  retp(sden, sdenp);
  endif;
if alpha == oldalf_; goto startx; endif;
znot = seqa(1,1,19) /20;
zn4 = (1-znot)^4;
zn5 = zn4 .* (1-znot);
a_ = afun(alpha);
alf2_ = 2 - alpha;
alf1_ = alpha - 1;
pialf = pi * alpha;
sp0 = dgamma(1/alpha) / pialf;
sppp0 = -dgamma(3/alpha) / pialf;
xp0 = 1/(alpha * a_);
xpp0 = xp0 * (1+alpha) / alpha;
xppp0 = xpp0 * (1+2*alpha)/alpha;
spzp1 = (a_^alpha) * dgamma(alpha) * sin(pialf / 2) / pi;
rp0 = -sp0 * xp0 + alf2_ * cpxp0_ + alf1_ * gpxp0_;
rpp0 = -sp0 * xpp0 + alf2_ * cpxpp0_ + alf1_ * gpxpp0_;
rppp0 = -sp0 * xppp0 - sppp0 * xp0^3 + alf2_ * cppp_ + alf1_ * gppp_;
rp1 = -spzp1 + alf2_ / pi;
q = zeros(6,1);
q[2] = rp0;
q[3] = rpp0/2;
q[4] = rppp0/6;
r = alf2_ * (s_ * (alf2_^seqa(1,1,4) - 1));
b = -sumc(q[2:4]) - sumc(r.*zn5);
c = rp1 - q[2:4]' * seqa(1,1,3) - 5 * sumc(r.*zn4);
q[5] = 5*b-c;
q[6] = c-4*b;
p_ = zeros(6,20);
i = 0;
do until i > 5;
i1 = i + 1;
j = 1;
do until j > 19;
j1 = j + 1;
p_[i1, j1] = p_[i1, j] + r[j] * (-znot[j])^(5-i);
j = j1;
endo;
i = i1;
endo;
p_ = q + c5i_ .* p_;
pd_ = p_[2:6,.] .* seqa (1,1,5);
pdd_ = pd_[2:5,.] .* seqa (1,1,4);
oldalf_ = alpha;

startx:
sign = 2 * (x .>= zeros(rows(x),1)) - 1;
xx = abs(x);
xa1 = 1 + a_ * xx;
xa1a = xa1 ^ (-alpha);
z = 1 - xa1a;
zp = alpha * a_ * xa1a ./ xa1;
zpp = -(alpha + 1) * a_ * zp ./ xa1;
x11 = 1/xa1a;
x1 = x11 - 1;
x1p = x11 .* x11;
x1pp = 2 * x11 .* x1p;
x2a1 = 1/sqrt(xa1a);
x2 = (x2a1 - 1) / a2_;
x2p = x2a1 ./ (2*a2_ * xa1a);
x2pp = 1.5 * x2p ./ xa1a;
cd = cden (x1);
cdp = - 2 * pi * x1 .* cd .* cd;
gd = gden(x2);
gdp = - x2 .* gd / 2;
j1 = trunc(20 * z) + 1;
rpz = poly(pd_[.,j1], z, 4);
rppz = poly(pdd_[.,j1], z, 3);
cgr = alf2_ * cd .* x1p + alf1_ * gd .* x2p - rpz;
sden = cgr .* zp;
sdenp = sign .* (cgr .* zpp + zp .* zp .* (alf2_ * (cd.*x1pp + x1p.*x1p.*cdp)
         + alf1_ * (gd .* x2pp + x2p .* x2p .* gdp) -rppz));
retp (sden, sdenp);
endp;

proc poly(a, x, k);
local sum, j;
sum = a[k+1,.]';
j = k;
do until j ==0;
sum = sum .* x + a[j,.]';
j = j - 1;
endo;
retp(sum);
endp;

proc dgamma(x);
@ Note -- GAUSS 3.2+ has adopted an improved gamma function similar
to this one, making this proc redundant. @
local t, ser, c, stp;
c = {76.18009172947146, -86.50532032941677, 24.01409824083091,
     -1.231739572450155, .001208650973866179, -.5395239384953d-5};
stp = 2.5066282746310005;
t = x + 5.5;
t = t^(x+.5)/exp(t);
ser = 1.000000000190015 + sumc(c./(x+seqa(1,1,6)));
retp(t * stp * ser / x);
endp;

proc(3) = symstb3(x,alpha);
@ June 1993: GAUSS version @
@ If you use this program, please cite @
@  J. Huston McCulloch, "Numerical Approximation of the Symmetric @
@  Stable Distribution and Density," OSU Econ Dept., Oct. 1994.  @
@ x may be a column vector; alpha is a scalar @
@ returns {sdisc,sden,sdenp} @
@ sdisc = complemented symmetric stable distribution @
@ sden = symmetric stable density @
@ sdenp = derivative of density wrt x @
declare s_, sqpi_, a2_, cpxp0_, gpxp0_, cpxpp0_, gpxpp0_, cppp_, gppp_, c5i_;
declare a_, alf2_, alf1_, pd_, p_, pdd_, sq2_;
declare first_ = 1;
declare oldalf_ = -999.99;
local r, q, znot, zn4, zn5;
local i, i1, j, j1;
local pialf, sp0, sppp0, xp0, xpp0, xppp0, spzp1;
local rp0, rpp0, rppp0, rp1, b, c;
local xx, xa1, xa1a, z, x1, x2, zp, x1p, x2p, rpz, sden, pos, sdisc, rz;
local sign, zpp, x11, x1pp, x2a1, x2pp, cd, cdp, gd, gdp, rppz, cgr, sdenp;
fn afun(alpha) = 2^(1/alpha) - 1;
fn cden(x) = 1/(pi*(1+x.*x));
fn gden(x) = exp(maxc(-700*ones(1,rows(x)) | -(x.*x)'/4)) / (2*sqpi_);
fn cfun(x) = .5 - atan(x) / pi;
fn gfun(x) = cdfnc (x/sq2_);
fn zfun(x, alpha, a) = 1 - (1 + a*x) ^ (-alpha);
fn zpfun(x, alpha, a) = alpha * a * (1+a*x) ^ (-alpha -1);
fn xfun(z, alpha, a) = ((1-z)^(-1/alpha) - 1) / a;

if(not first_); goto starta; endif;
@ this s is the gij of the writeup @
s_ =    {1.8514190959D2,   -4.6769332663D2,
         4.8424720302D2,   -1.7639153404D2,
        -3.0236552164D2,    7.6351931975D2,
        -7.8560342101D2,    2.8426313374D2,
         4.4078923600D2,   -1.1181138121D3,
         1.1548311335D3,   -4.1969666223D2,
        -5.2448142165D2,    1.3224487717D3,
        -1.3555648053D3,    4.8834079950D2,
         5.3530435018D2,   -1.3374570340D3,
         1.3660140118D3,   -4.9286099583D2,
        -4.8988957866D2,    1.2091418165D3,
        -1.2285872257D3,    4.4063174114D2,
         3.2905528742D2,   -7.3211767697D2,
         6.8183641829D2,   -2.2824291084D2,
        -2.1495402244D2,    3.9694906604D2,
        -3.3695710692D2,    1.0905855709D2,
         2.1112581866D2,   -2.7921107017D2,
         1.1717966020D2,    3.4394664342D0,
        -2.6486798043D2,    1.1999093707D2,
         2.1044841328D2,   -1.5110881541D2,
         9.4105784123D2,   -1.7221988478D3,
         1.4087544698D3,   -4.2472511892D2,
        -2.1990475933D3,    4.2637720422D3,
        -3.4723981786D3,    1.0174373627D3,
         3.1047490290D3,   -5.4204210990D3,
         4.2221052925D3,   -1.2345971177D3,
        -5.1408260668D3,    1.1090264364D4,
        -1.0270337246D4,    3.4243449595D3,
         1.1215157876D4,   -2.4243529825D4,
         2.1536057267D4,   -6.8490996103D3,
        -1.8120631586D4,    3.1430132257D4,
        -2.4164285641D4,    6.9126862826D3,
         1.7388413126D4,   -2.2108397686D4,
         1.3397999271D4,   -3.1246611987D3,
        -7.2435775303D3,    4.3545399418D3,
         2.3616155949D2,   -7.6571653073D2,
        -8.7376725439D3,    1.5510852129D4,
        -1.3789764138D4,    4.6387417712D3};
s_ = reshape (s_,19, 4);
sqpi_ = sqrt(pi);
sq2_ = sqrt(2);
a2_ = afun(2);
cpxp0_ = 1/pi;
gpxp0_ = 1/(4*sqpi_*a2_);
cpxpp0_ = 2/pi;
gpxpp0_ = gpxp0_ * 1.5;
cppp_ = 4/pi;
gppp_ = gpxpp0_ * 2.5 - 1/(32*sqpi_*a2_^3);
c5i_ = {1, 5, 10, 10, 5, 1};
first_ = 0;

starta:
if alpha > 2.000001;
  format /rdn 20,15;
  print "alpha too big: " alpha;
  end;
  endif;
if alpha >= 2;
  sdisc = gfun(x);
  sden = gden(x);
  sdenp = -sden .* x / 2;
  retp(sdisc, sden, sdenp);
  endif;
if alpha == oldalf_; goto startx; endif;
znot = seqa(1,1,19) /20;
zn4 = (1-znot)^4;
zn5 = zn4 .* (1-znot);
a_ = afun(alpha);
alf2_ = 2 - alpha;
alf1_ = alpha - 1;
pialf = pi * alpha;
sp0 = dgamma(1/alpha) / pialf;
sppp0 = -dgamma(3/alpha) / pialf;
xp0 = 1/(alpha * a_);
xpp0 = xp0 * (1+alpha) / alpha;
xppp0 = xpp0 * (1+2*alpha)/alpha;
spzp1 = (a_^alpha) * dgamma(alpha) * sin(pialf / 2) / pi;
rp0 = -sp0 * xp0 + alf2_ * cpxp0_ + alf1_ * gpxp0_;
rpp0 = -sp0 * xpp0 + alf2_ * cpxpp0_ + alf1_ * gpxpp0_;
rppp0 = -sp0 * xppp0 - sppp0 * xp0^3 + alf2_ * cppp_ + alf1_ * gppp_;
rp1 = -spzp1 + alf2_ / pi;
q = zeros(6,1);
q[2] = rp0;
q[3] = rpp0/2;
q[4] = rppp0/6;
r = alf2_ * (s_ * (alf2_^seqa(1,1,4) - 1));
b = -sumc(q[2:4]) - sumc(r.*zn5);
c = rp1 - q[2:4]' * seqa(1,1,3) - 5 * sumc(r.*zn4);
q[5] = 5*b-c;
q[6] = c-4*b;
p_ = zeros(6,20);
i = 0;
do until i > 5;
i1 = i + 1;
j = 1;
do until j > 19;
j1 = j + 1;
p_[i1, j1] = p_[i1, j] + r[j] * (-znot[j])^(5-i);
j = j1;
endo;
i = i1;
endo;
p_ = q + c5i_ .* p_;
pd_ = p_[2:6,.] .* seqa (1,1,5);
pdd_ = pd_[2:5,.] .* seqa (1,1,4);
oldalf_ = alpha;

startx:
pos = x .>= zeros(rows(x),1);
sign = 2 * pos - 1;
xx = abs(x);
xa1 = 1 + a_ * xx;
xa1a = xa1 ^ (-alpha);
z = 1 - xa1a;
zp = alpha * a_ * xa1a ./ xa1;
zpp = -(alpha + 1) * a_ * zp ./ xa1;
x11 = 1/xa1a;
x1 = x11 - 1;
x1p = x11 .* x11;
x1pp = 2 * x11 .* x1p;
x2a1 = 1/sqrt(xa1a);
x2 = (x2a1 - 1) / a2_;
x2p = x2a1 ./ (2*a2_ * xa1a);
x2pp = 1.5 * x2p ./ xa1a;
cd = cden (x1);
cdp = - 2 * pi * x1 .* cd .* cd;
gd = gden(x2);
gdp = - x2 .* gd / 2;
j1 = trunc(20 * z) + 1;
rz = poly(p_[.,j1], z, 5);
rpz = poly(pd_[.,j1], z, 4);
rppz = poly(pdd_[.,j1], z, 3);
cgr = alf2_ * cd .* x1p + alf1_ * gd .* x2p - rpz;
sden = cgr .* zp;
sdenp = sign .* (cgr .* zpp + zp .* zp .* (alf2_ * (cd.*x1pp + x1p.*x1p.*cdp)
         + alf1_ * (gd .* x2pp + x2p .* x2p .* gdp) -rppz));
sdisc = alf2_*cfun(x1) + alf1_*gfun(x2) + rz;
sdisc = sdisc .* pos + (1-sdisc) .* (1-pos);
retp (sdisc, sden, sdenp);
endp;

proc symstb1(x,alpha);
@ June 1993: GAUSS version @
@ If you use this program, please cite @
@  J. Huston McCulloch, "Numerical Approximation of the Symmetric @
@  Stable Distribution and Density," OSU Econ Dept., Oct. 1994.  @
@ x may be a column vector; alpha is a scalar @
@ returns sden = symmetric stable density @
@ Program symstb3 gives distribution fn as well @
@ Program symstb1 gives only density @
declare s_, sqpi_, a2_, cpxp0_, gpxp0_, cpxpp0_, gpxpp0_, cppp_, gppp_, c5i_;
declare a_, alf2_, alf1_, pd_, p_, pdd_;
declare first_ = 1;
declare oldalf_ = -999.99;
local r, q, znot, zn4, zn5;
local i, i1, j, j1;
local pialf, sp0, sppp0, xp0, xpp0, xppp0, spzp1;
local rp0, rpp0, rppp0, rp1, b, c;
local xa1, xa1a, z, x1, x2, zp, x1p, x2p, rpz, sden;
local x11, x2a1, x2pp, cd, gd;
fn afun(alpha) = 2^(1/alpha) - 1;
fn cden(x) = 1/(pi*(1+x.*x));
fn gden(x) = exp(maxc(-700*ones(1,rows(x)) | -(x.*x)'/4)) / (2*sqpi_);
fn zfun(x, alpha, a) = 1 - (1 + a*x) ^ (-alpha);
fn zpfun(x, alpha, a) = alpha * a * (1+a*x) ^ (-alpha -1);
fn xfun(z, alpha, a) = ((1-z)^(-1/alpha) - 1) / a;

if(not first_); goto starta; endif;
@ this s is the gij of the writeup @
s_ =    {1.8514190959D2,   -4.6769332663D2,
         4.8424720302D2,   -1.7639153404D2,
        -3.0236552164D2,    7.6351931975D2,
        -7.8560342101D2,    2.8426313374D2,
         4.4078923600D2,   -1.1181138121D3,
         1.1548311335D3,   -4.1969666223D2,
        -5.2448142165D2,    1.3224487717D3,
        -1.3555648053D3,    4.8834079950D2,
         5.3530435018D2,   -1.3374570340D3,
         1.3660140118D3,   -4.9286099583D2,
        -4.8988957866D2,    1.2091418165D3,
        -1.2285872257D3,    4.4063174114D2,
         3.2905528742D2,   -7.3211767697D2,
         6.8183641829D2,   -2.2824291084D2,
        -2.1495402244D2,    3.9694906604D2,
        -3.3695710692D2,    1.0905855709D2,
         2.1112581866D2,   -2.7921107017D2,
         1.1717966020D2,    3.4394664342D0,
        -2.6486798043D2,    1.1999093707D2,
         2.1044841328D2,   -1.5110881541D2,
         9.4105784123D2,   -1.7221988478D3,
         1.4087544698D3,   -4.2472511892D2,
        -2.1990475933D3,    4.2637720422D3,
        -3.4723981786D3,    1.0174373627D3,
         3.1047490290D3,   -5.4204210990D3,
         4.2221052925D3,   -1.2345971177D3,
        -5.1408260668D3,    1.1090264364D4,
        -1.0270337246D4,    3.4243449595D3,
         1.1215157876D4,   -2.4243529825D4,
         2.1536057267D4,   -6.8490996103D3,
        -1.8120631586D4,    3.1430132257D4,
        -2.4164285641D4,    6.9126862826D3,
         1.7388413126D4,   -2.2108397686D4,
         1.3397999271D4,   -3.1246611987D3,
        -7.2435775303D3,    4.3545399418D3,
         2.3616155949D2,   -7.6571653073D2,
        -8.7376725439D3,    1.5510852129D4,
        -1.3789764138D4,    4.6387417712D3};
s_ = reshape (s_,19, 4);
sqpi_ = sqrt(pi);
a2_ = afun(2);
cpxp0_ = 1/pi;
gpxp0_ = 1/(4*sqpi_*a2_);
cpxpp0_ = 2/pi;
gpxpp0_ = gpxp0_ * 1.5;
cppp_ = 4/pi;
gppp_ = gpxpp0_ * 2.5 - 1/(32*sqpi_*a2_^3);
c5i_ = {1, 5, 10, 10, 5, 1};
first_ = 0;

starta:
if alpha > 2.000001;
  format /rdn 20,15;
  print "alpha too big: " alpha;
  end;
  endif;
if alpha >= 2;
  sden = gden(x);
  retp(sden);
  endif;
if alpha == oldalf_; goto startx; endif;
znot = seqa(1,1,19) /20;
zn4 = (1-znot)^4;
zn5 = zn4 .* (1-znot);
a_ = afun(alpha);
alf2_ = 2 - alpha;
alf1_ = alpha - 1;
pialf = pi * alpha;
sp0 = dgamma(1/alpha) / pialf;
sppp0 = -dgamma(3/alpha) / pialf;
xp0 = 1/(alpha * a_);
xpp0 = xp0 * (1+alpha) / alpha;
xppp0 = xpp0 * (1+2*alpha)/alpha;
rp0 = -sp0 * xp0 + alf2_ * cpxp0_ + alf1_ * gpxp0_;
rpp0 = -sp0 * xpp0 + alf2_ * cpxpp0_ + alf1_ * gpxpp0_;
rppp0 = -sp0 * xppp0 - sppp0 * xp0^3 + alf2_ * cppp_ + alf1_ * gppp_;
spzp1 = (a_^alpha) * dgamma(alpha) * sin(pialf / 2) / pi;
rp1 = -spzp1 + alf2_ / pi;
q = zeros(6,1);
q[2] = rp0;
q[3] = rpp0/2;
q[4] = rppp0/6;
r = alf2_ * (s_ * (alf2_^seqa(1,1,4) - 1));
b = -sumc(q[2:4]) - sumc(r.*zn5);
c = rp1 - q[2:4]' * seqa(1,1,3) - 5 * sumc(r.*zn4);
q[5] = 5*b-c;
q[6] = c-4*b;
p_ = zeros(6,20);
i = 0;
do until i > 5;
i1 = i + 1;
j = 1;
do until j > 19;
j1 = j + 1;
p_[i1, j1] = p_[i1, j] + r[j] * (-znot[j])^(5-i);
j = j1;
endo;
i = i1;
endo;
p_ = q + c5i_ .* p_;
pd_ = p_[2:6,.] .* seqa (1,1,5);
oldalf_ = alpha;

startx:
xa1 = 1 + a_ * abs(x);
xa1a = xa1 ^ (-alpha);
z = 1 - xa1a;
zp = alpha * a_ * xa1a ./ xa1;
x11 = 1/xa1a;
x1 = x11 - 1;
x1p = x11 .* x11;
x2a1 = 1/sqrt(xa1a);
x2 = (x2a1 -1) /a2_;
x2p = x2a1 ./ (2*a2_ * xa1a);
cd = cden (x1);
gd = gden(x2);
j1 = trunc(20 * z) + 1;
rpz = poly(pd_[.,j1], z, 4);
sden = (alf2_ * cd .* x1p + alf1_ * gd .* x2p - rpz) .* zp;
retp (sden);
endp;