Signed differences of counts are common yet challenging to model while preserving both parsimony and interpretability. The Skellam distribution is the classical choice, and Skellam regression remains standard for score differences; zero–inflated variants and ad hoc mixtures partly address excess ties and asymmetric tails. We develop a Beta–Gamma Skellam (BG–Skellam) regression that reparameterizes the competing Poisson means via a total intensity \(\lambda \) and a balance fraction \(\xi \) , assigns log/logit links to \(\mathbb {E}[\lambda \mid x]\) and \(\mathbb {E}[\xi \mid x]\) , and induces flexible, finite–variance Skellam mixtures through Gamma/Beta compounding. The construction nests Skellam regression and zero–inflated Skellam as limits and is identifiable under standard analytic–transform conditions. Estimation proceeds by maximum likelihood, evaluating the BG–Skellam kernel stably via pgf inversion with short quadrature; for completeness, we also record in an appendix a Bayesian formulation with Bessel augmentation as a computational roadmap. As an illustration, we analyze English Premier League full–time goal differentials: covariates summarizing recent attacking/defensive form and overall scoring intensity improve held–out log predictive density and better calibrate ties and small–margin outcomes than Skellam and zero–inflated Skellam benchmarks, while retaining a transparent exposure–balance interpretation through \((\lambda ,\xi )\) and interpretable covariate effects on \(\mathbb {E}[\lambda \mid x]\) and \(\mathbb {E}[\xi \mid x]\) . BG–Skellam offers a unified, regression–ready framework for signed discrete outcomes.
