Logit panel with Simulated Maximum Likelihood. Convergence question with fminunc

I'm testing a very basic logit panel model. The setup is as follows: We observe a binary variable where i is individual and t is time. For simplicity two periods. It is assumed that the true parameters are and and ε is an independent logit error. It's assumed that is individual specific but the same accross time.

I coded this, the problem is that when passing the function thorugh fminunc, the algorithm runs only one or two iterations and returns exactly the same values I provide as initial parameters. I coded this in other languages, but in this case I cannot see the problem.

I've tried different algorithms, other initial values, etc. The info after fminunc is pretty much the same:

Iteration Func-count f(x) Step-size optimality

0 4 1386.49 2.42e+07

1 64 1386.25 6.57442e-15 2.62e+07

Local minimum possible.

fminunc stopped because the size of the current step is less than the value of the step size tolerance.

My exact code is:

clear;

%%to 2 and the algorithm to Mersenne Twister

rng(66,'twister');

% goal: estimate the panel logit

%% DGP

T = 2; % maximum number of periods per agent

N = 1000; % number of agents

x = randn(N, T); % covariates: observables

beta_aux = 2+randn(N, 1); % random slope; agent specific but same in time

beta_i = [beta_aux, beta_aux];

beta_0 = 1;

z = beta_0 + beta_i .* x; % Calculate linear predictors

p = 1 ./ (1 + exp(-z)); %Calculate choice probs for each individual at every time (because independence)

% Generate binary outcomes based on probabilities

y = binornd(1, p); % N x T matrix of binary outcomes

%%

options = optimset('Display', 'iter', 'MaxIter', 1000);

negLogLik=@(params) logit_loglik(params,x,y,N,T);

params_init=[0, 0, 0];

[params_hat, fval] = fminunc(negLogLik, params_init, options);

%% Functions

function loglik = logit_loglik(params, x, y,N,T)

%unpack params

beta_0=params(1);

mu_beta=params(2);

aux_sigma=params(3);

sum_py=zeros([N,T]);

sigma_beta=exp(aux_sigma); %To ensure sigma is greater than 0

%Montecarlo draws

M=10000;

%We first need to get the unconditional choice probabilities by integrating

%over values of beta_i

for m = 1:M

beta_aux=mu_beta+sigma_beta*randn(N, 1);

beta_i = [beta_aux,beta_aux];

z = beta_0 + beta_i .* x;

% Calculate probabilities for all individuals, time periods, and MC simulations

py = 1 ./ (1 + exp(-z));

sum_py=sum_py+py;

end

%approximated integral of the choice probability Pr(y=1) for each

%person/time/choice

integrated_py=sum_py/M;

integrated_PY=log(y.*(integrated_py) + (1 - y).* (1 - integrated_py));

%We sum first over time

Lt=sum(integrated_PY,2);

%Then over individuals

loglik=-sum(Lt);

end

Answer 1

[Torsten]

Your objective function depends on random values. This is not allowed for deterministic optimization as done by "fminunc". Convergence of these solvers depends on the differentiability of the objective function that gets lost in this case.

To estimate distribution parameters by maximum likelihood estimate, use "mle":

https://uk.mathworks.com/help/stats/mle.html

[Marcos]

Appreciate your answer. That makes totally sense. I will proceed as suggested

[Marcos]

I will make a follow up question as documentation for mle() is not quite intuitive. Everything above is the same. I'm trying to use the mle() function but it is not intuitive at all. mle() requires a single data vector. In my case, my negative log likelihood function needs x and y. I vectorized x, y to pass them to in the function but its not woking at all.

I don't see the proble. Is there any other maximization function that allow me to do Maximum Simulated Likelihood as I am proposing?

clear;

%%to 2 and the algorithm to Mersenne Twister

rng(66,'twister');

% goal: estimate the panel logit

%% DGP

T = 2; % maximum number of periods per agent

N = 1000; % number of agents

x = randn(N, T); % covariates: observables

beta_aux = 2+randn(N, 1); % random slope; agent specific but same in time

beta_i = [beta_aux, beta_aux];

beta_0 = 1;

z = beta_0 + beta_i .* x; % Calculate linear predictors

p = 1 ./ (1 + exp(-z));

% Generate binary outcomes based on probabilities

y = binornd(1, p); % N x T matrix of binary outcomes

% Combine x and y into a single matrix for mle

data = [x(:); y(:)];

%%

negLogLik=@(p1,p2,p3,data) logit_loglik(p1,p2,p3,data);

params_init=[0, 0, 0];

% Set options for the optimization

options = statset('MaxIter', 1000, 'Display', 'iter');

dummy_data = zeros(N*T, 1);

% Estimate parameters using mle with custom negative log-likelihood function

[params_hat, ~] = mle(data, 'nloglf', negLogLik, 'Start', params_init, 'options', options);

% Display the estimated parameters

disp(params_hat);

%% Functions

function loglik = logit_loglik(p1,p2,p3, data)

x=[data(1:1000), data(1001:2000)];

y=[data(2001:3000), data(3001:4000)];

N = size(x, 1);

T = size(x, 2);

%unpack params

beta_0=p1;

mu_beta=p2;

aux_sigma=p3;

sum_py=zeros([N,T]);

sigma_beta=exp(aux_sigma); %To ensure sigma is greater than 0

%Montecarlo draws

M=10000;

%We first need to get the unconditional choice probabilities by integrating

%over values of beta_i

for m = 1:M

beta_aux=mu_beta+sigma_beta*randn(N, 1);

beta_i = [beta_aux,beta_aux];

z = beta_0 + beta_i .* x;

% Calculate probabilities for all individuals, time periods, and MC simulations

py = 1 ./ (1 + exp(-z));

sum_py=sum_py+py;

end

%approximated integral of the choice probability Pr(y=1) for each

%person/time/choice

integrated_py=sum_py/M;

integrated_PY=log(y.*(integrated_py) + (1 - y).* (1 - integrated_py));

%We sum first over time

Lt=sum(integrated_PY,2);

%Then over individuals

loglik=-sum(Lt);

end

[Torsten]

I don't have the background knowledge to follow what you are trying to do.

For "mle", usually realizations of a random variable are known (thus one vector of values). You guess that these realizations follow a certain distribution (from which usually the pdf is known). "mle" then tries to fit the distribution parameters such that the maximum likelihood function attains its maximum.