Simulating election outcomes from betting market data

Because betting market odds can be expressed as probabilities, they can be used in Monte Carlo simulations to reflect the election outcome that punters think most likely. However, there are a couple of technical issues that need to be considered. One is simple; the other, less so.

Bookmakers make their money by ensuring that whomever wins, their outlays will always be less than their takings. A bookmaker adjusts her odds as bets are laid to maintain a margin for herself. Fair enough, bookmakers need to cover their costs and make a profit. This is variously known as the bookmaker's overround, margin, vigorish, vig or juice. To get fair probabilities from the bookmaker's odds, we need to correct for the bookmaker's margin. 

The more challenging technical issue is what is known as the Favourite-Longshot Bias (FLB). First reported by R M Griffith in the American Journal of Psychology in 1949, FLB notes that on average short odds under-estimate the probability of winning and that long odds over estimate this probability. As a consequence, bets placed on higher-valued odds (the longshots) offer worse rates of return than bets placed on lower-valued odds (favourites). This tendency in betting markets has been empirically validated many times.

It is evident in our data. In almost every seat, sportsbet offers odds of \$101 on the United Australia Party (UAP) winning the seat. This is an individual seat win probability of 0.0099009900990. In a 151 seat Parliament, the most likely outcome for these same odds across every seat is the UAP winning one or two seats in the next election. Sorry, but I would be surprised if the UAP won any of these seats, with those odds. Consistent with FLB, these odds of \$101 overstate the UAPs probability of winning.

If I don't adjust the bookmaker's odds for FLB, the Monte Carlo simulations ends up with an implausibly large number of Greens, independents and other minor parties winning seats in the House of Representatives, as evident in the following charts.

So what adjustments should I make? One option is to arbitrarily ignore longshot odds over a certain value (say \$20). This is the approach I took for the 2019 election.

While this yields a more plausible simulation, the cut-off is arbitrary. Another option is to transform the odds so that the higher valued odds are made even higher to correct for FLB. I have considered two transformations. First multiplying the raw-odds by their square root (this is the same as raising the raw-odds to the power of 1.5). Second, squaring the raw-odds (which is the same as raising the raw-odds to the power of 2).

The next set of charts come from the Monte Carlo simulation where the raw odds were multiplied by their square root. The results are similar to the charts above where longshot odds over \$20 were ignored.

The final treatment I considered was squaring the raw-odds, before converting them to probabilities and standardizing these probabilities so that they sum to one in each seat.

I am still deciding which transformation of the raw-odds is best. While the squared approach yields probabilities for the number of seats won by the Greens and Others more in line with the current parliament, the variance of the probability distributions for all four groups has been reduced. The bias-variance trade-off  suggests that we need to at least consider the possibility that we might be over-fitting the betting market data. Nonetheless, for the moment, I plan to use this squared-odds approach for managing the favourite-longshot bias (FLB) within the Monte Carlo simulations. Otherwise, the number of seats won by the Greens and others still seems too high. But I will review this approach from time to time.

For the individual seat outcome probabilities, I will take the less aggressive approach of multiplying the raw-odds by their square root before calculating the probabilities and adjusting for the bookmaker's margin.

If you have alternate/better treatments for FLB, please drop a note in the comments below and argue your case.

Caveat

A couple of people have noted that I am using independent draws in the above Monte Carlo simulation. They argued, because voting across electorates is correlated, I should have used dependent draws, or correlated draws. It is a fair points, but this requires substantial work to understand and model the dependency structure. I have added this to my to do list. One of my interlocutors provided a couple of useful links:

• From Andrew Gelman's blog: https://statmodeling.stat.columbia.edu/2020/10/24/reverse-engineering-the-problematic-tail-behavior-of-the-fivethirtyeight-presidential-election-forecast/
• From 538: https://fivethirtyeight.com/features/election-update-why-our-model-is-more-bullish-than-others-on-trump/
  

Monte Carlo simulation

At the start of the campaign, I ran a Monte Carlo simulation of the election. I thought it would be interesting to run it again with the latest Newspoll data. This simulation is based on a 49 per cent two-party preferred (TPP) vote share for the Coalition.

But first some caveats. My Monte Carlo model is very simple. It does not make adjustments for the sophomore effect (a second election bounce for first-term members running for re-election), nor the retirement effect (a loss of party support when a long-term member retires). The model also makes use of an experimental application of the Dirichlet distribution to randomise state swings within a constrained national swing.

