Massey ratings computation
in SageMath

This converted Jupyter notebook shows how to do the Massey rating computations for the first basic example shown in this Observable notebook but using SageMath to do the computations.

Even if you don't have SageMath, you can find a live version of these computation on the SageCell Server.

The setup

Recall that we have three teams (let's call them $A$, $B$, and $C$) that play each other twice each. We can represent the games using a directed graph:

In this image, each edge corresponds to a game. There is an edge $e_k$ from team $i$ to team $j$ if team $i$ defeated team $j$ in game $k$ and the edge is labeled with margin of victory.

Written out as a table we might have

game	winner	loser	margin
0	A	B	2
1	A	B	3
2	A	C	3
3	B	C	1
4	C	A	1
5	C	B	2

We can record these games in a win/loss matrix $X=(x_{kj})$ together with a vector $\vec{y}=(y_k)$ of margins of victory. Recall that the rows of $X$ are indexed by the games, the columns are indexed by the teams, and $$x_{kj} =\begin{cases} 1 & \text{if team } j \text{ won game } k \\ -1 & \text{if team } j \text{ lost game } k \\ 0 & \text{ if team } j \text{ did not play in game } k. \end{cases}$$

In addition, $y_k$ denotes the margin of victory in game $k$.

Written in code, the matrix $M$ and vector $\vec{y}$ look like so:

X = Matrix([
    [1,-1,0],
    [1,-1,0],
    [1,0,-1],
    [0,1,-1],
    [-1,0,1],
    [0,-1,1]
])
y = Matrix([[2,3,3,1,1,2]]).transpose()

We are interested in a least squares solution of $$M\vec{r} \approx \vec{y},$$ for then, $$r_i - r_j = y_k.$$ In other words, the difference in the ratings of two teams is indicative of the score difference. To find this solution, we write

$$M = X^T X \text{ and } \vec{p} = X^T \vec{y}$$

and then solve $$M\vec{r} = \vec{p}.$$

In code, we can find $M$ and $\vec{r}$ as follows:

XT = X.transpose()
M = XT*X
n = M.dimensions()[0]
for i in range(n):
    M[n-1,i] = 1
p = XT*y
p[n-1,0] = 0
[M,p]

[
[ 4 -2 -2]  [ 7]
[-2  4 -2]  [-6]
[ 1  1  1], [ 0]
]

We can then solve the system:

M.solve_right(p)

[ 7/6]
[  -1]
[-1/6]

This indicates that team $A$ is the strongest, followed by team $C$ and finally by team $B$.