Maximum Likelihood estimation with fminunc

Hello,

I would like to estimate two parameters using the maximum likelihood method. The log-likelihood function for a Poisson distribution is ∑D(t,x)ln[E(t,x)m(t,x,ϕ)]-E(t,x)m(t,x,ϕ)-ln[D(t,x)!] (the sum over x and t, at first the x from 1 to 106 for the fixed t, then the same for fixed t+1 on so on.), where ϕ represents the parameters k1 and k2 to be estimated. Given is E(t,x) and D(t,x) from the data. My idea:

EInitial_xt = xlsread('Exposures.xls'); 
E_xt = EInitial_xt(: ,5); %populatiion by age and year

DInitial_xt = xlsread('Deaths.xls');
D_xt = DInitial_xt(: ,5); %Death year, age



m = 106; % quantity of age

n = 181;

x = zeros(m, 1);

for i= 1:m

x(i, 1) = i;

end



organizedtableE = zeros(m, n);
for j = 1:n
    
    rowIndex = m * j;

    age = EInitial_xt(rowIndex, 1); 
    
    organizedtableE(1, j+1) = age;

end

for i = 1:m

    organizedtableE(i+1, 1) = i-1;
  
end



for k = 1:length(E_xt)
    age = mod(k-1, m) + 1; 
    yearIndex = ceil(k / m);  

    
    if yearIndex <= n
        organizedtableE(age+1, yearIndex+1) = E_xt(k);
    end
end

organizedtableD = zeros(m, n);
for j = 1:n
    
    rowIndex = m * j;

    year = DInitial_xt(rowIndex, 1); 
    
    organizedtableD(1, j+1) = year;

end

for i = 1:m

    organizedtableD(i+1, 1) = i-1;
  
end


for k = 1:length(E_xt)
    age = mod(k-1, m) + 1; 
    yearIndex = ceil(k / m);  

    
    if yearIndex <= n
        organizedtableD(age+1, yearIndex+1) = D_xt(k);
    end
end

cleanE = organizedtableE(2:end, 2:end);
cleanD = organizedtableD(2:end, 2:end);

[a, b] = size(cleanE);

for i = 1:181

options=optimoptions(@fminunc,'Algorithm','quasi-newton');

x0=[0, 0];

  fun = @(k) sum_l(k, (i-1)*a+1, (i)*a, EInitial_xt, D_xt, E_xt); % (i-1)*a+1 is the startIndex of the loop in the lower function and i*a is the endIndex, so both form the sum length.

    [X, fval] = fminunc(fun, x0, options);

    results(i, :) = [X, fval];

end




function ml = l(k, x, D, E)
k1 = k(1);
k2 = k(2);

ml = (-1)*(D*log(E*log(1+exp(k1+(x-52.5)*k2)))-E*log(1+exp(k1+(x-52.5)*k2))-log(sqrt(2*pi*D))-D*log(D/exp(1)));

end

function sum_ml = sum_l(k, startIndex, endIndex, EInitial_xt, D_xt, E_xt);

sum_ml = 0 ;

    for j = startIndex:endIndex

        sum_ml = sum_ml + l(k, EInitial_xt(j, 2), D_xt(j, 1), E_xt(j, 1)); % The idea here is that you have the sum by the recursive function, since E and D depend on the time t and the age x.

    end

end

The error that occurs:

Error using fminusub

Objective function is undefined at initial point. Fminunc cannot continue.

Error in fminunc (line 499)

[x,FVAL,GRAD,HESSIAN,EXITFLAG,OUTPUT] = fminusub(funfcn,x, ...

Error in B_Replikation (line 166)

[X, fval] = fminunc(fun, x0, options);

The problem is that 67 values come out, but not the full 106. The error message is that the function is not defined at this point, but I don't understand why. Also, the values for k2 are too low when compared to the values from the OLS method. What is wrong with the code?

Thanks in advance for your help.

Kind regards

[Dyuman Joshi]

Values/definitions of variables are missing.

Please provide the values so that we can run your code and recreate the error.

[Dave]

Do you need the excel sheet with the data as well?

[Dyuman Joshi]

Yes, that would be helpful as well.

[Torsten]

I don't have much experience with maximum likelihood problems, but why do you have to estimate the two parameters 181 times ?

[Dave]

Because I need to estimate the two parameters in every given t.

[Torsten]

Are you sure you need to estimate the values for each year separately ?

Answer 1

[Torsten]

D_xt(8162) = 0, and the log of 0 gives NaN. This causes problems for "fminunc".

So you should either delete rows where E_xt = 0 or D_xt = 0 or replace the value 0 by 1, e.g.

You can find the corresponding indices using

idx = find(D_xt==0)

[Dave]

It worked now, thank you very much!