Custom maximum likelihood models

When no built-in model fits your data-generating process, writing the likelihood gives you estimates and standard errors anyway.

Maximum likelihood estimation picks parameters that make the observed data most probable. statsmodels' GenericLikelihoodModel (in statsmodels.base.model) turns a log-likelihood into a full model with the usual results object.

Subclass it and implement loglikeobs(params), returning the log-likelihood of each observation as an array (or nloglikeobs, its negative). Then:

fit(start_params=..., method="bfgs", maxiter=...) optimises the total likelihood,
standard errors come from the inverse Hessian, computed numerically unless you supply score and hessian,
summary(), conf_int(), t_test(), wald_test() and cov_type="HC0" (sandwich) standard errors all work.

Implementing per-observation log-likelihoods enables the robust (sandwich) covariance, which uses each observation's score contribution.

Practical points:

Reparameterise constrained parameters: estimate log(sigma) instead of sigma so the optimiser can't step into invalid values; transform back and use the delta method for its standard error.
Good starting values (for example OLS estimates) prevent many convergence failures.
Use scipy.stats log-densities (stats.norm.logpdf, stats.t.logpdf) to write likelihoods accurately.
Name parameters with exog_names or the extra_params_names argument so the summary is readable.

The same machinery underlies the built-in discrete and count models.

Example

import numpy as np
from scipy import stats
from statsmodels.base.model import GenericLikelihoodModel

class TRegression(GenericLikelihoodModel):
    """Linear regression with Student-t errors (robust to outliers)."""
    def loglikeobs(self, params):
        beta, log_scale, log_df = params[:-2], params[-2], params[-1]
        resid = self.endog - self.exog @ beta
        return stats.t.logpdf(resid, df=np.exp(log_df), scale=np.exp(log_scale))

rng = np.random.default_rng(14)
x = rng.uniform(0, 10, 300)
y = 1 + 2 * x + stats.t.rvs(df=3, size=300, random_state=rng)
X = np.column_stack([np.ones(300), x])

mod = TRegression(y, X, extra_params_names=["log_scale", "log_df"])
res = mod.fit(start_params=[0, 1, 0, 1], method="bfgs", disp=0)
print(res.params.round(2))            # [0.91 1.99 0.   1.15]
print(res.bse[:2].round(3))           # [0.139 0.024]
print(round(np.exp(res.params[-1]), 1))  # 3.2 (estimated degrees of freedom)

Remember

Implement loglikeobs returning per-observation log-likelihoods
Reparameterise constrained parameters (log scale) and give good start values
You get standard errors, tests and robust covariances for free

Watch out for

Returning the total log-likelihood from loglikeobs instead of one value per observation.

Check yourself

Why estimate log(sigma) instead of sigma in a custom likelihood?

Official docs: www.statsmodels.org/stable/examples/notebooks/generated/generic_mle.html

Example

Remember

Watch out for

Check yourself

Builds on

Unlocks