The most likely result is Labor would win 80 seats and the Coalition 65. I have allocated six seats to others (two from Labor and four from the Coalition). This is much closer than it was at the start of the election campaign (Labor 87 to Coalition 58) and reflects a significant tightening of from Newspoll over the last three months.

The full set of probability outcomes for Labor and the Coalition from the Monte Carlo simulation is as follows. There is some possibility of a Labor minority government. Much of the distribution comes from the risk of a polling failure.

In terms of the Cube Rule, this suggests Labor would only win 78 seats from a national TPP vote for Labor of 51 per cent. Again I have allocated six seats to others.

My adjusted state swings were as follows. These are the swings to the Coalition in percentage points for each state from that state's 2016 TPP vote share for the Coalition.

The seat-by-seat TPP results follow. This includes seats which would be won by the other party on a two-candidate preferred basis.

Monte Carlo simulation of the 2019 Australian election

In September last year, I cobbled together a simple Monte Carlo simulation of the Federal election using:

• 2016 seat-by-seat two-party preferred (TPP) outcome data,  
• an estimate of state-by-state TPP polling outcomes from the Newspoll Quarterly,
• an estimate of the national TPP vote from my latest Bayesian poll aggregation

Back then, when the Coalition's TPP vote share was 45 per cent in the opinion polls, I estimated the most likely election outcome was a Labor win with 95 seats. I thought the Coalition would get 51 seats, and 5 seats would go to others. I noted that my Monte Carlo model did not account for the retirement effect nor the sophomore effect.

Today, I re-ran that model with updated inputs:

• the latest Coalition TPP vote share estimate (47.4 per cent) 
• the number of enrolled voters by state
• the estimate of TPP vote share at the 2016 election, applied to the new electoral boundaries
• the latest state-by-state TPP estimates from the Newspoll Quarterly
• I added Wentworth to the seats that would be held by a minor party or independent (bringing the total to six). This seat would come from the Coalition's total.

The headline result was a Labor win with 87 seats. The Coalition got 58 seats in the simulation. And 6 seats were allocated to minor parties and independents.

This compares reasonably well with the Cube Rule, which I use as a back-of-the-envelope validation check.

My results are reasonably comparable with Antony Green's election calculator (although I give three more seats to Labor than Antony). For a Labor TPP vote share of 52.6 per cent, Antony's calculator estimates Labor would win 84 seats, and the Coalition would win 61 seats. Like me, Antony allocates six seats to independents and other parties. I think the key difference is that my model automatically updates for state swings; whereas the way in which I used Antony's model was a single national swing.

My adjusted state swings have Queensland, WA and NSW as the key trouble spots for the Coalition. Victoria looks the least profitable for Labor. Note: I used the national swing for Tasmania, ACT and the NT, because I did not have any other data. (Note: my model adjusts the state swings from the Newspoll Quarterly such that the state swings are always consistent with the national swing. It does this in a way that allows some variability across simulation runs/ This variability can be seen in the next chart which provides a kernel density estimate for each state swing across the 100,000 simulation runs.)

Based on the two-party preferred data, I plot a two-party probability outcome for each seat. This is an artificial construct, as a number of seats were a two-candidate preferred outcome at the 2016 election where one of the candidates did not come from the major parties. This table includes the six seats that I have allocated to others in the totals above.

Update

I have quite different state-by-state TPP estimates compared with William Bowe from the Poll Bludger. I am looking into these differences, and may need to update this page when I better understand the source of Bowe's estimates.

Monte Carlo simulation of elections

Between elections, the Australian Election Commission (AEC) redraws the electoral boundaries to ensure each seat has a similar number of voters. Now that this redistribution process has been completed, I can use the new seats to model election outcomes.

The first thing I needed was the recalculated margins for each seat. For this data I used Wikipedia. Antony Green has also undertaken these calculations. For the seats that had not been redistributed, we have original polling outcome data from the AEC. This base, expressed as margins, looks something like this.

With this base, I have built a Monte Carlo simulation. In a Monte Carlo simulation we sample from probability distributions many thousands of times to identify the range of possible outcomes. These are then analysed to identify the probabilities for different events occurring. The model needs to consider those factors that can see the results vary.

