Data and code for election 2016

I have made most of my data and code base available on Google Drive. Please note, this is my live code base, which I play with quite a bit. So, there will be times when it is broken or in some stage of being edited.

What I have not made available is the Excel spreadsheets into which I initially place my data. These live in the (hidden) raw-data directory. However, the collated data for the Bayesian model lives in the intermediate directory, visible from the above link.

The program that collates and organises the spreadsheets into a single CSV input file for the Bayesian analysis is TPP-step1.py. There are two intermediate input files (at the moment): TPP-3-stage1.csv and TPP-all-stage1.csv. The first of these input files covers the most recent three months. The second of these input files is for all polls since the 2013 election.

The Bayesian model itself lives in the file TPP-step2.R.

The code for producing the plots lives in TPP-step3.py.

The files that begin with the letter 'z' are bash shell scripts.

The most recent set of charts live in the graphs directory. I don't keep historical charts.

There are a handful of helper programs that live in the bin directory.

There are a couple files that I am working on in respect of a primary votes model. This is still a long way from finished.

If you see an error in my code or data, please drop me a line (comments below or email address in right hand column), and let me know. I can only improve with your help.

Further explorations in non-linearity

On the weekend I began exploring a non-linear model for aggregating polling. My test case produced nice looking graphs; but the results were in some large part an artifact of the priors I had chosen for the model (not good).

In the comments to that post, I suggested that part of the problem may have been how I had defined the model.

Thinking about this some more, I think the problem was in this statement: << beta-for-pollster * (poll-result - minimum-poll-result) >>. Rather than the minimum, it should be a central tendency of some type (mid-point, mean, median, etc.).

I have now redefined the model to use the midpoint for each polling house (defined as the (min+max)/2 for that house). The use of house-specific mid-points acknowledges that the raw x scores include the alpha house effect I am trying to estimate. The revised model is:

house-effect-for-pollster = alpha-for-pollster + beta-for-pollster * (poll-result - pollster-mid-point)

The priors to the beta effect have been completely freed-up, so that they are uninformative. As a result, this model works better internally than the previous model.

However, I have added a constraint such that the aggregation assumes the "fizziness" factors also sum to zero. This may be a bollocks assumption; and will need some further analysis.

So - with the caveat that this is still very early exploratory analysis - the results of this new model follow.

On a plain reading, the above charts suggest that the current Rudd effect may actually be better for Labor than other aggregations have found. However, I am not convinced this is anything more than an artifact of the model. I am particularly concerned about the inclusion of data from polling houses that do not have data points across the highs and lows of the Gillard period, and across the Rudd and Gillard periods. Such data don't fully occupy the model, which can result in artifacts.

If we limit the analysis to data from Nielsen and Newspoll for a comparison using data that does fully occupy the model, we get the following.

