Consider the following static panel model with two-way fixed effects:

$y_{it}=β z_{it}+ δ_i+δ_t+ϵ_{it}$

I am interested in the estimate of $\beta$. In my sample, $N=600,000$, $T=100$ and $I=6,000$. The model thus contains 6,100 dummy variables to be estimated. $y_{it}$ is a rare outcome variable with around 4,000 non-zero values. I ran models with different transformations: a 0-1 dummy, a simple count (from 0 to 20), a count of log(1+values), or the values scaled over population. (I am aware none of these are perfect.)

Given the large number of dummy variables, I have so far estimated models via OLS. $\hat{\beta}$ is always negative and highly statistically significant. I was asked to estimate this model using estimators for count data (poisson, negative binomial) as well as logit. For logit, Fernandez-Val and Weidner (2016) propose a logit estimator that allows for two-way fixed effects, but this does not converge.

It is often said (and written on blogs) that nonlinear models are not well-suited for panel data. I have not seen a discussion for settings with many multi-dimensional fixed effects. In this case, where the number of cross-sections is much larger than $T$, how should I think about incidental parameter bias? Which other reference is there to tell me what makes sense and what doesn't?

Thanks a lot

Reference:

Fernandez-Val, Ivan and Weidner, Martin, (2016), Individual and time effects in nonlinear panel models with large N, T, Journal of Econometrics, 192, issue 1, p. 291-312.