Hmm, I need to noodle on this. I really appreciate the thought you’ve put in!
I think what I possibly fail to understand is how a distributional assumption on u_i is any different than placing a distributional assumption on a “standard” intercept (see above where there is \beta_{0,i} and u_i in the model). Is it because we can just define u_i = \beta_{0,i} + v_i, where v_i is the unobserved heterogeneity and \beta_{0,i} is the constant? Maybe I’m overthinking things.
Update: I think what you’re doing is using \alpha in your code as a “global” constant (kind of capturing the global mean across groups of y), but augmenting the intercept by u_i, which captures the unobserved heterogeneity. A distributional assumption is made on \alpha, too, but again we treat this as common to all groups. This is starting to make more sense, but I need to think about it some more. It also occurs to me that if the u_i are fixed rather than random effects, then the u_i will be correlated with a subset, if not all, of the predictors. This requires a different solution, likely in the form of a transformation or IV or both.
Were you able to recover the true parameters in your simulation based on the DGP?