乱数の平均で答える
金融の多くの問いにはきれいな式がありません。ポートフォリオが1年で20%失う確率や、ペイオフが経路全体に依存するオプションの価格などです。モンテカルロ法は力技で答えます。求めたい量を平均 E[f(X)] として書き、X を多数シミュレーションし、それぞれに f を適用して平均を取ります。大数の法則により、標本平均は真の期待値に収束します。確率は「はい／いいえ」の指示関数の期待値にすぎないので、「イベントが起きたシミュレーションの割合」がその推定値になります。
import math
import numpy as np
# Estimate pi: the quarter circle x^2 + y^2 <= 1 covers pi/4 of the unit square
rng = np.random.default_rng(42)
x, y = rng.random(100000), rng.random(100000)
inside = x ** 2 + y ** 2 <= 1 # True/False for every point
print(round(4 * inside.mean(), 4)) # 3.1378 close to 3.1416, with 100,000 points価格と確率のシミュレーション
株をシミュレーションするには、標準正規のショックを引き、前のトピックの対数正規の式に通します。n 個の最終価格が1つの配列式で得られます。あとは価格分布に関するどんな問いも、配列に対する平均になります。80未満で終わる確率は np.mean(prices < 80)、期待価格は prices.mean()、損失の分位点は np.quantile(...) です。生成器 rng は関数の中で作らず引数として渡し、呼び出し側が再現性を管理できるようにします。
def simulate_terminal_prices(rng, s0, mu, sigma, t, n):
z = rng.standard_normal(n)
return s0 * np.exp((mu - 0.5 * sigma ** 2) * t + sigma * math.sqrt(t) * z)
prices = simulate_terminal_prices(np.random.default_rng(1), 100, 0.05, 0.2, 1, 200000)
print(round(float(np.mean(prices < 80)), 4)) # 0.1027 P(price below 80 after a year), about 10%
print(round(float(np.mean(prices > 130)), 4)) # 0.1216 P(price above 130)オプションの価格付け
ヨーロピアン・コールは満期に max(S - K, 0) を支払います。リスク中立の世界ではすべての資産が無リスク金利 rate で成長するので、公正価格は割り引いた期待ペイオフです。mu = rate で最終価格をシミュレーションし、それぞれのペイオフを計算して平均し、exp(-rate * t) を掛けます。この単純なオプションには閉じた式（ブラック＝ショールズ、下の数字では10.4506）があり、シミュレーションはそれに近づきます。モンテカルロの価値は、式のないペイオフにも同じ数行がそのまま使えることです。
def mc_call_price(rng, s0, strike, rate, sigma, t, n):
prices = simulate_terminal_prices(rng, s0, rate, sigma, t, n) # risk-neutral: the price grows at the risk-free rate
payoff = np.maximum(prices - strike, 0) # max(S - K, 0) for every simulated path
return math.exp(-rate * t) * payoff.mean() # discounted average payoff
price = mc_call_price(np.random.default_rng(1), 100, 100, 0.05, 0.2, 1, 200000)
print(round(price, 4)) # 10.3993 the exact Black-Scholes price is 10.4506: the difference is Monte Carlo error# More points, less error (1 / sqrt(n)): the estimate of pi with the same seed
for n in (100, 10000, 1000000):
g = np.random.default_rng(42)
x, y = g.random(n), g.random(n)
print(n, round(4 * float((x * x + y * y <= 1).mean()), 4))
# 100 3.36
# 10000 3.1488
# 1000000 3.1424- モンテカルロ：量を平均
E[f(X)]と書き、多数の抽出をシミュレーションし、fを適用して平均を取る - 確率は比較の平均：
np.mean(prices < 80)。生成器は引数で渡し、実行を再現可能にする - オプション価格 =
exp(-rate * t) * mean(max(S - K, 0))（価格はmu = rateでシミュレーション）。誤差も報告する
次はモンテカルロ誤差と収束です。推定がどれだけ正確か、少ない計算でより良い精度を得る方法を学びます。
ラボ：採点される内容
estimate_piで正方形にダーツを投げて円周率を求める25点estimate_probabilityでイベントの発生頻度を数える25点simulate_terminal_pricesで株価がどこに行き着くかを模擬する25点mc_call_priceでシミュレーションによりオプションを価格付けする25点