Because I am committed to presenting poll analysis in a transparent and unbiased fashion, the next two code snippets show how I package up the data for the model (in R), and the model itself (in JAGS).

    # - prepare the data ...
    data$y <- data[ , y] / 100 # vote share between 0-1
    data$variance <- data$y * (1 - data$y) / data$Sample #pq/n
    data$standardError <- sqrt(data$variance)
    data$samplePrecision <- 1 / data$variance
    HOUSES <- levels(data$House)
    HOUSECOUNT <- length(levels(data$House))
    NUMPOLLS <- length(data$y)
    cat('Number of polls: '); print(NUMPOLLS)
    midpoints <- rep(NA, HOUSECOUNT)
    for(i in seq_len(HOUSECOUNT)) {
     midpoints[i] <- (max(data[data$House==HOUSES[i], 'y']) + 
      min(data[data$House==HOUSES[i], 'y'])) / 2
     cat('Midpoints: '); cat(i); cat(' '); cat(HOUSES[i]); cat(': '); print(midpoints[i])
    }
    
    # - remember these dots 
    dotsFile <- paste('./files/', fPrefix, 'original.csv', sep='')
    write.csv(data, file=dotsFile)

    # - manage dates
    day0 <- min(data$Date) - 1  # walk starts from earliest date
    if(endDate == TODAY)
        endDate <- max(data$Date)
    PERIOD <- as.numeric(endDate - day0) # length of walk in days
    DISCOUNTINUITYDAY <- as.numeric(as.Date(discontinuity) - day0)
    cat('Discontinuity: '); print(DISCOUNTINUITYDAY)
    tPrefix <- paste(format(day0+1, '%d-%b-%Y'), ' to ', format(endDate, '%d-%b-%Y'), sep='')
    data$day <- as.numeric(data$Date - day0)

    # - do the MCMC thing ...
    
    parameters <- list(PERIOD = PERIOD,
        HOUSECOUNT = HOUSECOUNT,
        NUMPOLLS = NUMPOLLS,
        DISCOUNTINUITYDAY = DISCOUNTINUITYDAY, 
        NEWSPOLL = which(levels(data$House) == 'Newspoll'),
        y = data$y,
        x = data$y,
        day = data$day,        house = as.integer(data$House),
        samplePrecision = data$samplePrecision,
        midpoints=midpoints
    )
                
    jags <- jags.model(textConnection(model),
        data=parameters,
        n.chains=4,
        n.adapt=n.adapt
    )

    # - burn in
    update(jags, n.iter=n.update) # burn-in the chains

    # - capture results
    jags.capture <- c('walk', 'alpha', 'beta', 'discontinuityValue')
    coda.samples <- coda.samples(jags, jags.capture, n.iter=n.iter, thin=n.thin)
    coda.matrix <- as.matrix(coda.samples)

model {
    ## Derived from Simon Jackman's original model 
    
    ## -- observational model
    for(poll in 1:NUMPOLLS) { 
        # note: x and y are the original polling series
        houseEffect[poll] <- alpha[house[poll]] + 
         beta[house[poll]]*(x[poll]-midpoints[house[poll]]) 
        mu[poll] <- walk[day[poll]] + houseEffect[poll]
        y[poll] ~ dnorm(mu[poll], samplePrecision[poll]) 
    }
            
    ## -- temporal model
    for(i in 2:PERIOD) { # for each day under analysis ...
        day2DayAdj[i] <- ifelse(i==DISCOUNTINUITYDAY, 
            walk[i-1]+discontinuityValue, walk[i-1])
        walk[i] ~ dnorm(day2DayAdj[i], walkPrecision)
    }
    sigmaWalk ~ dunif(0, 0.01)            ## uniform prior on std. dev.  
    walkPrecision <- pow(sigmaWalk, -2)   ##   for the day-to-day random walk
    walk[1] ~ dunif(0.4, 0.6)             ## uninformative prior
    discontinuityValue ~ dunif(-0.2, 0.2) ## uninformative prior

    ## -- house effects model
    for(i in 2:HOUSECOUNT) { ## vague normal priors for house effects
        alpha[i] ~ dunif(-0.1,0.1)
        beta[i] ~ dunif(-1,1) 
    }
    alpha[NEWSPOLL] <- -sum(alpha[2:HOUSECOUNT]) ## sum to zero
    beta[NEWSPOLL] <- -sum(beta[2:HOUSECOUNT]) ## sum to zero
}

How much was Kevin Rudd worth?

I was a little surprised when I saw Simon Jackman suggest that Kevin Rudd had moved the two-party preferred voting intention by seven percentage points in Labor's favour. It was not consistent with my own analysis and only one pollster (Morgan) has data that supports a seven point movement. Data from all the remaining pollsters suggest the "Rudd Effect" was less than seven percentage points.

Now don't get me wrong, I have enormous respect for Professor Jackman. I purchased and read his 600 page text, Bayesian Analysis for Social Sciences. It is a tour de force on Bayesian statistics. I cannot recommend this book enough. His understanding and knowledge in this area far surpasses my own. Unashamedly, I have used Jackman's approach as the basis for my own aggregation efforts.

