%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 'dstable' : computes a stable distribution density
%
% see limitations in 'verify_limits' function below
%
%
% Example:
% density=dstable(-3:0.1:3,1.7,0.5,1,0)
% plot(-3:0.1:3,density)
%
%
% Written by Jeremie Cosmao -- jcosmao@monaco.edu -- June 2007
% Hedge Fund Research Institute (HFRI) -- http://www.hfri.org
% International University of Monaco -- http://www.monaco.edu
% 
% inspired by the fBasics package of R CRAN, written by Diethelm Wurtz - Institut fur Theoretische Physik
%
% This program is free software; you can redistribute it and/or modify it under the terms of the
% GNU Library General Public License as published by the Free Software Foundation; either
% version 2 of the License, or (at your option) any later version. This library is distributed in
% the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied
% warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU Library General Public License for more details.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function ret = dstable(x, alpha, beta, gamma_, delta)

if not(verify_limits(alpha, beta, gamma_))
    ret=0;
    return
end

x = (x - delta) / gamma_;
zeta_ = -beta * tan(pi*alpha/2);

for z=1:length(x)
    if (abs(x(z) - zeta_) < 0.001)
        theta0 = (1/alpha) * atan( beta * tan(pi*alpha/2));
        ret(z) = gamma(1+1/alpha)*cos(theta0) / (pi*(1+zeta_^2)^(1/(2*alpha)));
    else
        if (x(z) > zeta_)
            ret(z) = stablecalc(x(z), alpha, beta);
        end
        if (x(z) < zeta_)
            ret(z) = stablecalc(-x(z), alpha, -beta);
        end
    end
end

ret=ret/gamma_;
end


function [ret]=stablecalc(xarg, alpha, beta)

zeta_ = -beta * tan(pi*alpha/2);
theta0 = (1/alpha) * atan( beta * tan(pi*alpha/2));

% two integrals rather than one, in order not to miss the mass
% (xmin it the point where @fun is minimum)
xmin = fminbnd(@fun,-theta0, pi/2);
integral1 = -quadl(@fun,-theta0,xmin);
integral2 = -quadl(@fun,xmin,pi/2);

% interv=-theta0:0.001:pi/2;
% figure(26); plot(interv,fun(interv)); hold all; drawnow;

ret=alpha / (pi*abs(alpha-1)*(xarg-zeta_)) * (integral1 + integral2);

    function ret = fun(x)
        zero = abs(alpha*(theta0+x))<10^-20 | abs(x-pi/2)<10^-20 | sin(alpha*(theta0+x))<0;
        ret(zero)=0;
        x=x(not(zero));
        if length(x)>0
            v = (cos(alpha*theta0))^(1/(alpha-1)) *  (cos(x)./sin(alpha*(theta0+x))).^(alpha/(alpha-1)) .*  (cos(alpha*theta0+(alpha-1)*x)./cos(x));
            
            g = (xarg-zeta_)^(alpha/(alpha-1)) * v;
            gval = g .* exp(-g);

            ret(not(zero)) = -gval;
        end
    end
end


function ret=verify_limits(alpha, beta, gamma_)
ret=true;
if alpha<1.1 | alpha>2
    warning('alpha out of bounds (1.1<alpha<2)'); ret=false;
end
if beta<-1 | beta>1
    warning('beta out of bounds (-1<beta<1)'); ret=false;
end
if gamma_<0
    warning('gamma out of bounds (gamma>0)'); ret=false;
end
end