/* This file contains proc SMSTRG, which performs a
linear regression with symmetric stable disturbances,
preceded by a short driver program to demonstrate its
use.  See the initial comments in SMSTRG for details.

Included in this file are all procs called by SMSTRG,
namely SRLIKC, SYMSTB1, POLY, DGAMMA, STAB, FRCTL and MYOLS.
Also included is RNDSSTA, which is used by the demonstration program.
Numerous globals are generated, all of which are "uglified"
by the addtion of a terminal underscore (_) to reduce
the likelihood of a conflict. */

/* demonstrates smstrg */
new;
n = 1000;
alpha = 1.5;
x = seqa(0,1,n)/n;
y = 10 + 3 * x + 2 * rndssta(alpha, n);
{alphahat, chat, betahat, cov, its, maxlik, lrstat} = smstrg(y,x,1,.0001,2,1);
end;

proc(7) = smstrg(y,x,const,tol, vu, printout);
/* Performs linear regression y = X * beta + epsilon with
   symmetric stable residuals by quasi-Newton method, using
   OPG matrix.
   Requires procs symstb1, myols, stab, srlikc.
   Returns {alpha, c, beta, cov, its, maxlik, lrstat}.
   Converges when d log L (adjusted for splits) < tol .
   const = 1 to add a constant, which will be beta[1]
   vu = 0 for nothing on screen
      = 1 for results only on screen
      = 2 for progress on screen
      = 3 for above with pauses
   printout = 0 for nothing
            = 1 for printed results
   cov = cov(alpha, ln(c), beta'), as computed from numerical Hessian.
   If alpha is at (or within numerical increment of) .84
   or 2.00, alpha is treated as restricted, and cov[.,1] =
   cov[1,.] is returned as 0.  */
format /rdn 10,4;
local n, xx, k, beta, se, s2, dw, cln, thetaols, dth, view, theta;
local maxlik, its, alpha, c, glik, lrstat, tstat;
local h, dt, i, j, cov, hessian, var;
local e, n1, n2, ytrim, xtrim, thetatrm, trimlik, dt1, dt2;
local grad, g, hstar, d, newth, lambda, dlik, newlik;
local alfmin, alfmax, wald;
local delta, zeta, betas;
local f1, f2, aye, i1, lik;
local dmax, sqn;
declare xxglob_, yglob_;
declare lnden_;
alfmin = 0.84;
alfmax = 2.00;
n = rows(y);
xx = x;
if const == 1;
  xx = ones(n,1)~xx;
  endif;
k = cols (xx);
yglob_ = y;
xxglob_ = xx;
if printout;
  lprint "SYMMETRIC STABLE REGRESSION (SMSTRG)";
  lprint;
  lprint "OLS preliminary estimates";
  endif;
{beta, se, s2, dw} = myols(y,xx,0,printout);
var = s2*(n-k)/n;
cln = ln(sqrt(var/2));
@ theta is a column vector in smstrg@
thetaols = 2|cln|beta;
glik = -(n/2)*(ln(2*pi*var)+1);
theta = thetaols;
maxlik = glik;
if printout;
  lprint "  log L = " glik;
  endif;
sqn = sqrt(n);
dth = 1/sqn|1/sqn|se; @ (used for Hessian) @
e = y - xx * beta;
e = e~y~xx;
e = sortc(e,1);
n1 = round(n/4);
n2 = n + 1 - n1;
if n2-n1+1 < k+2; goto start; endif;
e = e[n1:n2,.];
ytrim = e[.,2];
xtrim = e[.,3:2+k];
if printout > 1;
  lprint "Trimmed OLS preliminary estimates";
  endif;
{beta, se, s2, dw} = myols(ytrim,xtrim,0,printout > 1);
e = y - xx * beta;
{alpha, betas, c, delta, zeta} = stab(e,0);
cln = ln(c);
thetatrm = alpha|cln|beta;
trimlik = srlikc(thetatrm);
if printout;
  lprint;
  lprint "Quantile preliminary estimates";
  lprint "  alpha = " alpha;
  lprint "  c =     " c;
  lprint "  log L = " trimlik;
  endif;
if trimlik > glik;
  theta = thetatrm;
  maxlik = trimlik;
  endif;

start:
its = 1;
if vu > 1;
  print its~theta';
  if vu > 2;
    print "waiting "; wait;
    endif;
  endif;