However, I suspect he has not noticed that the data since the second ascension of Keven Rudd violates a number of the linear assumptions implicit in his model. In particular, some of the house effects before and after Kevin are radically different. I blogged on this under the rubric: When models fail us. As I noted previously, the violation of the underpinning assumptions results in the model producing incorrect results.

Revisiting the discontinuity model I initially used following Rudd's restoration, I have treated the Morgan, Galaxy and Essential data before and after the restoration as different series. I have also centred the aggregation on the assumption that the house effects for Newspoll and Nielsen sum to zero (this may turn out to be problematic, but it is sufficient for the moment). Notwithstanding, some remaining doubts, I think this approach overcomes many of the problems my earlier discontinuity model had. I will cut to the results before reviewing the R and JAGS code.

The key finding is that Kevin was worth 5.6 percentage points in Labor's two party preferred vote share.

Turning to the house effects, we can see some of the variability in the pre-Rudd (PR) and after-Rudd (AR) values.

The revised model follows. In the first code block is the R code for managing the Morgan sample size and for separating the relevant polls into pre-Rudd (PR) and after-Rudd (AR) series. The second code block has the JAGS code. (As an aside, I have been playing with Stan lately, and might make a switch down the track).

# fudge sample size for Morgan multi - adjustment for observed over-dispersion
output.data[output.data[, 'House'] == 'Morgan multi', 'Sample'] <- 1000

# treat before and after for Morgan, Galaxy and Essential as different series
output.data$House <- paste(as.character(output.data$House), 
    ifelse(as.character(output.data$House) %in% c('Essential', 'Morgan multi', 'Galaxy'),
        ifelse(output.data[, 'Date'] >= as.Date(discontinuity), ' AR', ' PR'), ''), 
    sep='')
l <- levels(factor(output.data$House))
n <- which(l == 'Newspoll')
l[n] <- l[1]
l[1] <- 'Newspoll' # Newspoll is House number one in the factor ...
output.data$House <- factor(output.data$House, levels=l)

model {
    ## Based on Simon Jackman's original model 
    
    ## -- observational model
    for(poll in 1:NUMPOLLS) { 
        y[poll] ~ dnorm(walk[day[poll]] + houseEffect[house[poll]], samplePrecision[poll]) 
    }
            
    ## -- temporal model
    for(i in 2:PERIOD) { # for each day under analysis ...
        day2DayAdj[i] <- ifelse(i==DISCOUNTINUITYDAY, walk[i-1]+discontinuityValue, walk[i-1])
        walk[i] ~ dnorm(day2DayAdj[i], walkPrecision)
    }
    sigmaWalk ~ dunif(0, 0.01)            ## uniform prior on std. dev.  
    walkPrecision <- pow(sigmaWalk, -2)   ##   for the day-to-day random walk
    walk[1] ~ dunif(0.01, 0.99)           ## uninformative prior
    discontinuityValue ~ dunif(-0.2, 0.2) ## uninformative prior

    ## -- sum-to-zero constraint on house effects 
    for(i in 2:HOUSECOUNT) { ## vague normal priors for house effects
        houseEffect[i] ~ dnorm(0, pow(0.1, -2))
    }
    #houseEffect[NEWSPOLL] <- -sum(houseEffect[2:HOUSECOUNT])  ## all sum to zero
    houseEffect[NEWSPOLL] <- -houseEffect[NIELSEN]   ## Newspoll and Nielsen sum to zero
    #houseEffect[NEWSPOLL] <- 0 ## centred on Newspoll as zero
}

How I convert national TPP estimates into likely election outcomes

This is a short methodology discussion on how I generate a possible House of Representatives outcome from a national two-party preferred (TPP) poll estimate. This is not the most sophisticated approach possible (and it will not be the ultimate way I generate these predictions as I refine my models). But, it is what I do at the moment.

