How do I fit a double Poisson regression model? $$
X_k \sim \mathcal{Pois}(\lambda_k) \\
Y_k \sim \mathcal{Pois}(\mu_k) \\
\\
\ln \lambda_k = \alpha_{i(k)} + \beta_{j(k)} + \gamma + \eta / 2 \\
\ln \mu_k = \alpha_{j(k)} + \beta_{i(k)} + \gamma - \eta / 2
$$
Imagine the above problem for modelling basketball points,
Where $X_k$ represents the total points scored by the home team in game $k$ and $Y_k$ represents the total points scored by the away team in game $k$ (assume $X$ and $Y$ are independent for simplicity).
How would one fit these models using Maximum-likelihood estimation?
$(\alpha_1,...,\alpha_n)$ and $(\beta_1,...,\beta_n)$ are two vectors of coefficients to be estimated.
$\alpha_{i(k)}$ is an 'attacking' coefficient for team $i$ (the 'home' team) in game $k$.
$\beta_{j(k)}$ the 'defending' coefficient for team $j$ (the 'away' team) in game $k$.
$\gamma$ is a constant, could be interpreted as the average points scored per game.
$\eta$: I'm not too sure, possibly a binary indicator variable or maybe another constant. Something like home advantage which could take the form $(X-Y)$ but I think that would violate the model ?

You can fit this as an ordinary Poisson glm with a log link function, treating each $X_k$ and $Y_k$ scored in each game as separate observations, including the attacking and defending teams as factors, and including the home-/away-field advantage/disadvantage via a binary factor with a sum-to-zero contrast.  If you have $m$ teams, you need to impose some constraints on the $m$ $\alpha_i$'s and $\beta_i$'s to make the model identifiable, for example, via treatment or sum-to-zero contrasts.