newit:
h = .001;
dt1 = h * dth;
dt2 = - dt1;
alpha = theta[1];
if alpha + dt1[1] > alfmax; dt1[1] = alfmax - alpha; endif;
if alpha + dt2[1] < alfmin; dt2[1] = alfmin - alpha; endif;
dt = dt1-dt2;
dt1 = diagrv(zeros(k+2,k+2),dt1);
dt2 = diagrv(zeros(k+2,k+2),dt2);
grad = zeros(k+2,1);
g = zeros(n,k+2);
j = 1;
do until j > k+2;
grad[j] = srlikc(theta + dt1[.,j]);
g[.,j] = lnden_;
grad[j] = (grad[j] - srlikc(theta + dt2[.,j]))/dt[j];
g[.,j] = (g[.,j] - lnden_) / dt[j];
j = j + 1;
endo;
hstar = g'*g;
if (theta[1] >= alfmax and grad[1] >= 0) or
   (theta[1] <= alfmin and grad[1] <= 0) ;
  d = 0|inv(hstar[2:k+2,2:k+2])*grad[2:k+2];
  if vu > 1; print "restricted d "; print d'; endif;
  else;
  d = inv(hstar)*grad;
  if vu > 1; print "regular d "; print d'; endif;
  endif;
dmax = maxc(abs(d)./dth);
newth = theta + d;
lambda = 1;
if newth[1] > alfmax;
  lambda = (alfmax-theta[1])/(newth[1]-theta[1]);
  if vu > 1; print "foreshortened d. lambda = "; print lambda; endif;
  newth = theta + lambda * d;
  elseif newth[1] < alfmin;
  lambda = (alfmin-theta[1])/(newth[1]-theta[1]);
  if vu > 1; print "foreshortened d (low alf) "; print lambda; endif;
  newth = theta + lambda * d;
  endif;

newlik:
its = its + 1;
newlik = srlikc(newth);
if vu > 1;
  print (its|theta|newlik)';
  if vu > 2; wait; endif;
  endif;
if newlik < maxlik and dmax > 1.e-6;
  if vu > 1; print "overshoot"; endif;
  lambda = lambda /2;
  newth = theta + lambda * d;
  goto newlik;
else;
  theta = newth;
  dlik = newlik - maxlik;
  maxlik = newlik;
  if dlik > tol * lambda;
     goto newit;
     endif;
  endif;

hessian:
alpha = theta[1];
cln = theta[2];
c = exp(cln);
beta = theta[3:k+2];
if vu; print "Computing Hessian"; endif;
h = .01;
dt = h * dth;
f1 = zeros(k+2,1);
f2 = f1;
hessian = zeros(k+2,k+2);
aye = eye(k+2);
i1 = 1;
if (alpha + abs(dt[1]) > alfmax or alpha - abs(dt[1]) < alfmin);
  i1 = 2;
  endif;
i = i1;
do until i > k+2;
f1[i] = srlikc(theta + dt[i] * aye[.,i]);
f2[i] = srlikc(theta - dt[i] * aye[.,i]);
i = i + 1;
endo;
i = i1;
do until i > k+2;
j = i1;
do until j > k+2;
@ if (i < 3 and j > 2) or (j < 3 and i > 2);
  /* for any symmetric distribution, information matrix is
  block diagonal */
  goto nextj;
else@
if i == j;
  hessian[i,j] = (2*maxlik - f1[i] - f2[i])/(-dt[i]^2);
else;
  lik = srlikc(theta + dt[i] * aye[.,i] - dt[j] * aye[.,j]);
  hessian[i,j] = (lik - f1[i] -f2[j] + maxlik) / (-dt[i]*dt[j]);
