The requirements for this Monte Carlo simulation are taken from Giordano, Fox, & Horton, “A First Course in Mathematical Modeling”, Fifth Edition, Cengage Learning, 2013: Chapter 5.3, page 199, Project 1.

NOTE: The R code used for Blackjack simulation can be found in the attached Appendix.

Project 1 in Chapter 5.3 asks us to construct a Monte Carlo simulation of the card game called Blackjack. We are told that in addition to the basic rules of Blackjack, the simulation must adhere to the following requirements:

Bets of $2 will be placed on each hand by the “player”, with no limit on the number of $2 bets that can be placed.
A win of a hand will earn the player $2 above and beyond their initial $2 bet. A loss of a hand will result in the player losing the entirety of their $2 bet. If a hand results in a tie between the player and the dealer, no money is either won or lost.
The game will start using two decks of cards. Hands will be played until all cards from the two decks have been used.
When the two decks of cards are exhausted, the player’s running total winnings/losses for all hands played with the two decks of cards are saved.
When two decks of cards are exhausted, a new set of two decks is initiated and play resumes, with the player’s winnings/losses tally being accummulated anew.
A total of 12 sets of two decks of cards should be “played”.
When all 12 set of two decks have been played, we must average the winning/losses from all 12 iterations to determine the overall performance of the player.
The player cannot see any of the dealer’s cards

The simulation must therefore be capable of generating multiple two-deck sets of cards for purposes of enabling the required 12 iterations of two-deck play. Furthermore, the simulation must also carefully manage the dealing of cards to both the player and the dealer, with a variety of constraints being checked as each card is dealt.

For example, since the required simulation requires the play of nearly 12 * 2 * 52 = 1248 cards, it is simply not practical to require a player to interactively play the blackjack game. In fact, the requirements for the simulation state quite clearly that we must choose a predefined strategy to play for the player and then play that strategy throughout the entire simulation. As such, we assume that the player will always accept another card if their current hand totals 15 or less, and will always stand whenever their hand totals 16 or more.

Similarly, per the instructions the dealer will always accept another card if their hand totals 16 or less and will always stand whenever their hand totals 17 or more.

The shuffled 2-deck set of cards is then “played” until all cards have been used. A total of 12 iterations of ‘play’ are simulated, with the player’s winnings from each simulation stored in a numeric vector. Finally, the vector of winnings is averaged to determine the player’s average winnings per hand, with the results displayed onscreen.

iteration	res.winnings
1	-14
2	16
3	-8
4	2
5	6
6	-2
7	14
8	16
9	8
10	0
11	-8
12	2

## Average Winnings per Hand =  2.66666666666667

The results shown above represent a fairly limited sample of possible player winnings relative to the betting strategy described in the requirements for this project. To get a better sense of what a player’s likely winnings would be we can run a large number of iterations (e.g., 1000) of the above simulation and then calculate the average winnings for all iterations.

Summary statistics for the player’s average winnings from 1000 iterations of playing 12 sets of 2 decks of Blackjack with the player standing if the value of their hand is 16 or greater are shown below:

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -12.670  -4.333  -2.833  -2.816  -1.167   4.667

Based on the Blackjack rules we were required to adhere to, the player’s average losses when playing 12 sets of 2 decks of cards is approximately ($2.81). As such, a player should expect to lose money if they adhere to the “stand at 16 or greater” strategy.

A histogram of the average winnings of 1000 iterations of playing 12 sets of 2 decks of Blackjack shows that the distribution of the average winnings is approximately normal, and also reinforces the conclusion that the player should on average expect to lose money when adhering to the “stand at 16 or greater” strategy.

We can attempt to pinpoint the player strategy that minimizes the expected average losses by evaluating “stand” values ranging from 11 to 20 via further simulation iterations:

stand	avg_winnings
11	-5.148667
12	-3.723333
13	-3.248833
14	-2.930833
15	-2.589167
16	-2.683167
17	-3.011333
18	-4.290167
19	-7.613667
20	-12.870833

The table above shows that the average player loss gradually declines as they increase their Blackjack “stand at” value from 11 to 15. However, increasing the “stand at” hurdle to any value above 15 will result in average losses greater than would be expected if the player adhered to a “stand at” value of 15. As such, if the player wishes to minimize their average losses, they should adhere to a Blackjack strategy of “stand at 15”.

————————————————————————————————————————–

Appendix

Complete R Code for Blackjack Simulation

We start by defining a collection of functions that will be used to facilitate the “play” of both the dealer and the player.

#########################################################################
# Function sum_cards: sums the values of the cards in a blackjack hand. Aces
# are counted as either 11 or 1 depending on the sum total of non-ace cards in the hand
#
# Parameters:
# - hand: a blackjack hand that has been dealt to either a player or dealer. 
#
# NOTE: Blackjack hands in this simulation consist solely of integer values since the
#       cards do not need to be displayed graphically. Suits of cards also are irrelevant
#       for similar reasons
# ========================================================================

