%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 'stablemlefit' computes the Maximum Likelihood estimates of a stable distribution
%
%
%
% Example:
% data=stabrnd(1.7,0,1,0,500,1);
% stabmlefit(data)
%
% 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] = stabmlefit(data)

tic
params=stabfit(data,0);
strcat('initial params:',num2str(params,5))

hist(data,min(data):(max(data)-min(data))/30:max(data))
hold all
a.fval=0;
outfun(params, a, 0);

%warning off
[x,fval,exitflag,output] = fminsearch(@sumlog,params, optimset('OutputFcn',@outfun))
%warning on

ret=x;
toc

    function ret=sumlog(par)
        l=dstable(data, par(1),par(2),par(3),par(4));
        if sum(l==0)>0
            ret=Inf;
        else
            ret=-sum(log(l));
        end
    end

    function stop = outfun(x, optimValues, state)
        %x
        stop=false;
        xx=dstable(data,x(1),x(2),x(3),x(4));
        plot(data,xx*length(data)*((max(data)-min(data))/30),'.'); title(num2str([x optimValues.fval],3));
        drawnow
    end
end