endif;
nextj:
j = j + 1;
endo;
i = i + 1;
endo;
@ This numerical Hessian may not be precisely symmetrical @
hessian = (hessian + hessian')/2;

if i1 == 1;
  cov = inv(-hessian);
else;
  cov = inv(-hessian[2:k+2,2:k+2]);
  cov = zeros(k+2,1)~(zeros(1,k+1)|cov);
endif;
se = sqrt(diag(cov));
wald = beta./se[3:k+2];

if printout;
  lprint;
  lprint "Symmetric Stable ML Estimates";
  lprint "  alpha = " alpha " se " se[1];
  lprint "  log c = " cln " se " se[2];
  lprint "  c =     " c;
  lprint "  log L = " maxlik;
  lprint "  LR stat " 2*(maxlik-glik);
  if const;
    lprint "  beta[1] is intercept";
    endif;
  lprint;
  lprint "      j      beta[j]      se     Wald ratio";
  lprint seqa(1,1,k)~beta~se[3:k+2]~wald;
  lprint;
  endif;
if vu;
  print;
  print "Symmetric Stable ML Estimates";
  print "  alpha = " alpha " se " se[1];
  print "  c =     " c;
  print "  log c = " cln " se " se[2];
  print "  log L = " maxlik;
  print "  LR stat " 2*(maxlik-glik);
  if const;
    print "  beta[1] is intercept";
    endif;
  print;
  print "      j      beta[j]      se     Wald ratio";
  print seqa(1,1,k)~beta~se[3:k+2]~wald;
  endif;
if vu == 3;
  wait;
  endif;
retp(alpha, c, beta, cov, its,maxlik,2*(maxlik-glik));
endp;

proc srlikc(theta);
/* theta = alpha|cln|beta
   receives y as yglob_, XX as XXglob_ (globals)
   passes likelihood contributions back as global lnden_
   theta is a column vector in this "c" version of srlik */
local alpha, c, beta, e;
alpha = theta[1];
c = exp(theta[2]);
beta = theta[3:rows(theta)];
e = yglob_ - xxglob_ * beta;
lnden_ = ln(symstb1(e/c, alpha)/c);
retp(sumc(lnden_));
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;

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 rndssta(alpha, r);
/* Returns rX1 vector of iid standard symmetric stable pseudo-random
   numbers with characteristic exponent alpha in [.1,2], and skewness
   parameter beta = 0, using method of Chambers, Mallows and Stuck
   (JASA 1976, 340-4).
   Encoded in GAUSS by J. Huston McCulloch, Ohio State University Econ Dept.
   (mcculloch.2@osu.edu), 12/95, directly from the CMS equation (4.1).
   Each r.v. has log characteristic function
        log E exp(ixt) = -abs(t)^alpha
   When alpha = 2, the distribution is Gaussian, with variance 2.
   When alpha = 1, the distribution is standard Cauchy.
   If alpha > 1, Ex = 0.  In all cases, the median is 0.
   This proc uses 2r uniform pseudo-random numbers from RNDU.  */

local phi, w;
if (alpha < .1) or (alpha > 2);
  format /rdn 10,4;
  print "Alpha must be in [.1,2] for proc RNDSSTA." ;
  print "Actual value is " alpha;
  if alpha < .1 and alpha > 0;
    print "If alpha < .1, probability of overflow becomes non-negligible.";
    endif;
  stop;
  endif;
w = - ln(rndu(r,1));
phi = (rndu(r,1) - .5) * pi;
if feq(alpha,2);
  retp(2 * sqrt(w) .* sin(phi));
elseif feq(alpha,1);
  retp(tan(phi));
endif;
retp((cos((1-alpha) * phi) ./ w) ^ (1/alpha - 1) .* sin(alpha * phi)
             ./ cos(phi) ^ (1/alpha));
endp;

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;

proc(4) = myols(y,x,const,out);
@ I had trouble getting proc OLS to work, and
  so wrote this. @
@ returns beta, se, s2, dw @
@ const = 1 for constant @;
@ out = 1 for printed output @;
local xx, n, k, xpxi, beta, yhat, e, ee, s2, se, dw, t, de;
n = rows(y);
if rows(x) /= n;
  print "rows(y) /= rows(x) " n rows(x);
  stop;
  endif;
xx = x;
if const;
  xx = ones(n,1)~x;
  endif;
k = cols(xx);
xpxi = inv(xx'*xx);
beta = xpxi * xx' * y;
yhat = xx * beta;
e = y - yhat;
ee = e'*e;
s2 = ee / (n-k);
se = sqrt(s2 * diag(xpxi));
de = e[1:n-1] - e[2:n];
dw = de'de / ee;
if out;
  format /rdn 5,0;
  print "  n " n;
  print "  k " k;
  format /rdn 10,4;
  t = beta./se;
  print "      j      beta[j]      se        t";
  print seqa(1,1,k)~beta~se~t;
  print;
  print "  dw " dw;
  format /rdn 5,0;
  lprint "  n " n;
  lprint "  k " k;
  format /rdn 10,4;
  lprint;
  lprint "      j      beta[j]      se        t";
  lprint seqa(1,1,k)~beta~se~t;
  lprint;
  lprint "  dw " dw;
  lprint;
endif;
retp(beta, se, s2, dw);
endp;

/* End of file smstrg */