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 ?
 A: 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.
A: It sounds like you could use bivariate Poisson distribution (Karlis and Ntzoufras, 2003) to model this data, especially since it was already used in similar context. The distribution is defined in terms of three random variables
$$
Z_i \sim \mathcal{P}(\lambda_i), \quad i = 1,2,3
$$
where
$$\begin{align}
X &= Z_1 + Z_3 \\
Y &= Z_2 + Z_3 \\
(X, Y) &\sim \mathcal{BP}(\lambda_1, \lambda_2, \lambda_3)
\end{align}$$
with probability mass function
$$
f(x) = \exp \left\{-(a+b+c)\right\} \frac{a^x}{x!} \frac{b^y}{y!} \sum_{k=0}^{\min(x,y)}{x \choose k} {y \choose k} k! \left( \frac{c}{ab} \right)^k
$$
In such a case, you would model $\lambda_1$ and $\lambda_2$ parameters in terms of your $\alpha$ and $\beta$ parameters and $\lambda_3$ as the joint term. $\lambda_3$ would influence the degree of dependence between $X$ and $Y$. So there would be some slight modifications of your model.
There is an R package for bivariate Poisson GLMs, PMF and random generator are available through R's extraDistr package (disclosure: I'm the author), or you could define the model in a probabilistic programming language of your choice.

Karlis, D. and Ntzoufras, I. (2003). Analysis of sports data by using bivariate Poisson models.
Journal of the Royal Statistical Society: Series D (The Statistician), 52(3), 381-393
