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 ?

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
