確率変数と期待値

まだ値が決まっていない数

確率変数は、サイコロの目、明日のリターン、満期時の株価のように、値が偶然に左右される量です。離散確率変数は、取りうる値のリストを持ち、それぞれに確率があります。値と確率をあわせたものが分布です。確率は負にならず、合計は1です。Python では、同じ長さの2つのリスト values と probs で表せます。公正なサイコロは values = [1, 2, 3, 4, 5, 6] と probs = [1/6] * 6 です。

import math

values = [1, 2, 3, 4, 5, 6]
probs = [1 / 6] * 6

# A valid distribution: no negatives, and the probabilities add up to 1
print(all(p >= 0 for p in probs))             # True
print(math.isclose(sum(probs), 1.0))          # True (compare floats with isclose)

# A stock tomorrow: down 5%, up 2% or up 8%
returns = [-0.05, 0.02, 0.08]
chances = [0.2, 0.5, 0.3]
print(math.isclose(sum(chances), 1.0))        # True

期待値：確率で重み付けした平均

期待値 E[X] は、多数回繰り返したときの平均です。各値にその確率を掛けて足します：sum(v * p)。サイコロの E[X] = 3.5 は決して出ない数ですが、これは長期平均であって、1回の結果の予測ではないので問題ありません。期待値は線形です：E[a*X + b] = a*E[X] + b。関数の期待値 E[g(X)] は、まず各値に g を適用します：sum(g(v) * p)。注意：一般に E[g(X)] は g(E[X]) と等しくありません。オプションのペイオフ max(S - K, 0) が典型例です。

def expected_value(values, probs):
    return sum(v * p for v, p in zip(values, probs))

die = [1, 2, 3, 4, 5, 6]
fair = [1 / 6] * 6
print(round(expected_value(die, fair), 4))            # 3.5

# E[g(X)]: apply g to each value first. Call payoff with strike 100:
prices = [90, 100, 110, 120]
chances = [0.1, 0.4, 0.3, 0.2]
mean_price = expected_value(prices, chances)
print(mean_price)                                      # 106.0
print(expected_value([max(s - 100, 0) for s in prices], chances))  # 7.0  = E[max(S - K, 0)]
print(max(mean_price - 100, 0))                        # 6.0  = max(E[S] - K, 0), not the same

分散：どれだけ散らばっているか

期待リターンが同じでもリスクが大きく異なる2つの投資があり得ます。分散は平均のまわりの散らばりを測ります：Var(X) = E[(X - mu)^2] で、E[X^2] - E[X]^2 としても計算できます。標準偏差はその平方根で、X と同じ単位です。リターンならボラティリティです。サイコロでは Var = 35/12 = 2.9167 です。どちらも X の関数の期待値から作れるので、expected_of(g, values, probs) を1つきちんと書けば、平均も2次モーメントもその先も得られます。

def expected_value(values, probs):
    return sum(v * p for v, p in zip(values, probs))

def expected_of(g, values, probs):
    return expected_value([g(v) for v in values], probs)

def variance(values, probs):
    mean = expected_of(lambda v: v, values, probs)
    return expected_of(lambda v: v * v, values, probs) - mean ** 2

die, fair = [1, 2, 3, 4, 5, 6], [1 / 6] * 6
print(round(variance(die, fair), 4))                   # 2.9167

returns, chances = [-0.05, 0.02, 0.08], [0.2, 0.5, 0.3]
print(round(expected_value(returns, chances), 4))      # 0.024
print(round(variance(returns, chances), 6))            # 0.002044

シミュレーションと大数の法則

サンプルを引いて期待値を推定することもできます。rng.choices(values, weights=probs, k=n) は n 個の標本を返し、n が大きくなると標本平均は E[X] に近づきます。これが大数の法則で、モンテカルロ法による価格付けの核心です。実行が再現でき、他のコードに影響しないよう、自分専用の random.Random(seed) を使いましょう。標本が少ないと平均はぶれます。誤差は n の平方根の逆数でしか減らないため、標本を100倍にしても精度は10倍にしかなりません。

import random

die, fair = [1, 2, 3, 4, 5, 6], [1 / 6] * 6

def sample_mean(values, probs, n, seed):
    rng = random.Random(seed)                    # private, repeatable generator
    draws = rng.choices(values, weights=probs, k=n)
    return sum(draws) / n

for n in (10, 1000, 100000):
    print(n, round(sample_mean(die, fair, n, 1), 4))
# 10 3.1
# 1000 3.598
# 100000 3.5017       (closing in on 3.5)

print(sample_mean(die, fair, 1000, 1) == sample_mean(die, fair, 1000, 1))  # True: same seed, same answer

• 分布は values と probs（負なし、合計1）。E[X] = sum(v * p)
• E[g(X)] = sum(g(v) * p) は g(E[X]) と異なる。Var = E[X^2] - E[X]^2
• 標本平均は E[X] に収束する（大数の法則）。自分専用の random.Random(seed) を使う

次は NumPy 配列と dtype です。コースの残りの土台となる、高速で型付きの配列です。

ラボ：採点される内容

• expected_value で結果を確率で重み付け平均する25点
• expected_of でペイオフの期待値を求める25点
• variance で散らばりを測る25点
• sample_mean でシミュレーションにより期待値を推定する25点