I always had problems matching llf of the linear models.
I ran into this again, and finally figured out that the scale estimate is MLE based with denominator nobs, instead of the denominator in results.scale which is (nobs - k_vars)
>>> stats.norm.logpdf(res.resid, scale=np.sqrt((res.resid**2).mean())).sum()
Even though the scale estimate is different, I think llf is exactly the same as Stata reports (IIRC)
At several places I needed a normal linear model that does not concentrate out the scale/sigma_squared estimate. We should get a basic generic Normal Model class for this.
just to clarify what I think I was thinking when I opened this:
The scale estimate in the linear models uses the unbiased least squares, while the concentrated loglike uses the definition of the MLE. So llf is not the loglike at [params, scale] but at [params, scale_MLE]
related discussion in #1170