Supreme Court Puzzle

This is a simulation in R for the solution to Fivethirtyeight’s Riddler Classic from April 14th 2017, found here.

Imagine that U.S. Supreme Court nominees are only confirmed if the same party holds the presidency and the Senate. What is the expected number of vacancies on the bench in the long run?

You can assume the following:

  • You start with an empty, nine-person bench.
  • There are two parties, and each has a 50 percent chance of winning the presidency and a 50 percent chance of winning the Senate in each election.
  • The outcomes of Senate elections and presidential elections are independent.
  • The length of time for which a justice serves is uniformly distributed between zero and 40 years.

The following code usually gives an answer between 0.68 and 0.74. Running it for longer gives something around 0.72.

library(tidyverse)
justices <- seq(from = 0, to = 0, length.out = 9)
parties <- c('Democrat', 'Republican')

pres <- sample(parties, 1)
senate <- sample(parties, 1)

fillSeats <- function(justices, year, senate, pres){
  if (year == 0){
    return(justices)
  }

  else{
    if (pres == senate){
      for (num in 1:length(justices)){
        if (justices[num] == 0){
          justices[num] <- sample(0:40, 1)      
        }
        
      }
    }
    justices
  }
  
}


sim <- function(time){
  justices <- seq(from = 0, to = 0, length.out = 9)
  histVacancies <- c()
  histMeanVacancies <- c()
  period <- 0:time
  
  for (year in period){
    justices <- justices - 1
    justices <- replace(justices, which(justices == -1), 0)
    
    justices <- fillSeats(justices, year, senate, pres)
    vacancies <- sum(justices == 0)
    
    if (year %% 2 == 0){
      senate <- sample(parties, 1)    
    }
    
    if (year %% 4 == 0){
      pres <- sample(parties, 1)    
    }
#    print(c(senate, pres))

    histVacancies <- c(histVacancies, vacancies)
    meanVacancies <- mean(histVacancies)
    histMeanVacancies <- c(histMeanVacancies, meanVacancies)
    
  }
  histMeanVacancies
}

vec <- sim(10000)
print(tail(vec, 1))
## [1] 0.7422258

And here’s a graph. I cut off the limits on the y-axis for ease of view.

df <- data.frame(year = 0:(length(vec)-1), meanVacancies = vec)

ggplot(df, aes(x = year, y = meanVacancies)) +
  geom_line() +
  scale_y_continuous(limits = c(0.5,1))