sum_cards <- function(hand) {
  
  aces <- FALSE
  ace_locs <- numeric()
  sum <- 0
  
  # if the hand is only 2 cards, just sum and exit
  if (length(hand) == 2) {
    sum <- hand[1] + hand[2]
    return(sum)
  }
  
  # find aces
  for (k in 1:length(hand)) {
    if(hand[k] == 11) {
      aces <- TRUE
      ace_locs <- c(ace_locs, k)
    }
  }
  
  # if no aces were found, sum all items and return
  if (aces == FALSE) {
    for (i in hand) {
       sum <- sum + i
    }
    return(sum)
  }
  
  # otherwise aces were found so first sum items that are NOT aces
  not_aces <- subset(hand, !(hand %in% c(11)))
  for (i in not_aces) {
    sum <- sum + i
  }
  
  # if sum of non-aces > 10 all aces must be treated as 1's
  if (sum > 10) {
    for (i in ace_locs) {
      sum <- sum + 1
    }
    return(sum)
  }
  
  # if sum of non-aces <= 9, need to check count of aces. 
  # If > 1 one can be treated as 11 while other must be treated as 1
  if (sum <= 9 & length(ace_locs) > 1) {
    sum <- sum + 11
    for (i in 1:(length(ace_locs) - 1)) {
      sum <- sum + 1
    }
    return(sum)
  }
  
  # if sum of non-aces == 10, need to check count of aces.
  # if == 1, then add 11 to sum and exit. if > 1, treat all aces as 1's
  if (sum == 10 & length(ace_locs) == 1) {
    sum <- sum + 11
    return(sum)
  } else {
    for (i in 1:length(ace_locs)) {
      sum <- sum + 1
    }
    return(sum)
  }
  
}

#########################################################################
# Function player_h: implements standard blackjack rules logic for the player's hand
# Parameters:
# - p_hand: the initial 2 cards dealt to the player
# - cards: a vector containing a pre-shuffled set of cards
# - i: an index variable that indicates the position of the card most recently dealt
#
# NOTE: Player strategy will be to take another card any time hand is 15 or less. This 
#       implies that the player will stand with any hand of 16 or greater.
# ========================================================================

player_h <- function(p_hand, cards, i, p_stand){

  ncards <- length(cards)
  stand <- FALSE
  
  while (stand != TRUE) {
    # sum current hand
    current <- sum_cards(p_hand)
  
    # if card index == number of cards in decks, time to stop
    if (i == ncards) {
      stand <- TRUE
      
    # else if hand is < p_stand value, take another card  
    } else if (current < p_stand) {
      i <- i + 1
      p_hand <- c(p_hand, cards[i])
    } else {
      stand <- TRUE
    }
    
  } # end while stand != TRUE
  
  return(c(current, i))
}

#########################################################################
# Function dealer_h: implements standard blackjack rules logic for the dealer's hand
# Parameters:
# - d_hand: the initial 2 cards dealt to the dealer
# - cards: a vector containing a pre-shuffled set of cards
# - i: an index variable that indicates the position of the card most recently dealt
#
# NOTE: Dealer MUST stand on hand total of 17 or greater, so any time dealer hand is <= 16
#       another card must be drawn
# ========================================================================

dealer_h <- function(d_hand, cards, i){
  
  ncards <- length(cards)
  stand <- FALSE
  
  while (stand != TRUE) {
    # sum current hand
    current <- sum_cards(d_hand)
    
    # if card index == number of cards in decks, time to stop
    if (i == ncards) {
      stand <- TRUE
      
    # else if hand is <= 16, take another card
    } else if (current <= 16) {
      i <- i + 1
      d_hand <- c(d_hand, cards[i])
    } else {
      stand <- TRUE
    }
    
  } # end while stand != TRUE

  return(c(current, i))
}

#########################################################################
# Function play_bj: plays blackjack until 'cards' are exhausted and tallies player's
# winnings and losses
#
# Parameters:
# - cards: a pre-shuffled set of cards
# - bet: the dollar amount of bets placeable by the player for each blackjack hand
#
# NOTE: need to sample WITHOUT replacement from cards. Since they are 
# pre-shuffled before calling this function, just deal from top of deck and 
# use and index to keep track of which cards haven't been dealt
# ========================================================================

