Monte Carlo Simulation

Probability

Monte Carlo Method

The Monte Carlo method is a general term for methods that use random numbers to estimate the expectation and variance of random variables. Monte Carlo is apparently the name of a district in the Principality of Monaco (where casinos are popular, with the Monte Carlo Casino being the first casino built). Since the expectation and variance:

\[ \small \begin{align*} &E[X] = \int_D xp(x)dx \\ &V[X] = \int_D (x-E[x])^2p(x)dx \end{align*} \]

are calculated by integrating over the domain \(\small D\) of the probability distribution, Monte Carlo methods are often considered a form of numerical integration. In particular, integrals with many dimensions are difficult to evaluate using analytical formulas or numerical integration using discretization, so the Monte Carlo method is often used.

 For example, if a random variable \(\small x\sim N(\mu,\sigma^2)\) that follows a normal distribution is given, when you want to find the expected value of a function \(\small f(x)\) whose value depends on it, you can approximate

\[ \small E\left[f(x)\right] \approx \frac{1}{n} \sum_{i=1}^n f(x_i) \]

by expressing the sample of the random variable \(\small x\) as \(\small x_1,\cdots,x_n\). It is easy to guess that as the number of samples is increased, the probability distribution of the samples will converge to the original probability distribution, and so converge to

\[ \small \lim_{n\rightarrow \infty} \frac{1}{n} \sum_{i=1}^n f(x_i) \rightarrow E\left[f(x)\right]. \]

In probability theory and statistics, the degree of certainty of this convergence is discussed in detail, but we will not go into that here.

 The Monte Carlo method requires generating random numbers that follow a specified probability distribution \(\small F(x)\). There are various methodologies for this depending on the probability distribution, but a universal method is to generate \(\small x_i=F^{-1}(u_i)\) using uniformly distributed random numbers \(\small u_i\sim U[0,1]\) (functions for generating this are implemented in all programming languages). This is called the inverse function method. It is generally difficult to generate multivariate random numbers, so they are often generated using a method based on Bayesian statistics known as the Markov chain Monte Carlo method. Multivariate normal random numbers can be generated relatively easily, and the method for generating them will be described later.

Random Walk and Brownian Motion

 The Monte Carlo method is also often used to simulate stochastic processes, and in financial engineering it is used to estimate the present value of complex financial products by simulating stock prices, exchange rates and etc., using random numbers. A stochastic process is a general term for random variables \(\small x(t)\) that are generated as a function of time. For example, if \(\small w_1,\cdots,w_n\) follow the standard normal distribution \(\small N(0,1)\), the stochastic process represented as:

\[ \small \begin{align*} &w(t) = w_t \\ &r(t) = \sum_{i=1}^t w_i, \quad t=1,2,\cdots, n \end{align*} \]

is called a white noise process and a random walk, respectively.

 The concept of random walk as a continuous-time stochastic process is called Brownian motion. Brownian motion is a stochastic process whose variance follows a normal distribution proportional to time \(\small t\), and

\[ \small B(t)\sim N(\mu t, \sigma^2t) \]

is a probability distribution. To simplify, let \(\small \mu=0,\sigma^2=1\) and discretize time into \(\small n\) equal parts, and it can be expressed in the form:

\[\small B(t)=\lim_{n\rightarrow \infty} \sqrt{\frac{t}{n}} \sum_{i=1}^n w_i =\lim_{\Delta t\rightarrow 0} \sum_{i=1}^n w_i \sqrt{\Delta t}. \]

By using this type of discrete approximation, Brownian motion can also be simulated on a computer. It is easy to misunderstand, but note that

\[ \small B(t) = \frac{t}{n}\sum_{i=1}^n w_i = \sum_{i=1}^n w_i \Delta t \]

does not become Brownian motion because it converges to 0 when \(\small n\rightarrow \infty\).

Multivariate Normal Random Numbers

When performing a Monte Carlo simulation on multiple random variables, if the random variables are independent and uncorrelated with each other, random numbers can be generated individually to perform the simulation. However, when dealing with multiple random variables, it is rare that they can be treated as independent variables with no correlation. Here, we consider a method for generating random numbers that follow a multivariate normal distribution that reflects the correlation coefficient.

 Suppose the correlation coefficient of \(\small m\) normally distributed random variables \(\small x_1,x_2,\cdots,x_m\) is given by \(\small \rho_{ij},i,j=1,\cdots,m\). When expressed in matrix form:

