I'm testing a very basic logit panel model in matlab. The setup is as follows: We observe a binary variable $y_{it} = 1(\beta_0 + \beta_{it}x_{it} + \varepsilon_{it} > 0)$ where i is individual and t is time. For simplicity two periods. It is assumed that the true parameters are $\beta_0$ and $ \beta_{it} \sim \mathcal{N}(2, 1)$ and $\varepsilon$ is an independent logit error. It's assumed that $\beta_i$ is individual specific but the same across time.
I coded this in matlab. 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
```