Code covered by the BSD License

# STBL: Alpha stable distributions for MATLAB

Functions to handle one-dimensional alpha stable distributions

stblcdf(x,alpha,beta,gam,delta,varargin)
```function F = stblcdf(x,alpha,beta,gam,delta,varargin)
%F = STBLCDF(X,ALPHA,BETA,GAM,DELTA) returns the cdf of the stable
% distribtuion with characteristic exponent ALPHA, skewness BETA, scale
% parameter GAM, and location parameter DELTA, at the values in X.  We use
% the parameterization of stable distribtuions used in [2] - The
% characteristic function phi(t) of a S(ALPHA,BETA,GAM,DELTA)
% random variable has the form
%
% phi(t) = exp(-GAM^ALPHA |t|^ALPHA [1 - i BETA (tan(pi ALPHA/2) sign(t)]
%                  + i DELTA t )  if alpha ~= 1
%
% phi(t) = exp(-GAM |t| [ 1 + i BETA (2/pi) (sign(t)) log|t| ] + i DELTA t
%                                 if alpha = 1
%
% The size of P is the size of X.  ALPHA,BETA,GAM and DELTA must be scalars
%
%P = STBLCDF(X,ALPHA,BETA,GAM,DELTA,TOL) computes the cdf to within an
% absolute error of TOL.  Default for TOL is 1e-8.
%
% The algorithm works by computing the numerical integrals in Theorem
% 1 in [1] using MATLAB's QUADV function.
%
% If abs(ALPHA - 1) < 1e-5,  ALPHA is rounded to 1.
%
%
% References:
%
% [1] J. P. Nolan (1997)
%     "Numerical Calculation of Stable Densities and Distribution
%     Functions"  Commun. Statist. - Stochastic Modles, 13(4), 759-774
%
% [2] G Samorodnitsky, MS Taqqu (1994)
%     "Stable non-Gaussian random processes: stochastic models with
%      infinite variance"  CRC Press
%
%
%

if nargin < 5
error('stblcdf:TooFewInputs','Requires at least five input arguments.');
end

% Check parameters
if alpha <= 0 || alpha > 2 || ~isscalar(alpha)
error('stblcdf:BadInputs',' "alpha" must be a scalar which lies in the interval (0,2]');
end
if abs(beta) > 1 || ~isscalar(beta)
error('stblcdf:BadInputs',' "beta" must be a scalar which lies in the interval [-1,1]');
end
if gam < 0 || ~isscalar(gam)
error('stblcdf:BadInputs',' "gam" must be a non-negative scalar');
end
if ~isscalar(delta)
error('stblcdf:BadInputs',' "delta" must be a scalar');
end

if nargin > 6
error('stblcdf:TooManyInputs','Accepts at most six input arguments.');
elseif isempty(varargin)
tol = 1e-8;
elseif isscalar(varargin{1})
tol = varargin{1};
else
end

% Warn if alpha is very close to 1 or 0
if (1e-5 < abs(1 - alpha) && abs(1 - alpha) < .02) || alpha < .02
warning('stblcdf:ScaryAlpha',...
'Difficult to approximate cdf for alpha close to 0 or 1')
end

%========= Compute CDF =============%

% Check to see if you are in a simple case, if so be quick, if not do
% general algorithm
if alpha == 2                  % Gaussian distribution
x = (x - delta)/gam;       % Standardize
F = .5*(1 + erf(x/2));     % ~ N(0,2)

elseif alpha==1 && beta == 0   % Cauchy distribution
x = (x - delta)/gam;       % Standardize
F = 1/pi * atan(x) + .5;

elseif alpha == .5 && abs(beta) == 1 % Levy distribution
x = (x - delta)/gam;  % Standardize
F = zeros(size(x));
if beta > 0
F(x > 0) = erfc(sqrt(1./(2*x(x>0))));
F(x <= 0) = 0;
else % beta < 0
F(x < 0) = 1 - erfc(sqrt(-1./(2*x(x<0))));
F(x >= 0) = 1;
end

elseif abs(alpha - 1) > 1e-5         % Gen. Case, alpha ~= 1

xold = x;  % Save for possible later use
% Standardize in (M) parameterization ( See equation (2) in [1] )
x = (x - delta)/gam - beta * tan(alpha*pi/2);
F = zeros(size(x));
% Compute CDF
zeta = -beta * tan(pi*alpha/2);
theta0 = (1/alpha) * atan( beta * tan(pi*alpha/2) );
A1 = alpha*theta0;
c1 = (alpha > 1) + (alpha < 1)*(1/pi)*(pi/2 - theta0);
A2 = cos(A1)^(1/(alpha-1));
exp1 = alpha/(alpha-1);
alpham1 = alpha - 1;
V = @(theta) A2 * ( cos(theta) ./ sin( alpha*(theta + theta0) ) ).^exp1.*...
cos( A1 + alpham1*theta ) ./ cos(theta);

if any(x > zeta)
xshift = (x(x>zeta) - zeta).^(alpha/(alpha - 1));
% shave off end points of integral to avoid numerical instability
% when calculating V
F( x > zeta ) = c1 + sign(1-alpha)/pi * ...
end

if any(abs(x - zeta) < 1e-8)
F(abs(x - zeta) < 1e-8) = (1/pi) * (pi/2 - theta0);
end

if any( x < zeta)
% Recall with -xold, -beta, -delta
F(x < zeta) =...
1 - stblcdf(-xold(x < zeta),alpha,-beta,gam,-delta);
end
elseif beta > 0                     % Gen. Case, alpha = 1, beta >0
x = (x - (2/pi) * beta * gam * log(gam) - delta)/gam; % Standardize
piover2 = pi/2;
twooverpi = 2/pi;
oneoverb = 1/beta;
% Use logs to avoid overflow/underflow
logV = @(theta) log(twooverpi * ((piover2 + beta *theta)./cos(theta))) + ...
( oneoverb * (piover2 + beta *theta) .* tan(theta) );

xterm = (-pi*x/(2*beta));
F = (1/pi)*quadv(@(theta) exp(-exp(xterm + logV(theta))),...
-pi/2+1e-12,pi/2-1e-12,tol);

else                           % alpha = 1, beta < 0
F = 1 - stblcdf(-x,1,-beta,gam,-delta,tol);
end

F = max(real(F),0); % in case of small imaginary or negative resutls from QUADV

end

```