Negative covariance between goals in bivariate Poisson model I'd like to approximate the covariance between goals in football dataset assuming they come from bivariate Poisson distribution. In previous works, which model the distribution of goals by bivariate Poisson process (i.e. http://tolstoy.newcastle.edu.au/R/e8/help/att-6544/karlisntzuofras03.pdf by Karlis and Ntzoufras or https://pdfs.semanticscholar.org/fca7/0c6bd99b082759c08c035a8ecc6a2cac15de.pdf by Koopman and Lit) the authors state, that the covariance is positive and equal approximately 0.1.
However, when I use the sample covariance estimator function on goal counts, I obtain negative values (-0.05 to -0.18, depending on the football league and time period). I am using the same dataset as Koopman and Lit in above paper.
I'm using the following formula for sample covariance:
$cov(X,Y)= \frac{1}{n-1}\sum (X_i-\overline{X})(Y_i-\overline{Y})$
Where do the differences come from? 
 A: In my own opinion, the best model is offered by creating a bivariate Poisson distribution where the likelihood of the teams scoring goals is directly proportional to the teams' rankings inflated by a power of a number, which will be calibrated based on real results.
In this case I have two teams, A and B. To keep things simple they play at a neutral venue. I can talk later about home advantage but it is a different issue. Team A has a ranking of 1600 points and Team B has a ranking of 800 points.
Now: Team A will have: An expected number of goals in the match dependent on the average number of goals scored in a match, divided by 2. In this case, I will inflate its ranking by one observation from the normally distributed inflation factors that I get from one random observation: Let us say in this case that the value I get is 100%. If average goals per match = 2,8, then: Lambda of Team A will be: 2,8*1600/(1600+800)=1,87 goals, Team B will get the rest, i.e. 0,93 goals.
I create a 10x10 table, and I insert all scores, based on: Score (x:y) follows B.P. (Lambda1, Lambda2) where Lambda1=1,87, Lambda2=0,93.
Then, I can estimate the probability of a draw, home or away win by just counting the scores corresponding to the specific result. 
I have tested my model and it performs better than most, the mean of the normal distribution for inflation would be around 185,8% and the s.d. would be around 110%... 
I can send you a spreadsheet you can see for yourself... 
