Thread Subject:

Help on fitting a distribution to data with fminunc

Hi All,

I want to fit a skewed-t distribution to some return data (a single series of 2000 observations) in order to obtain the maximum likelihood estimates of the degrees of freedom (nu) and skewness (lambda) parameters.

When I run the code below I get the following:

"fminunc stopped because it cannot decrease the objective function
along the current search direction.

<stopping criteria details>


nu =

  1.0e-007 *

    0.1512
         0


lam =

    1.5154"


The nu parameter is too low and I don't know if I am on track,
Any suggestions on improving the code in order to get some reliable parameter estimates would be greatly appreciated.

Thank you for your time!


%%%%%%%%%%%%%%%%%code%%%%%%%%%%%%%%%%%%%%%%%

function [nu,lambda]=EstSktParmsc(data)
%% Optimisation Parameters%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Tol = 1e-7;
T = size(data,1);
thita=[6;0];
% Initial estimates
x0 = [thita];
%% Optimisation Options Set - fmincon
options = optimset('fminunc');
options = optimset(options , 'Algorithm' , 'trust-region-reflective');
options = optimset(options , 'Display' , 'iter');
options = optimset(options , 'Diagnostics' , 'on');
options = optimset(options , 'LargeScale' , 'off');
options = optimset(options , 'MaxFunEvals' , 5000);
options = optimset(options , 'MaxIter' , 5000);
options = optimset(options , 'TolFun' , 1e-12);
options = optimset(options , 'TolCon' , 1e-12);
options = optimset(options , 'TolX' , 1e-12);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Optimisation Subroutine%%%%%%%%%%%
[nu,lambda]=fminunc(@skewtdis_llf,x0,options);
%% Objective function(Log was obtained from Andrew Patton's site)%%%%%%%%
    function [LL, likelihoods] = skewtdis_llf(data, theta)
theta=[6;0];
        nu = theta(1);
lambda = theta(2);
[T,k] = size(data);
logc = gammaln((nu+1)/2) - gammaln(nu/2) - 0.5*log(pi*(nu-2));
c = exp(logc);
a = 4*lambda.*c.*((nu-2)./(nu-1));
logb = 0.5*log(1 + 3*lambda.^2 - a.^2);
b = exp(logb);
find1 = (data<(-a./b));
find2 = (data>=(-a./b));
LL1 = logb + logc - (nu+1)/2.*log(1+1./(nu-2).*((b.*data(find1)+a)./(1-lambda)).^2);
LL2 = logb + logc - (nu+1)/2.*log(1+1./(nu-2).*((b.*data(find2)+a)./(1+lambda)).^2);
LL = sum(LL1) + sum(LL2);
LL = -LL;

if nargout>1
     likelihoods=zeros(size(data));
     likelihoods(find1)=LL1;
     likelihoods(find2)=LL2;
end
    end
end

"Chris" wrote in message <ju974t$fpr$1@newscl01ah.mathworks.com>...
> Hi All,
>
> I want to fit a skewed-t distribution to some return data (a single series of 2000 observations) in order to obtain the maximum likelihood estimates of the degrees of freedom (nu) and skewness (lambda) parameters.
===============

Hard to say, but here are a few issues,

(1) You are calling fminunc seemingly in disregard of the input/output syntax in the documentation. In particular, lambda is one of your unknown variables, yet you are returning it as the fval output argument. All unknowns should be bundled together in a single vector.

(2) The unknown loglikelihood parameters theta are passed as the 2nd argument to skewtdis_llf. FMINUNC expects it to be the first argument. A solution would be to wrap it in an anomyous function
 fminunc(@(theta) skewtdis_llf(data,theta),x0,options).

(3) It is not clear how the "data" variable gets passed to skewtdis_llf in the code you've shown. Similarly, your output is a variable "lam", which appears nowhere in the code. This leads me to think we're not seeing the code you're actually running, making it hard for us to help you troubleshoot...


In any case, when you have problems with FMINUNC or the other solvers, you should call them with all output arguments, including exitflag, so that the output conditions can be studied.

"Matt J" wrote in message <ju996k$oc1$1@newscl01ah.mathworks.com>...
>
> Hard to say, but here are a few issues,
>
> (1) You are calling fminunc seemingly in disregard of the input/output syntax in the documentation. In particular, lambda is one of your unknown variables, yet you are returning it as the fval output argument. All unknowns should be bundled together in a single vector.
>
> (2) The unknown loglikelihood parameters theta are passed as the 2nd argument to skewtdis_llf. FMINUNC expects it to be the first argument. A solution would be to wrap it in an anomyous function
> fminunc(@(theta) skewtdis_llf(data,theta),x0,options).
>
> (3) It is not clear how the "data" variable gets passed to skewtdis_llf in the code you've shown. Similarly, your output is a variable "lam", which appears nowhere in the code. This leads me to think we're not seeing the code you're actually running, making it hard for us to help you troubleshoot...
>
>
> In any case, when you have problems with FMINUNC or the other solvers, you should call them with all output arguments, including exitflag, so that the output conditions can be studied.

Dear Matt,

thank you for your help,

I think I have addressed the above mentioned inconsistencies and now I have to do the fine tuning.

The code now has as follows:

function [x,fval,exitflag]=estskewt(data)

%% Objective function
    function [LL, likelihoods] = skewt_llf(pars, data)
nu = pars(1); %Degrees of Freedom parameter
lambda = pars(2); %Skewness parameter
[T,k] = size(data);

logc = gammaln((nu+1)/2) - gammaln(nu/2) - 0.5*log(pi*(nu-2));
c = exp(logc);
a = 4*lambda.*c.*((nu-2)./(nu-1));
logb = 0.5*log(1 + 3*lambda.^2 - a.^2);
b = exp(logb);

find1 = (data<(-a./b));
find2 = (data>=(-a./b));
LL1 = logb + logc - (nu+1)/2.*log(1+1./(nu-2).*((b.*data(find1)+a)./(1-lambda)).^2);
LL2 = logb + logc - (nu+1)/2.*log(1+1./(nu-2).*((b.*data(find2)+a)./(1+lambda)).^2);
LL = sum(LL1) + sum(LL2);
LL = -LL;

if nargout>1
     likelihoods=zeros(size(data));
     likelihoods(find1)=LL1;
     likelihoods(find2)=LL2;
end
    end
%% Optimisation Parameters & Starting Values
nu=5;
lambda=0;
pars=[nu;lambda];
parsL=[0;-15];
parsU=[Inf;15];
% Initial estimates
x0 = [5;0];
%% Optimisation Options Set - fmincon
optopt=optimset('Algorithm', 'interior-point', 'MaxIter',100,'Display','Iter','LargeScale','on')

%% Optimisation Subroutine
[x,fval,exitflag]=fmincon(@skewt_llf,pars, [],[],[],[],parsL,parsU,[],optopt,data);
pars=x;
[LL,likelihoods]=skewt_llf(pars,data);

end

After running it with my data (15 iterations) what I get is:
"Local minimum possible. Constraints satisfied.

fmincon stopped because the size of the current step is less than
the default value of the step size tolerance and constraints were
satisfied to within the default value of the constraint tolerance.

<stopping criteria details>


parameters =

    2.0001
   -0.0314


fval =

 -5.8077e+003


exitflag =

     2"

"Chris" wrote in message <ju9ma4$j7j$1@newscl01ah.mathworks.com>...
>
> After running it with my data (15 iterations) what I get is:
> "Local minimum possible. Constraints satisfied.
>
> fmincon stopped because the size of the current step is less than
> the default value of the step size tolerance and constraints were
> satisfied to within the default value of the constraint tolerance.
>
> <stopping criteria details>
>
>
> parameters =
>
> 2.0001
> -0.0314
>
>
> fval =
>
> -5.8077e+003
>
>
> exitflag =
>
> 2"
=================

Ok, but what exactly is the issue now? FMINCON is telling you that it thinks it succeeded.

Incidentally, by returning "likelihoods" as the 2nd output argument from your objective function, you are preventing yourself from turning on the 'GradObj' option and supplying your own calculated gradient. Indeed, it looks like computing the analytical gradient wouldn't be a hard calculation and would reduce the time FMINCON spends on finite differences.

"Matt J" wrote in message <juekmp$86n$1@newscl01ah.mathworks.com>...
> "Chris" wrote in message <ju9ma4$j7j$1@newscl01ah.mathworks.com>...

>
> Incidentally, by returning "likelihoods" as the 2nd output argument from your objective function, you are preventing yourself from turning on the 'GradObj' option and supplying your own calculated gradient. Indeed, it looks like computing the analytical gradient wouldn't be a hard calculation and would reduce the time FMINCON spends on finite differences.

Thanks a million, I think it is ok, parameter estimates make sense.