//////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// '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_)

pi=3.1415926535897931;

if ~(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_;
endfunction



function [ret]=stablecalc(xarg, alpha_, beta_)

    function ret = fun(x)
    
        pi=3.1415926535897931;

        zero = abs(alpha_*(theta0+x))<10e-20 | abs(x-pi/2)<10e-20 | sin(alpha_*(theta0+x))<0;
        ret(zero)=0;
        x=x(~(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(~(zero)) = -gval';
        end
    endfunction


  function [ret,g,ind] = fun2(x,ind)
    ret=1000+fun(x);
    ind=4;
    g=numdiff(fun,x);
  endfunction
  
  
pi=3.1415926535897931;

zeta_ = -beta_ * tan(pi*alpha_/2);
theta0 = (1/alpha_) * atan( beta_ * tan(pi*alpha_/2))

disp(-theta0)
disp(pi/2)
disp((-theta0+pi/2)/2)
disp(alpha_)
disp(beta_)
disp(delta_)
disp(gamma_)
disp(fun(-theta0:0.01:pi/2))
[a1,b1,c1]=fun2(0.5,0)
disp([a1,b1,c1])
//optim(fun2,'b',-theta0,pi/2,(-theta0+pi/2)/2)
// two integrals rather than one, in order not to miss the mass
// (xmin it the point where @fun is minimum)

[b,xmin]=optim(fun2,'b',-theta0,pi/2,(-theta0+pi/2)/2);


//disp('xmin'); disp(xmin)

integral1 = -intg(-theta0,xmin,fun);
integral2 = -intg(xmin,pi/2,fun);

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

endfunction

    

function ret=verify_limits(alpha_, beta_, gamma_)
ret=1;
if alpha_<1.1 | alpha_>2
    warning('alpha_ out of bounds (1.1<alpha_<2)'); ret=0;
end
if beta_<-1 | beta_>1
    warning('beta_ out of bounds (-1<beta_<1)'); ret=0;
end
if gamma_<0
    warning('gamma out of bounds (gamma>0)'); ret=0;
end
endfunction