The first thing I did was to secure an estimate of the TPP vote for each seat in the 2010 federal election adjusted for boundary changes and redistributions since the 2010 election. I have shamelessly used Antony Green's pendulum for this purpose.

Secondly, building on some of the assumptions Antony made, I thought about how I should handle the treatment of independents (which sit a little outside of the mechanics of a TPP estimate). My current approach is based on the following assumptions (which are not dissimilar to Poliquant's approach):

1. Bob Katter will win Kennedy
2. Andrew Wilkie will win Denison (Labor polling in the Oz 26/06/12, ReachTEL 29/6/12)
3. Adam Bandt will lose Melbourne. [Note: this assumption is a little speculative. It rests on the Liberals changing their preference strategy from their 2010 approach (which they have said they will do). This would  see Liberal preferences flow to Labor:Greens at 67:33 in 2013 rather than the 20:80 flow in  2010. At the 2010 election Labor won 38.1 and the Greens 36.2 per cent of the primary vote. I am not aware of any subsequent polling in the Federal seat of Melbourne; but the Greens lost the 2012 by-election in the related State seat of Melbourne (where Liberals preferenced Labor ahead of Greens)].
4. Rob Oakshott will lose Lyne (Newspoll 24/10/11; ReachTEL 25/8/2011, 20/6/2012)
5. Tony Windsor will lose New England (Newspoll 24/10/11; ReachTEL 19/6/2012) 
6. Peter Slipper's seat of Fisher will be a normal Coalition/Labor contest next election 
7. Craig Thomson's seat of Dobell will be a normal Coalition/Labor contest next election
8. Tony Crook will re-contest O'Connor for the Coalition in a normal Coalition/Labor contest 

I made some tweak's to the pendulum to give effect to these assumptions.

Next I made a quick estimate of the number of seats by calculating the swing from the previous election and summing the probabilities of a win for each of the 150 seats if that swing was applied. The R-code for this function follows. As you can see, it is a short piece of code. The heavy lifting is done by the sum(pnorm(...)) functions in the middle of this code.

seatCountFromTPPbyProbabilitySum <- function(pendulumFile='./files/AntonyGreenTPP.csv', 
    pendulum, LaborTPP) {

    ALP.Outcome.2010 <- 50.12
    swing <- LaborTPP - ALP.Outcome.2010

    if(missing(pendulum)) {
        pendulum <- read.csv(pendulumFile, stringsAsFactors=FALSE)
        pendulum$ALP_TPP <- as.numeric(pendulum$ALP_TPP)
    }
    
    # Note: sd in next line comes from analysis of federal elections since 1996 ...
    ALP = round( sum( pnorm(pendulum$ALP_TPP + swing, mean=50, sd=3.27459) ) )
    
    pc <- pendulum[pendulum$OTHER == 'OTHER', ]
    OTHER = round( sum( pnorm(100 - pc$ALP_TPP - swing, mean=50, sd=3.27459) ) )

    COALITION = 150 - ALP - OTHER # Just to ensure it all adds to 150.

    # return a data frame - makes it easier to ggplot later
    results <-                data.frame(Party='Other',         Seats=OTHER)
    results <- rbind(results, data.frame(Party='Coalition',     Seats=COALITION))
    results <- rbind(results, data.frame(Party=factor('Labor'), Seats=ALP))
    
    return(results)
}

Update: I have updated the model to better manage how I treat Denison.

seatCountFromTPPbyProbabilitySum <- function(pendulumFile='./files/AntonyGreenTPP.csv', 
    pendulum, LaborTPP) {

    ALP.Outcome.2010 <- 50.12
    swing <- LaborTPP - ALP.Outcome.2010

    if(missing(pendulum)) {
        pendulum <- read.csv(pendulumFile, stringsAsFactors=FALSE)
        pendulum$ALP_TPP <- as.numeric(pendulum$ALP_TPP)
    }
    
    # Note: sd in next few lines comes from analysis of federal elections since 1996 ...
    pc <- pendulum[pendulum$OTHER == 'OTHER', ]
    other.raw <- sum( pnorm(100 - pc$ALP_TPP - swing, mean=50, sd=3.27459) )
    OTHER <- round( other.raw )
    carry <- other.raw - OTHER

    # this approach typically favours Labor (probably the right way to go)
    ALP <- round( carry + sum( pnorm(pendulum$ALP_TPP + swing, mean=50, sd=3.27459) ) )
    
    COALITION <- 150 - ALP - OTHER # Just to ensure it all adds to 150.

    # return a data frame - makes it easier to ggplot later
    results <-                data.frame(Party='Other',         Seats=OTHER)
    results <- rbind(results, data.frame(Party='Coalition',     Seats=COALITION))
    results <- rbind(results, data.frame(Party=factor('Labor'), Seats=ALP))
    
    return(results)
}

From this function we can plot a likely election outcome for a given a swing.

To get a more nuanced understanding of a potential election outcome, I undertake a simple Monte Carlo simulation (typically with 100,000 iterations). This is not a Bayesian MCMC approach. It's just a plain old fashioned MC simulation. The R-code for this procedure is more substantial.

storeResult <- function(N, pendulum, individualSeats=FALSE) {
    # Use of R's lexical scoping

    # entry sanity checks ...
    stopifnot(is.numeric(N))
    stopifnot(is.data.frame(pendulum))
    stopifnot(N > 0)
    seatCount <- nrow(pendulum)
    stopifnot(seatCount > 0)

    # sanity checking variables
    count <- 0
    finalised <- FALSE

    # where I store the house wins ...
    ALP <- rep(0, length=seatCount)
    COALITION <- rep(0, length=seatCount)
    OTHER <- rep(0, length=seatCount)
    
    CUM_ALP <- rep(0, length=seatCount)
    CUM_COALITION <- rep(0, length=seatCount)
    
    # where I keep the seat-by-seat wins
    seats <- data.frame(seat=pendulum$SEAT, state=pendulum$STATE, Labor=ALP, 
        Coalition=COALITION, Other=OTHER)

    rememberSim <- function(simResult) {

        # - sanity checker
        stopifnot(!finalised)
        stopifnot(count < N)
        count <<- count + 1
        stopifnot(length(simResult) == seatCount)
    
        # - overall result
        a <- table(simResult)
        ALP[ a[names(a)=='ALP'] ] <<- ALP[ a[names(a)=='ALP'] ] + 1
        COALITION[ a[names(a)=='COALITION'] ] <<- 
            COALITION[ a[names(a)=='COALITION'] ] + 1
        OTHER[ a[names(a)=='OTHER'] ] <<- OTHER[ a[names(a)=='OTHER'] ] + 1
        
        # - seat by seat result
        if(individualSeats) {
            seats$Labor <<- ifelse(simResult == 'ALP', seats$Labor + 1,
                seats$Labor)
            seats$Coalition <<- ifelse(simResult == 'COALITION', 
                seats$Coalition + 1, seats$Coalition)
            seats$Other <<- ifelse(simResult == 'OTHER', seats$Other + 1, 
                seats$Other)
        }
    }
    
    finalise <- function() {
        # sanity checker
        stopifnot(!finalised)
        stopifnot(count == N)
        
        ALP <<- ALP / N
        COALITION <<- COALITION / N
        OTHER <<- OTHER / N
        
        if(individualSeats) {
            seats$Labor <<- seats$Labor / N
            seats$Coalition <<- seats$Coalition / N
            seats$Other <<- seats$Other / N
        }
        
        for(i in 1:seatCount) {
            CUM_ALP[i] <<- 1 - sum(ALP[1:i])
            CUM_COALITION[i] <<- 1 - sum(COALITION[1:i])
        }
                
        finalised <<- TRUE
    }        
    
    results <- function() {
        stopifnot(finalised)
        data.frame(seatsWon=1:nrow(pendulum), Labor=ALP, Coalition=COALITION, 
            Other=OTHER)
    }
        
    cumResults <- function() {
        stopifnot(finalised)
        data.frame(seatsWon=1:nrow(pendulum), Labor=CUM_ALP, Coalition=CUM_COALITION)
    }

    winProbabilities <- function() {
        stopifnot(finalised)
        win <- (floor(seatCount/2) + 1):seatCount
        list(Labor = sum(ALP[win]), Coalition = sum(COALITION[win]))
    }

    seatResults <- function() {
        stopifnot(finalised)
        stopifnot(individualSeats)
        seats
    }
    
    list(rememberSim=rememberSim, finalise=finalise, results=results, cumResults=cumResults,
        seatResults=seatResults, winProbabilities=winProbabilities)
}

## -- similate one Federal election
simulateNationaLResult <- function(pendulum, swing) {
    rawPrediction <- pendulum$ALP_TPP + swing
    probabilisticPrediction <- rawPrediction + rnorm(nrow(pendulum), mean=0, sd=3.27459)
    ifelse(probabilisticPrediction >= 50, 'ALP', pendulum$OTHER)
}

## -- run N simulations of one Federal election outcome
simulateOneOutcome <- function(N=100000, pendulumFile='./files/AntonyGreenTPP.csv',
    pendulum, LaborTPP, individualSeats=FALSE) {

    ALP.Outcome.2010 <- 50.12
    swing <- LaborTPP - ALP.Outcome.2010

    if(missing(pendulum)) {
        pendulum <- read.csv(pendulumFile, stringsAsFactors=FALSE)
        pendulum$ALP_TPP <- as.numeric(pendulum$ALP_TPP)
    }
    
    r <- storeResult(N, pendulum, individualSeats)
    for(i in 1:N) r$rememberSim ( simulateNationaLResult(pendulum, swing) )
    r$finalise()

    invisible(r)
}

From this simulation, there are a few plots I can make:

I am currently working on a state-level frame for converting a series of state TPP estimates to a national outcome for the House of Representatives.

Cube Law

I have spent the past few days playing with Bayesian statistics, courtesy of JAGS (which is a Markov chain Monte Carlo (MCMC) engine where the acronym stands for Just Another Gibbs Sampler).

The problem I have been wrestling with is what the British call the Cube Law. In first past the post voting systems, with a two-party outcome, the Cube Law asserts that the ratio of seats a party wins at an election is approximately a cube of the ratio of votes the party won in that election. We can express this algebraically as follows (where s is the proportion of seats won by a party and v is the proportion of votes won by the party). Both s and v lie in the range from 0 to 1.

My question was whether the relationship held up under Australia's two-party-preferred voting system. For the record, I came across this formula in Simon Jackman's rather challenging text: Bayesian Analysis for the Social Sciences.

My first challenge was to make the formula tractable for analysis. I could not work out how Jackman did his analysis (in part because I could not work out how to generate an inverse gamma distribution from within JAGS, and it did not dawn on me initially to just use normal distributions). So I decided to pick at the edges of the problem and see if there was another way to get to grips with it. There are a few ways of algebraically rearranging the Cube Law identity. In the first of the following equations, I have made the power term (relabeled k) the subject of the equation, In the second, I made the proportion of seats won the subject of the equation.

In the end, I decided to run with the second equation, largely because I thought it could be modeled simply from the beta distribution which provides diverse output in the range 0 to 1. The next challenge was to construct a linking function from the second equation to the beta distribution. I am not sure whether my JAGS solution is efficient or correct, but here goes (constructive criticism welcomed).

    model {
        # likelihood function
        for(i in 1:length(s)) {
            s[i] ~ dbeta(alpha[i], beta[i]) # s is a proportion between 0 and 1
            alpha[i] <- theta[i] * phi
            beta[i]  <- (1-theta[i]) * phi
            theta[i] <- v[i]^k / ( v[i]^k + (1 - v[i])^k )   # Cube Law
        }
        # prior distributions
        phi ~ dgamma(0.01, 0.01)
        k ~ dnorm(0, 1 / (sigma ^ 2))   # vaguely informative prior
        sigma ~ dnorm(0, 1/10000) I(0,) # uninformative prior, positive
    }

The results were interesting. I used the Wikipedia data for Federal elections since 1937. And I framed the analysis from the ALP perspective (ALP TPP vote share and the ALP proportion of seats won).

The mean result for k was 2.94. The posterior distribution for k had a 95% credibility interval between 2.282 and 3.606. The median in the posterior distribution was 2.939 (pretty well the same as the mean; and both were very close to the magical 3 of the Cube Law). It would appear that the Federal Parliament, in terms of the ALP share of TPP vote and seats won operates pretty close to the Cube Law. The distribution of k, over 4 chains each with 50,000 iterations of the MCMC was:

The files I used in this analysis can be found here.

Technical follow-up: Simon Jackman deals with the Cube Law with what looks like an equation from a classical linear regression of logits (logs of odds). The core of this regression equation is as follows:

By way of comparison, the k in my equation is algebraically analogous to the β1 in Jackman's equation. Our results are close: I found a mean of 2.94, Jackman a mean of 3.04. In my equation, I implicitly treat β0 as zero. Jackman found a mean of -0.11. He uses the β0 to asses bias in the electoral system. Nonetheless, the density kernel I found for k (below) looks very similar to the kernel Jackman found for his β1 on page 149 of his text. [This last result may surprise a little as my data spanned the period 1937 to 2010, while Jackman's data spanned a shorter period: 1949 to 2004].

