Using Simulation to Model the Optimal Martingale Betting Strategy in Roulette

Introduction

Roulette is a game characterized by pure chance and probability, and I came across it shortly after taking a course in engineering probability. This prompted me to look into the game and its different strategies, and to try and make a model that will get me the best results.

Challenge

The problem I wanted to solve is this: Given an initial stack, what is the maximum/optimal initial Martingale bet size? 

There are many different subsets of problems that exist within this one, but we will narrow down our focus to a specific problem:
Given an initial stack, what is the maximum initial Maringale bet size with a 90% certainty of positive return?

Context: In solving this problem, bets will only be made on outside bets (betting on a color) and I will be simulating under American Roulette conditions. Outside bets have a slightly less than 50% chance of winning and the payout for an outside bet win is equal to the size of the bet. In this way, we can model each game of roulette as a Bernoulli random variable with a discrete probability distribution. I will be using Python to run simulations.

Approach

1. Data Generation

Here are the steps I took in solving this problem with simulation:

1. Create a function that simulates one game of roulette and returns the new stack size. Function header and input parameters are as follows:

def roulette(stack_size, bet_size, win_probability):

...

return updated_stack_size

2. Create a function that simulates a chain of roulette games under the Martingale strategy. The bet size will double if a loss occurs. The function will return a list of values representing the stack size after each game of roulette. The simulation stops when one win occurs and the ending stack size is greater than the initial stack size, or if the current bet size becomes greater than the remaining stack size. Function header and input parameters are as follows:

def martingale(stack_size, bet_size, win_probability):

...

return [stack_sizes]

3. Create a function that will simulate many Martingale chains at different bet values. In my case, I want to simulate Martingale chains with initial stack = 100, bet values from 0 to 100 at intervals of 0.1, and 20,000 simulated chains at each bet value. At each bet value, I want to determine the conversion rate (percentage of chains that end in positive return). The function will return two lists: a list with the bet values simulated on and a list with corresponding conversion rates. Function header and input parameters are as follows:

def conversion_rates_sim(stack_size = 100, bet_values = np.linspace(0,100,1001), num_sims = 20000):

...

return [bet_values], [conversion_rates]

4. Plot the generated simulation data.

2. Data Analysis

Looking at the data we generated, we can see a unique behavior. The data splits into two trends from bet sizes between 0 and 30, and there is a piecewise behavior after bet sizes of 30. We are going to look into bet sizes between 0 and 30, as that range of bet sizes covers the conversion rates between 0.7 and 1. Specifically, I am simulating 1329 bet sizes linearly spaced between bet size 0 and 33.2, with 20,000 simulations at each bet size. To generate data for this range, we are going to call:

conversion_rates_sim(stack_size = 100, bet_values = np.linspace(0,33.2,1329), num_sims = 20000)

This is the graph of the resulting data:

Here, we can see a close up of the graph behavior in our first range. The data seems to fit two trend lines, but I will be taking the average trend to model the true value. I will be fitting a second-degree polynomial regression line over the data. The results of the model are below:

With these results, the model I will use to approximate the function of bet size vs. conversion rate is the following:

let
x = bet size
y = conversion rate

y = -0.014156x² + 0.00016576x + 0.9991425

Work in progress. To be continued!