用 Python 學蒙地卡羅方法/模擬

如果說有一種方法能像奇異博士那樣看到 14,000,605 種未來，那大概就是蒙地卡羅方法了。

在接觸蒙地卡羅方法之前，我曾以為它是一種神奇的演算法。但深入了解後才發現，蒙地卡羅方法其實是一種簡單、直觀且實用的工具！

本文將介紹蒙地卡羅方法，並使用 Python 進行實際演練，從計算 π 值到模擬台股的長期投資，幫助大家從理論到實務全面掌握！

本文環境 #

• Python 3
• numpy

蒙地卡羅(Monte Carlo method)方法簡介 #

蒙地卡羅方法，也稱為蒙地卡羅模擬(Monte Carlo Simulation)，是一種基於機率與統計理論的計算方法，常用於預測具有不確定性事件的可能結果。

雖然名稱聽起來或許有些神秘，但蒙地卡羅方法的核心其實是利用隨機數或特定的機率分布，通過大量的模擬實驗來探究各種可能性，從而全面了解潛在的結果範圍。

使用蒙地卡羅方法時通常需要 2 個步驟（又或者說滿足以下 2 個條件，就等同使用蒙地卡羅方法）：

1. 定義事件的機率分佈：針對要模擬的事件，設定其隨機變數的機率分佈（例如常態分佈、均勻分佈等）。
2. 進行多次隨機模擬：基於所定義的機率分佈，對事件進行大量的隨機模擬，統計各種可能結果的出現頻率，從而預測事件的結果範圍。

舉賭博為例，假設一個賭徒本金為 100 塊，勝率為 30%，贏了可以得到 2 塊，輸了本金少 1 塊，那麼我們就可以用蒙地卡羅方法模擬 1,000 次賭徒進行 100 場賭博之後的可能結果。

在上述例子中，需要設定的機率為勝率 30%，需要模擬的事件為 100 場賭博，次數為 1,000 次，以 Python 程式碼表示的話，如下所示：

# 模擬 1,000 次
for _ in range(1_000):
    capital = 100
    round_results = [100]

    # 模擬 100 次賭博
    for _ in range(100):
        if np.random.rand() <= 0.3:
             # 贏的情況
            capital += 2
        else:
            # 輸的情況
            capital -= 1
        round_results.append(capital)

    results.append(round_results)

上述程式碼中的變數 results 會儲存 1,000 個長度為 100 的 list，每 1 個 list 都是 100 次賭博的金額變化。

我們可以把結果畫出來，看一下賭徒 1,000 次實驗的情況：

import matplotlib.pyplot as plt

plt.figure(figsize=(10, 6))
for i in range(1000):
    plt.plot(results[i], color='blue', alpha=0.03, linewidth=2)

plt.title('Monte Carlo Simulation')
plt.xlabel('Game #')
plt.ylabel('Capital ($)')
plt.show()

從結果可以看到賭徒最好的情況是累積本金到 130 左右，最糟情況是本金賠到約 40 左右，而顏色越深的部分代表多次實驗結果重疊，從結果也可以看到，對賭徒而言很可能在經過 100 場賭博之後，他的本金會低於 100。

蒙地卡羅(Monte Carlo method)方法經典案例 — 預估圓周率 π #

蒙地卡羅方法有個經典案例 — 預估圓周率 π 。

它的原理是假設我們在 1 個方形範圍中撒豆子，只要豆子足夠多、足夠均勻，我們就可以藉由豆子的面積總和估算出方形範圍內的圓面積。

假設半徑為 1 ，可以得到圓面積為 π ，而方形面積為 4，因此可以得出圓面積佔比方形面積 π/4 。

根據前述資訊，我們可以推導出：

圓面積 = 方形面積 * π / 4

而方形面積可以用落於方形區域的豆子面積總和，圓面積則可以用落於圓型區域的豆子面積總合：

圓形區域總豆子數 = 方形區域總豆子數 * π / 4

最後再左右等式移動一下：

π = 4 * 圓形區域總豆子數 / 方形區域總豆子數

如此就能夠預估出 π 的值。

以 Python 程式碼表示的話（取自 Wikepedia），如下所示：

import numpy as np

# 生成10的5次方個隨機點
num_points = 10**5
points = np.random.rand(num_points, 2)