The biggest source of uncertainty I need to manage is polling uncertainty. It is not unusual for an aggregated opinion poll to be plus or minus two percentage points from the final election outcome. In the Monte Carlo model I have assumed that the actual election outcome will be normally distributed around the poll estimate with a standard deviation of one percentage point.

Another source of uncertainty is the way in which the swing in the individual seats is distributed around the national swing and the way in which this swing varies state-by-state. Historically, individual seat swings have been close to normally distributed around the the national swing with a standard deviation of 3 percentage points. They have also been close to normally distributed around state swings with a standard deviation of 2.5 percentage points.

For this model, I have used state swings, based on the most recent state-by-state Newspoll (which pre-dates the Morrison ascendancy), and then adjusted for a change in the aggregate two-party preferred (TPP) since the Morrison ascendancy. I draw random numbers from a Dirichlet distribution to achieve this adjustment. The state swings since the last election the model used can be seen in the following kernel density estimate plot. The chart is of the state swings to the Coalition in percentage points since the 2016 election. The largest swings against the government at the moment appear to be in WA and Queensland.
 

There are two key factors that I have not modeled. The first is the sophomore effect - a bump that first term members of Parliament get when running for re-election. The second factor is the retirement effect - a decline in the party vote in a seat following the retirement of long standing member for that party. Labor has a large number of first-term parliamentarians, and is likely to benefit from the sophomore effect at the next election, not withstanding it also has a number of retirees.

A further (and perhaps more critical) factor I have not modeled is the outcome in seats currently held by other parties. For this analysis I have simply assumed those seats will continue to be held by other parties.

In the current climate, with an estimated aggregate TPP of 45 per cent for the Coalition. The model predicts a substantial victory for Labor were an election held now. Based on a simulation run of 100,000, the model predicts Labor is most likely to win 95 seats, and the Coalition 51 seats.

While this is the most likely outcome, there are a cluster of possible outcomes for both parties. But there is little doubt, if an election were held now, a significant Labor majority would be the outcome.

Turning to the individual seat outcomes, these are charted below. In this chart, the seats where we have the Coalition at zero or 100 per cent probability are not sorted.

And finally, my rough and ready code for this exercise. Usual caveats apply: this is a work in progress.

# PYTHON: Monte-Carlo simulation of election outcomes
# -- NOTE: a number of data sources need to be updated 
#          in this code before it is run.

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import sys
sys.path.append( '../bin' )
plt.style.use('../bin/markgraph.mplstyle')


# --- version information
print('Python version: {}'.format(sys.version))


# --- Seat data
# Seat data sourced from 
# https://en.wikipedia.org/wiki/Pre-election_pendulum_for_the_next_Australian_federal_election
workbook = pd.ExcelFile('./Data/Seats.xlsx')
df = workbook.parse('seats')
df.index = df.Seat

Coalition_TPP_2016 = 0.5036
# ===> UPDATE HERE <===
Coalition_TPP_now  = 0.4500 
# TO DO - source Coalition_TPP_now directly from TPP aggregation
Swing_to_Coalition = Coalition_TPP_now - Coalition_TPP_2016

others = df[df['LNP TPP'].isnull()]
df = df[df['LNP TPP'].notnull()]
base = ((df['LNP TPP'] / 100.0) - 0.5) 
# NOTE: base < 0 is Labor; base > 0 is Coalition 


# --- State Data - note: TPP from the Coalition's perspective. 
states =     ['NSW',   'Vic',   'Qld',   'WA',    'SA',    'Tas',   'ACT',   'NT']

# voters from https://www.aec.gov.au/Enrolling_to_vote/Enrolment_stats/national/index.htm
# ===> UPDATE HERE <===
voters =     [5211182, 4094212, 3203789, 1615900, 1200395,  381409,  290654, 138581]
voters = pd.Series(voters, index=states)

# State 2016 TPP source: https://results.aec.gov.au/20499/Website/HouseTppByState-20499.htm
tpp_2016 =   [0.5053,  0.4817,  0.5410,  0.5466,  0.4773,  0.4264,  0.3887,  0.4294]