play_bj <- function(cards, bet, p_stand){
  
  # initialize the number of cards to be played
  ncards <- length(cards)
  
  # initialize index for cards
  i <- 0
  
  # initialize winnings accummulator
  winnings <- 0
  
  # now loop until all cards used
  while (i < ncards) {
    
    # New hand needed: check to see whether at least 4 cards remain unused
    if (i <= (ncards - 4)) {
      
      #### Deal 2 cards to player
      # increment card index
      i <- i + 2
      p_hand <- c(cards[i-1], cards[i])
      
      ### Deal 2 cards to dealer
      # increment card index
      i <- i + 2
      d_hand <- c(cards[i-1], cards[i])
      
    } else {
      # else cards are exhausted so exit
      return(winnings)
    }
    
    #### logic for player's hand
    p_res <- player_h(p_hand, cards, i, p_stand)
    
    # update card index
    i <- p_res[2]
    
    # if player didn't go over 21 then apply dealer hand logic
    if (p_res[1] <= 21) {
      d_res <- dealer_h(d_hand, cards, i)
      
      # update card index
      i <- d_res[2]
    }
    
    # ---------------------------------------------------------------
    # now compare player's hand to dealer's hand to find out who won
    
    # if player hand exceeded 21 subtract 2 from winnings
    if (p_res[1] > 21) {
      winnings <- winnings - 2
      
    # else if dealer went over 21 dealer has lost so add 2 to winnings  
    } else if (d_res[1] > 21) {
      winnings <- winnings + 2
      
    # else if player hand > dealer hand add 2 to winnings
    } else if (p_res[1] > d_res[1]) {
      winnings <- winnings + 2
      
    # else if player hand < dealer hand subtract 2 from winnings
    } else if (p_res[1] < d_res[1]) {
      winnings <- winnings - 2
    }
      
    
  } # end while (i < ncards)
  return(winnings)
}

The ‘main’ portion of the simulation generates and shuffles a 2-deck set of cards, with each face card represented by its numeric Blackjack value. For example, jacks, queens, and kings are all treated as ‘10’ while aces are initially treated as ‘11’. (NOTE: Aces can also be treated as having a value of ‘1’ depending upon the values of other cards present in any given Blackjack hand. The programming logic for determining the value of an ace in any given hand is contained within the sum_Cards function defined above).

#########################################################################
# This is the 'main' portion of the simulation

# NOTE: Suits of cards don't really matter at all. All we really need to track is the number remaining for each class of card, e.g., how many 2's, how many 3's, etc..

# set number of 2-deck sets to play
N <- 12

# set size of bets
bet <- 2

# set hand total for player to stand at
p_stand <- 16

# generate a deck of cards
deck <- rep(c(2:10, 10, 10, 10, 11), 4)

# create a 2 deck set of cards
two_d <- rep(deck,2)

# create a vector to store results of play
res <- data.frame(N = 1:12, winnings = numeric(N))

for (i in 1:N){
  # shuffle 2 new decks of cards
  s_decks <-sample(two_d,length(two_d), replace = FALSE)
  s_decks

  # play until 2 decks are exhausted
  res$winnings[i] <- play_bj(s_decks, bet, p_stand)
}

# create a table of winnings for each iteration
iteration <- 1:N
kable(data.frame(iteration, res$winnings))

# average the winnings and print to screen
avg_winnings <- mean(res$winnings)
message(paste("Average Winnings per Hand = ", avg_winnings))

R code for 1000 iterations of simulation

# set number of simulation iterations to be run
iters <- 1000

# create a vector to store results of play
res2 <- data.frame(N = 1:iters, avg_winnings = numeric(iters))

for (k in 1:iters) {
  for (i in 1:N){
    # shuffle 2 new decks of cards
    s_decks <-sample(two_d,length(two_d), replace = FALSE)
    s_decks

    # play until 2 decks are exhausted
    res$winnings[i] <- play_bj(s_decks, bet, p_stand)
  } # for i

  # save the average the winnings to the results vector
  res2$avg_winnings[k] <- mean(res$winnings)
  
} # for k

Evaluating “stand” values ranging from 11 to 20 via further simulation iterations:

set.seed(123)

# create a vector to store the average winnings for various player "stand" values
playw <- data.frame(stand = 11:20, avg_winnings = numeric(10))

p_stand <- 11

for (z in 1:10) {
  res4 <- data.frame(N = 1:iters, avg_winnings = numeric(iters))

  for (k in 1:iters) {
    for (i in 1:N){
      # shuffle 2 new decks of cards
      s_decks <-sample(two_d,length(two_d), replace = FALSE)
      s_decks

      # play until 2 decks are exhausted
      res$winnings[i] <- play_bj(s_decks, bet, p_stand)
    } # for i

  # average the winnings and print to screen
  res4$avg_winnings[k] <- mean(res$winnings)
  
  } # for k
  
  # save the average winnings to the results vector
  playw$avg_winnings[z] <- mean(res4$avg_winnings)
  
  # increase the player's stand value
  p_stand <- p_stand + 1
  
} # for z

kable(playw)

Blackjack Simulation via Monte Carlo Methods

James Topor

December 16, 2016

Appendix

Complete R Code for Blackjack Simulation