# 計算點到原點的距離
distances = np.sqrt(points[:,0]**2 + points[:,1]**2)

# 計算落在單位圓內的點的數量
points_inside_circle = np.sum(distances <= 1)

# 蒙地卡羅方法計算圓周率
pi_estimate = 4 * points_inside_circle / num_points

print(pi_estimate)

上述程式碼中，np.random.rand(num_points, 2) 就同時完成蒙地卡羅產生隨機數值與重複實驗的過程。

蒙地卡羅(Monte Carlo method)方法應用 — 投資結果預估 #

蒙地卡羅方法可以運用在很多領域，包含金融、經濟、生物、物理（粒子、量子等等）、機器學習等等。

本文再展示 1 個跟大家較實用的應用領域 — 投資結果預估。

假設台股近 20 平均年化報酬率 8% ，報酬率的標準差為 24.5%（常態分佈需要設定標準差），我們每年固定投資 10 萬，我們就可以用蒙地卡羅方法預估 30 年後的總資產，不過股市每年漲幅不一，我們可以用常態分佈方式隨機設定每年的漲幅。

import numpy as np
import matplotlib.pyplot as plt

def monte_carlo_simulation(
    annual_investment,
    annual_return_mean,
    annual_return_std,
    years,
    num_simulations,
):
    results = []
    for _ in range(num_simulations):
        portfolio_value = 0
        yearly_values = [portfolio_value]
        for _ in range(years):
            portfolio_value += annual_investment
            annual_return = np.random.normal(
                annual_return_mean,
                annual_return_std,
            )
            portfolio_value *= (1 + annual_return)
            yearly_values.append(portfolio_value)

        results.append(yearly_values)

    return np.array(results) / 10000  # Convert to units of 10k


# Example parameters
annual_investment = 100000   # $100,000 invested every year
annual_return_mean = 0.08    # 8% mean return
annual_return_std = 0.245    # 24.5% volatility
years = 30                   # 30-year horizon
num_simulations = 1000       # 1000 simulations

# Run simulation
simulated_results = monte_carlo_simulation(
    annual_investment,
    annual_return_mean,
    annual_return_std,
    years,
    num_simulations,
)

上述程式碼的 monte_carlo_simulation() 為蒙地卡羅模擬的部分，其中 np.random.normal(annual_return_mean, annual_return_std) 會從常態分佈中抽隨機數作為該年度的報酬率，剩下的部分則是複利的計算，此處也可以按照需求改為更好的計算方式。

然後，我們可以使用 10, 25, 50, 75, 90 分位數，將 30 年後的結果視覺化呈現：

# Calculate percentiles for grouping results
percentiles = np.percentile(
    simulated_results,
     [10, 25, 50, 75, 90],
    axis=0,
)

# Plot the results grouped by percentile ranges
plt.figure(figsize=(10, 6))

# Fill between the 10th and 90th percentiles
plt.fill_between(
    range(years + 1),
    percentiles[0],
    percentiles[4],
    color='lightblue',
    alpha=0.5,
    label='10th-90th Percentile',
)

# Fill between the 25th and 75th percentiles
plt.fill_between(
    range(years + 1),
    percentiles[1],
    percentiles[3],
    color='blue',
    alpha=0.5,
    label='25th-75th Percentile',
)

# Plot the median (50th percentile)
plt.plot(percentiles[2], color='black', label='50th Percentile (Median)')

plt.title('Monte Carlo Simulation of Annual $100,000 Investments (in 10k Units)')
plt.xlabel('Years')
plt.ylabel('Portfolio Value (10k $)')
plt.legend()
plt.show()

上述程式碼執行結果如下：

從上圖可以看到，在最佳情況下，我們很可能累積超過 2,500 萬的資產，就算在 75 分位數的最好情況下，也可以累積接近 1,500 萬的資產，而一般的情況，中位數也有超過 500 萬的總資產。

總結 #

由於現實中充滿各種不確定因素，蒙地卡羅方法為我們提供了一種有效的工具來預測或預估這些不確定性對結果的影響，不僅能幫助我們更好地認識可能的最終結果，還能讓我們為最好的結果和最壞的情況做好準備。

以上！

Enjoy!

References #

蒙地卡羅方法 - 維基百科，自由的百科全書

什麼是蒙特卡羅模擬？