# latest TPP estimate draws on 
# https://www.theaustralian.com.au/national-affairs/turnbull-axed-as-coalition-closed-the-gap-on-labor/news-story/487dd05cd4dc95693bd6c55b44bfbe88
# ===> UPDATE HERE <===
tpp_est_now= [0.4963,  0.4597,  0.5000,  0.4996,  0.5103, 0.4164,   0.3787,  0.4194]

# the multinomial vector for drawing the Dirichlet random numbers of state swings
alpha_scale = 10000000 
state_alpha = (tpp_est_now * voters * alpha_scale / voters.sum()).astype(int) 


# --- let's simulate ...
Monte_Carlo_N = 100000
# next line - preallocate space to speed up calculations
simulations = pd.DataFrame(np.zeros((len(base),Monte_Carlo_N)))
simulations.index = df.index
state_swings = pd.DataFrame(np.zeros((len(states),Monte_Carlo_N)))
state_swings.index = states

print('Commencing ', str(Monte_Carlo_N), ' simulation run ...')
for i in range(Monte_Carlo_N) :

    # -- progress indication
    if i % (Monte_Carlo_N // 20) == 0 :
        print(i)

    # -- polling uncertainty - polls often out by up +/- two percentage points
    pollingUncertainty = np.random.standard_normal(1) * 0.01 # = standard deviation
    #pollingUncertainty = 0.0
    
    # -- variable swing by state - use a dirichlet random to manage this element to ensure
    #    total Coalition vote is the same as the Coalition TPP for all eligible voters
    # NOTE: Comment out this section to use national swings rather than state swings
    # NOTE: drawing random numbers from the Dirichlet distribution is slow
    state_dirichlet = np.random.dirichlet(state_alpha) # proportion of Coalition cote in each state
    state_Coalition_tpp = state_dirichlet * (voters.sum() * Coalition_TPP_now) / voters
    state_swing_to_coalition = state_Coalition_tpp - tpp_2016
    state_swings[i] = state_swing_to_coalition # we will plot this
    Swing_to_Coalition = df.State.map(state_swing_to_coalition)
    
    # -- TO DO - retirement effect
    
    # -- TO DO - sophomore effect
    
    # -- variable swing seat-by-seat - normally distributed noise around 0
    # -- use a standard deviation of 0.03 for national swings
    # -- use a standard deviation of 0.025 for state swings
    # -- https://marktheballot.blogspot.com/2016/11/how-are-seat-swings-distributed-around.html
    # -- https://marktheballot.blogspot.com/2012/11/state-swings.html
    seatDistributedAroundSwing = np.random.standard_normal(len(base)) * 0.025 # = standard deviation
    
    # -- bring it all together ...
    simulations[i] = base + Swing_to_Coalition + pollingUncertainty + seatDistributedAroundSwing
    

print('Finished simulation ... analysing data ...')
sumCoalition = simulations[simulations >= 0].count()
sumLabor = len(base) - sumCoalition
simSummary = pd.concat({'Coalition': sumCoalition.value_counts(), 
    'Labor': sumLabor.value_counts()}, axis=1)
min_value = simSummary.index.min()
max_value = simSummary.index.max()+1
simSummary = pd.DataFrame(simSummary[['Labor', 'Coalition']], 
    index=range(min_value, max_value))
simSummary = simSummary / simSummary.sum() 
simSummary = simSummary.sort_index()


# -- seat count distributional plot
print('About to plot ...')
ax = simSummary.plot.bar()
ax.set_title('Election Outcome Probabilities for Coalition TPP: ' 
    + str(Coalition_TPP_now * 100.0))
ax.set_xlabel('Seats Won')
ax.set_ylabel('Probability') 
ticks = np.arange(min_value, max_value, 5)
ax.set(xticks=[x - ticks[0] for x in ticks], xticklabels=ticks)

fig = ax.figure
fig.set_size_inches(8, 4)
fig.tight_layout(pad=1)
fig.text(0.99, 0.01, 'marktheballot.blogspot.com.au',
        ha='right', va='bottom', fontsize='x-small', 
        fontstyle='italic', color='#999999') 
fig.savefig('./Graphs/SeatCountProbabilities.png', dpi=125) 
plt.close()


# -- most likely outcome plot
ml_coalition = simSummary[simSummary['Coalition'] == 
    simSummary['Coalition'].max()].index[0] 
ml_labor = len(base) - ml_coalition
ml_other = len(others)
ml_outcome = pd.Series(data=[ml_labor, ml_coalition, ml_other], 
    index=['Labor', 'Coalition', 'Other'])
    
ax = ml_outcome.plot.barh()
ax.set_title('Most likely Election Outcome for Coalition TPP: ' 
    + str(Coalition_TPP_now * 100.0))
ax.set_xlabel('Number of Seats Won by Party')
ax.set_ylabel('') 

for i in ax.patches:
    ax.text(x=1, y=i.get_y()+.16, s=str(i.get_width()), color='white')
    
fig = ax.figure
fig.set_size_inches(8, 4)
fig.tight_layout(pad=1)
fig.text(0.99, 0.01, 'marktheballot.blogspot.com.au',
        ha='right', va='bottom', fontsize='x-small', 
        fontstyle='italic', color='#999999') 
fig.savefig('./Graphs/SeatLikelyOutcome.png', dpi=125) 
plt.close()


# -- plot state swings to the Coalition - as a KDE
state_swings = state_swings * 100 # covert to percentage points
ax = state_swings.T.plot.kde()
ax.set_title('State Swing Kernel Density Estimates for Coalition TPP: ' 
    + str(Coalition_TPP_now * 100.0))
ax.set_xlabel('Swing to the Coalition in Percentage Points')
ax.set_ylabel('Density') 

fig = ax.figure
fig.set_size_inches(8, 4)
fig.tight_layout(pad=1)
fig.text(0.99, 0.01, 'marktheballot.blogspot.com.au',
        ha='right', va='bottom', fontsize='x-small', 
        fontstyle='italic', color='#999999') 
fig.savefig('./Graphs/Seat-StateSwingKDE.png', dpi=125) 
plt.close()


# -- plot base
base = base * 100 # covert to percentage points
base.sort_values(inplace=True)
ax = base.plot.barh(color='royalblue')
ax.set_title('2016 Coalition Margins Starting Point')
ax.set_xlabel('Percentage Points (Labor is <0; Coalition is >0)')
ax.set_ylabel('Seat') 

fig = ax.figure
fig.set_size_inches(8, 30)
fig.tight_layout(pad=1)
fig.text(0.99, 0.005, 'marktheballot.blogspot.com.au',
        ha='right', va='bottom', fontsize='x-small', 
        fontstyle='italic', color='#999999') 
fig.savefig('./Graphs/Seat-baseMargins.png', dpi=125) 
plt.close()


# -- plot individual seat outcomes
sumSeatCoalition = simulations[simulations >= 0].count(axis=1) 
sumSeatCoalition = sumSeatCoalition / Monte_Carlo_N
sumSeatLabor = 1.0 - sumSeatCoalition
seatSummary = pd.DataFrame(data={'Coalition': sumSeatCoalition, 'Labor': sumSeatLabor})
seatSummary = seatSummary[['Labor', 'Coalition']] # correct order for colours
seatSummary.sort_index(inplace=True)

ax = seatSummary.plot.barh(stacked=True, legend=False)
ax.set_title('Seat Win Probabilities for Coalition TPP: ' 
    + str(Coalition_TPP_now * 100.0))
ax.set_xlabel('Probability')
ax.set_ylabel('') 

fig = ax.figure
fig.set_size_inches(8, 30)
fig.tight_layout(pad=1)
fig.text(0.99, 0.005, 'marktheballot.blogspot.com.au',
        ha='right', va='bottom', fontsize='x-small', 
        fontstyle='italic', color='#999999') 
fig.savefig('./Graphs/SeatWinProbabilitiesNameOrder.png', dpi=125) 
plt.close()

seatSummary.sort_values(by='Coalition', inplace=True)
ax = seatSummary.plot.barh(stacked=True, legend=False)
ax.set_title('Seat Win Probabilities for Coalition TPP: ' 
    + str(Coalition_TPP_now * 100.0))
ax.set_xlabel('Probability')
ax.set_ylabel('') 

fig = ax.figure
fig.set_size_inches(8, 30)
fig.tight_layout(pad=1)
fig.text(0.99, 0.01, 'marktheballot.blogspot.com.au',
        ha='right', va='bottom', fontsize='x-small', 
        fontstyle='italic', color='#999999') 
fig.savefig('./Graphs/SeatWinProbabilitiesOutcomeOrder.png', dpi=125) 
plt.close()