Using fminsearch to estimate utility function with MLE

Hello.

I have a slight issue which I hope you guys can help me with. I want to estimate the alpha (a) of a simple utility function of the form

EU(P,V) = P*V^a

where P is probability of outcome, and V is magnitude of outcome. The response variable (riskChoice) is 1 for risky choice, and 0 for non-risky choice. I want to use maximum likelihood with a logit link-function to estimate the binary choice model

1/1(e^(w*(EUrisk^a-EUcertain^a))

where EUrisk(cert) is the expected value of the risky(certain) option at each timepoint t, and w is some weight. From the above, I am interested in estimating w and a, using MLE.

I have the following in MATLAB

    F = @(a,w) (1/(1 + exp(w*(euCert.^a - euRisk.^a))))';

    logF1 = @(a,w) log(1/(1 + (exp(w*(euCert.^a - euRisk.^a)))))';

    logF2 = @(a,w) log(1-F(a,w));

And the (negative, since I want a max.) log-likelihood looks like

    NLL = @(a,w) -(sum(riskChoice.*logF1(a,w)+(1-riskChoice).*logF2(a,w)));

Problem No. 1:

When evaluating NLL in a certain (theoretically realistic) point, such as NLL(0.2,2) i get a NaN response. Actually, I get NaN responses for almost all values I put in there. I'm thinking that I have made some mistake in writing up the model? I can't seem to figure out the mistake!

Problem No. 2:

When I use the code below, I get the following error:

    initParam = zeros(100,numel(euRisk));

    paramFrom = -10:0.1:10;

    allOptWs=zeros(100,2);
    modelEvi=zeros(100,1);

    for r = 1:100

        a=rand*10;
        w=rand*20;
        initWs = [a w];

        optWs = fminsearch(NLL,initWs);
        allOptWs(r,:) = optWs;
        modelEvi(r) = NLL(optWs);

    end

    Error using @(a,w)-(sum(riskChoice.*logF1(a,w)+(1-riskChoice).*logF3(a,w)))
    Not enough input arguments.

    Error in fminsearch (line 191)
    fv(:,1) = funfcn(x,varargin{:});

I think this might have something to do with fminsearch only accepting one array? If this suspicion is correct, how exactly do I go about making sure that both a and w are in one array?

Thanks in advance!

Answer 1

[Matt J]

Problem No. 1:

Maybe because NLL references logF3 instead of logF2? You haven't defined logF3 anywhere that I can see. My guess is that either the expression riskChoice.*logF1(a,w) or the expression (1-riskChoice).*logF3(a,w) is creating a 0*log(0) situation. Bear in mind that the log likelihood needs special handling whenever the likelihood is zero. Else you get log(0) situations.

In any case, it should be an easy matter to step through the computation of NLL, since you already know a particular (a,w) pair where it blows up, and find where the NaN is produced.

Problem No. 2:

You have written NLL as a function of two input arguments (a,w). The unknowns must be fed in as a vector,

    F = @(p) (1/(1 + exp(p(2)*(euCert.^p(1) - euRisk.^p(1)))))';
    logF1 = @(p) log(F(p));
    logF2 = @(p) log(1-F(p));

    NLL = @(p) -sum( riskChoice.*logF1(p) + (1-riskChoice).*logF2(p) ) ;

[Matt J]

Comment by Tobias:

Hey Matt, and thanks for the swift reply and clearing up my second issue.

However, F2 = F3 in my code. It was just a typo in the text abov.

I will try to take your advice and step through NLL to see if I get any 0*log0 situations.

best, Tobias

Answer 2

[Tobias]

So, I've implemented Matts suggestion, but without much of a result. When I search NLL with a loop going from 1 to 100 by integers, or 1 to 10 by 0.1 decimals ect. I get the following error message:

Warning: Rank deficient, rank = 0, tol =           inf. 
> In @(p)(1/(1+exp(p(2)*(euCert.^p(1)-euRisk.^p(1)))))'
  In @(p)log(1-F(p))
  In @(p)-sum(riskChoice.*logF1(p)+(1-riskChoice).*logF2(p))

Any ideas?

[Matt J]

Use element-wise division. In fact, use element-wise everything,

    >> F=vectorize(F)

    F =

    @(p)(1./(1+exp(p(2).*(euCert.^p(1)-euRisk.^p(1)))))

Answer 3

[Tobias]

Just wanted to say thanks for the help Matt. And add that I got some "close to zero" problems, solved using eps. Code is below.

    F = @(p) (1./(1 + exp(p(2)*(euCert.^p(1) - euRisk.^p(1)))));
    logF1 = @(p) log(F(p)+eps);
    logF2 = @(p) log(1-F(p)+eps);
    NLL = @(p) -sum( riskChoice.*logF1(p) + (1-riskChoice).*logF2(p) ) ;

Code for fminsearch:

    initParam = zeros(100,numel(euRisk));

    paramFrom = -10:0.1:10;

    allOptWs=zeros(100,2);
    modelEvi=zeros(100,1);

    for r = 1:100

        a=rand*3;
        w=(rand-0.5)*20;
        initWs = [a w];

        optWs = fminsearch(NLL,initWs);
        allOptWs(r,:) = optWs;
        modelEvi(r) = NLL(optWs);

    end

    [minModelEvi,minidx]=min(modelEvi);
    optWs=allOptWs(minidx,:)

best, tobias