\[ \small \Sigma = \left[ \begin{array}{cccc} 1&\rho_{12}&\cdots&\rho_{1m} \\ \rho_{21}&1&\cdots&\rho_{2m} \\ \vdots&\vdots&\ddots&\vdots \\ \rho_{m1}&\rho_{m2}&\cdots&1 \end{array} \right], \]

it is called a correlation matrix, and the \(\small m\)-dimensional standard normal distribution of a correlation matrix \(\small \Sigma\) is expressed as \(\small N_{m,\Sigma}(0,1)\). What we want to do is consider an algorithm that generates \(\small m\) random numbers that follow this probability distribution.

 It is known that this can be generated by calculating a matrix decomposition called the Cholesky decomposition. The Cholesky matrix is ​​a lower triangular matrix that satisfies

\[ \small \begin{align*} &\Sigma = LL^{T} \ &L = \left[ \begin{array}{cccc} l_{11}&0&\cdots&0 \\ l_{21}&l_{22}&\cdots&0 \\ \vdots&\vdots&\ddots&\vdots \\ l_{m1}&l_{m2}&\cdots&l_{mm} \end{array} \right]. \end{align*} \]

The Cholesky decomposition can be computed starting from the previous column, \(\small j\), using the following algorithm:

\[ \small \begin{align*} &l_{jj} = \sqrt{ 1 -\sum_{k=1}^{j-1} l_{jk}^2} \qquad \qquad \quad \;\; \\ &l_{ij} = \frac{1}{l_{jj}} \left( \rho_{ij} -\sum_{k=1}^{j-1} l_{ik} l_{jk} \right), \;\; i > j. \end{align*} \]

If we calculate

\[ \small x = \left[ \begin{array}{c} x_1 \\ x_2 \\ \vdots \\ x_m \end{array} \right] = \left[ \begin{array}{cccc} l_{11}&0&\cdots&0 \\ l_{21}&l_{22}&\cdots&0 \\ \vdots&\vdots&\ddots&\vdots \\ l_{m1}&l_{m2}&\cdots&l_{mm} \end{array} \right] \left[ \begin{array}{c} w_1 \\ w_2 \\ \vdots \\ w_m \end{array} \right] \]

using independent standard normal distributed random numbers \(\small w_1,w_2,\cdots, w_m\) and the Cholesky decomposition \(\small L\) of the correlation matrix, then \(\small x_1,x_2,\cdots,x_m\) will become random numbers that follow a \(\small m\)-dimensional standard normal distribution with correlation matrix \(\small \Sigma\). By repeating this process, it is possible to perform a Monte Carlo simulation of a multivariate normal distribution.

Example: Texas Hold’em

Those who are used to draw poker may have some time to get used to good betting strategies in Texas Hold’em, because estimating the odds of winning in Texas Hold’em compared to draw poker is a much different game. For example, if you can get a straight or a flush, you might think it’s easy to win because it’s a strong hand, but if you are using four community cards, there’s a high chance that other players will also get a flush or straight, so your chances of winning may be lower. On the other hand, if you have a pair of cards (One Pair) and the common cards are all different, you may have a very high chance of winning. Just getting used to estimating the probability of winning may make you a fairly strong player.

 However, if you just play blindly, it may not be easy to get a sense of the probability of winning. In my previous post, I explained that in Texas Hold’em, when there are four players, there are

\[ \small {}_{52}C_5\times{}_{47}C_2\times{}_{45}C_2\times{}_{43}C_2\times{}_{41}C_2 = 2.0595 \times 10^{18} \]

possible card combinations. For example, at the pre-flop stage there are

\[ \small {}_{50}C_5\times{}_{48}C_2\times{}_{46}C_2\times{}_{44}C_2 = 2.34 \times 10^{15} \]

combinations, and at the flop stage there are

\[ \small {}_{47}C_2\times{}_{45}C_2\times{}_{43}C_2\times{}_{41}C_2 = 7.924 \times 10^{11} \]

combinations. Therefore, it can be assumed that it is very difficult to calculate the exact probability of winning. You can see that trying to learn the probabilities intuitively by playing repeatedly is a very inefficient approach.

 However, even without calculating all combinations, if you randomly generate cards using Monte Carlo simulation and run the simulation a certain number of times, you can expect the probability to converge to a rough estimate. When actually trying it, it appears that the probability converges after performing simulations between 10,000 and 100,000 times. Therefore, it is relatively easy to estimate the probability of winning from the hand and community cards, even without calculating all combinations. To develop a sense for estimating the probability of winning based on combinations of cards in your hand and community cards, it may be much more efficient to actually look at and memorize the probabilities calculated in a simulation. The actual implemented program will be posted on the following page, so if you would like to try it, please do so.

Comments