I suspect the pedagogic point of this example in Jackman's text was the demonstration of a particular "improper" prior density  and the use of its conjugate posterior density. I suspect I could have used Jackman's approach with normal priors and posteriors. For me it was a useful learning experience looking at other approaches as a result of not knowing how get an inverse gamma distribution working in JAGS. Nonetheless, if you know how to do the inverse gamma, please let me know. 

2013 Federal Election Simulation

Following his success at predicting the US Presidential election for 2012, I decided to read Nate Silver's book: The Signal and the Noise. I am enjoying it so far.

It inspired me to throw together a quick and dirty Monte Carlo simulation of the Australian 2013 Federal Election based on the current state of national opinion polling (assuming the Coalition would win 52.5 per cent of the two party preferred vote). In effect, the simulation runs the election 100,000 times collating the range of potential outcomes given the current polling.

The result is a probability distribution of likely seats won for each party. The most likely outcome is 63 seats for the ALP and 85 seats for the Coalition. It has a more than 50 per cent chance of being between 62 and 65 seats for the ALP, and between 83 and 86 seats for the Coalition.

I am assuming that Wilkie would retain Denison and that Katter would retain Kennedy. For the other seats currently held by the greens and independents, I am assuming they would return to the major parties. I am following Poliquant's logic on this. I used Antony Green's pendulum for some of the underlying arithmetic.

This is a very naive model! Over coming months I plan to refine it. Refinements would include better treatment of the non-major-party seats, adding state-by-state opinion polling results, accounting for systemic bias from the various pollsters (what Simon Jackman calls "house effects"), and a sprinkling of Bayesian intelligence. As we get closer to the election I will look at ballot position and whether the seat contest has a retiring member or not.

In addition to the above national totals, the model could also provide seat-by-seat probabilities. I have not written code for this as yet, but it should not take too long.

Update
Well, I have now written some code to aggregate the individual seat results from each simulation run. I have made some minor modifications to the way in which the second model manages the seats held by others. I have also tidied the plots a little. Nonetheless, it is still a very naive model. I have a long way to go before I am comfortable that it is making robust predictions on the available data.

The updated charts follow for the new 100,000 elections simulation:

The R code for these models can be found here. As always, if you see any errors in the code, or ways I can improve the analysis, please drop me a line.