What is the probability of having all wealth concentrated into one single pocket after n trades? On Joe Rogan's podcast #1769, Dr. Jordan Peterson said something like this:

*

*Take a population of 10 people, each starting with 100 dollars.

*They will "trade" based on a coin flip. Whoever "loses" has to give the other person one dollar.

*In the long run, all wealth will be concentrated into one pocket.

At first, I couldn't believe this was true because of the equilibrium randomness would bring to this system (compared to an "informed trader" scenario). Each player would win as much as they would lose, resulting in a zero-sum endless game.
A player's expected wealth after $ n $ trades:
$$ E[w] = w_0 + \sum_{i=0}^n {1 \over 2} \times 1 + \sum_{i=0}^n {1 \over 2} \times -1$$
$$ E[w] = w_0 + {n \over 2} - {n \over 2} $$
$$ E[w] = w_0 $$
I ran a Python simulation, and in the very first run I got this:
It took 70431 individual trades for all wealth to be concentrated in one pocket.
Am I missing something here?

This is the code for the simulation:
https://gist.github.com/victorvalentee/1c04b2d16005f7372273b153fe9ece23

The answer by Aksakal already explains it enough why the wealth is ending up into one single pocket. In this answer we elaborate a bit further on the quantitative approach by Matt F.
One could compute the probability in terms of the probability of a specific player remaining in the game and consider the position of the players independent.
Then there are two ways

*

*The probability that only one player is in the game is equal to the probability that a player has reached the maximum amount of money. Let's call the probability for an independent player to win $p_{win}$.
The probability that none of the players have won is $$(1-p_{win})^n$$


*The probability that only one player is in the game is equal to the probability that no more than $n-2$ players have lost their money yet. Let's call this probability for an independent player to loose $p_{loose}$.
The probability that none of the players have won, is the probability that less than 9 have lost and is $$1-(p_{loose})^{n}-n(1-p_{loose})(p_{loose})^{n-1}$$
However, when I apply these formulas then I get small discrepancies. One way overestimates the number of trials before the game ends, and the other underestimates it.
The reason for the discrepancy is that the scores of the individual players are correlated. If one player has lost after $n$ games, then this reduces the probability for the other players to have lost after $n$ games.

In the image, we also plotted an inverse Gaussian distribution fitted with the method of moments, and it seems to agree reasonably well with the shape of the distribution determined from the simulations. So possibly an inverse gaussian could be used as a model for the distribution (but one would still need to figure out the mean and variance).
### function to simulate game
sim = function(k=100,n=5) {
  # x: is a variable that keeps track of values, 
  # the size of this variable will reduce 
  # once participants hit zero
  x = rep(k,n)    
  ### n_active: number of active participants
  ### these equals the lenght of 'x' but we do not 
  ### want to compute that length everytime
  n_active = n
  ### counts: keeping track of the number of games
  counts = 0
  while (n_active>1) {
    counts = counts + 1
    
    ### sample two players
    s = sample(1:n_active,2)
    ### add +1 and -1 to the wealth of the players
    x[s] = x[s] + c(1,-1)
    
    ### check for zero wealth and remove the participant from 'x' if neccesary
    if (x[s[1]]==0) {
      x = x[-s[1]]
      n_active = n_active -1 
    }
    if (x[s[2]]==0) {
      x = x[-s[2]]
      n_active = n_active -1 
    }
  }
  return(counts)
}

### simulations
set.seed(1)
y = replicate(10^3,sim())

### plot histogram
hist(y, breaks = seq(0,10^6,10^4), main = "histogram of simulations \n compared with approximations")


### add estimates based on an individual's probability to loose or win after n turns
n = seq(0,10^6,10^4)

p1 = pnorm(500,mean = 100,sd = sqrt(0.4*n))-pnorm(-500,mean = 100,sd = sqrt(0.4*n))
lines(n[-1],-diff(p1^5)*10^3, col = 2)

p2 = pnorm(0,mean = 100,sd = sqrt(0.4*n))*2
lines(n[-1],-diff(pbinom(1,5,p2))*10^3, col = 3)

### add inverse gaussian curve fitted with method of moments
lines(n,statmod::dinvgauss(n, mean(y), dispersion = var(y)/mean(y)^3)*10^7, col = 4)

legend(4*10^5, 80,
       c("estimate based on probability of winning", "estimate based on probability of loosing",
         "inverse Gauss distribution"),
       cex = 0.7, col = c(2,3,4), lty = 1)

