【前提】株価のモデル「幾何ブラウン運動」
突然ですが、幾何ブラウン運動とは以下の確率微分方程式
$$
dS(t) = \alpha S(t) dt + \sigma S(t) dW(t)
$$
に従う確率過程$S(t)$のことを指します。これは通常の指数型微分方程式$dS(t) = \alpha S(t)dt$に不確定性を表すブラウン運動$W(t)$の要素が足しあわされた確率過程になっており、株価のモデルとして用いられています。なおこのモデルは、時間$\Delta t$の間に一定の割合で株価が変動すると仮定し、$\Delta t\rightarrow0$とすることで導出可能です。この式の中で、$\alpha$は平均成長率、$\sigma$はボラティリティ(不確定性の度合い)を表します。
本記事の目的
突然幾何ブラウン運動について語り始めてすみませんでした。最近金融工学を勉強していたところ幾何ブラウン運動が代表的な株価モデルと登場したので、実際の現象に適用してみたいと思いました。そこで今回は、現実の株価の代表例である日経平均株価に幾何ブラウン運動をあてはめ、その平均成長率$\alpha$とボラティリティ$\sigma$を同定してやろうと思います。
日経平均株価のデータ元
今回使用するデータはこちらの日経平均プロファイルからダウンロードできます。この中の日次データをダウンロードし、一番最初と最後の行が文字化けしていたので、最初の行は列名 Date,End,Start,Max,Min を入れ、最後の行は削除しました。
今回は時刻$t$を最初の行からの経過日数、その日の株価$S(t)$を終値Endとみなし、幾何ブラウン運動へのフィッティングを行いました。
最尤推定
$\alpha, \sigma$の推定に最尤推定法を用いるので、ここではその数学的議論を行います。
伊藤の公式により、以下のような関係式を得ます。
$$
d(\log S(t)) = \left(\alpha - \frac{1}{2}\sigma^2 \right)dt + \sigma dW(t)
$$
ちなみにこの$d(\log S(t))$のことを、対数リターンというらしいです。これを離散データ$t=[t_0, t_1, ..., t_N],\ S = \left[S_0, S_1, ..., S_N\right]$に対して、$\Delta t_i = t_{i+1} - t_i$を用いて書き直すと
$$
R_i = \log\frac{S_{i+1}}{S_i} = \left(\alpha - \frac{1}{2}\sigma^2 \right)\Delta t_i + \sigma\left(W(t_{i+1}) - W(t_i)\right)
$$
となります。簡単のため$\mu = \alpha - \frac{1}{2}\sigma^2$としておきます。ブラウン運動の定義より、$W(t_{i+1}) - W(t_i)\sim\mathcal{N}(0, \Delta t_i) $ですから、対数リターン全体としては
$$
R_i \sim \mathcal{N}\left(\mu\Delta t_i,\ \sigma^2 \Delta t_i\right)
$$
となります。従って$R_i$の確率密度は、$i=0, 1, ..., N-1$について
$$
f\left(R_i|\alpha, \sigma^2\right) = \frac{1}{\sqrt{2\pi\sigma^2 \Delta t_i}} \exp\left(-\frac{\left(R_i - \mu\Delta t_i\right)^2}{2\sigma^2 \Delta t_i}\right)
$$
となります。
従って尤度関数は
$$
L\left(\mu, \sigma^2\right) = \prod_{i=0}^{N-1}f\left(R_i|\alpha, \sigma^2\right)
$$
となり、対数尤度関数は
$$
\log L\left(\mu, \sigma^2\right) = -\frac{1}{2}\sum_{i=0}^{N-1}\left(\log(2\pi\sigma^2\Delta t_i) + \frac{\left(R_i - \mu\Delta t_i\right)^2}{2\sigma^2 \Delta t_i}\right)
$$
となります。これを最大化する$\hat{\mu}, \hat{\sigma}^2$が最尤推定値なので、それぞれのパラメタで対数尤度関数を偏微分したものが0であるとして
$$
\hat{\mu} = \frac{\sum_{i=0}^{N-1} R_i}{\sum_{i=0}^{N-1}\Delta t_i}
$$
$$
\hat{\sigma}^2 = \frac{1}{n}\sum_{i=0}^{N-1} \frac{(R_i - \hat{\mu}\Delta t_i)^2}{\Delta t_i}
$$
と求まり、$\hat{\alpha}$については
$$
\hat{\alpha} = \hat{\mu} + \frac{1}{2}\hat{\sigma}^2
$$
と求まります。
ソースコード
以上の議論を日経平均株価に対して適用したものが以下のコードです。jp_ave_stock.pyとnikkei_stock_average_daily_jp.csvを同じディレクトリに配置して実行すると動きます。
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import scipy.stats as stats
# 日経平均株価シートのロード
df = pd.read_csv("nikkei_stock_average_daily_jp.csv")
df['Date'] = pd.to_datetime(df['Date'], format="%Y/%m/%d")
start_date = df['Date'].iloc[0]
df['ElapsedDays'] = (df['Date'] - start_date).dt.days
elapsed_days = df['ElapsedDays'].to_numpy()
end_values = df['End'].to_numpy()
# 対数リターンを計算
log_returns = np.diff(np.log(end_values))
delta_t = np.diff(elapsed_days)
n = len(delta_t)
# パラメータの最尤推定
mu_hat = np.sum(log_returns) / np.sum(delta_t)
sigma_hat = np.sqrt(np.sum((log_returns - mu_hat * delta_t)**2 / delta_t) / n)
alpha_hat = mu_hat + 0.5 * sigma_hat**2
print(f"期待成長率(/日): {alpha_hat:.5f}, ボラティリティ(/日): {sigma_hat:.5f}")
# プロット
fig = plt.figure(figsize=(15, 6))
# ax1: 日経平均株価とその値
ax1 = fig.add_subplot(121)
ax1.plot(elapsed_days, end_values, label=r'$S(t)$')
ax1.plot(elapsed_days, end_values[0] * np.exp(alpha_hat*elapsed_days), label=r'$S(0)\exp(\alpha t)$')
ax1.set_title(rf"ASV(Average Stock Value) vs. exponential fitting ($\alpha$={alpha_hat:.5f})")
ax1.set_xlabel("elapsed days counted from 2022/1/4")
ax1.set_ylabel("Japanese ASV (Yen)")
ax1.legend()
ax1.grid(True)
# ax2: ヒストグラムと正規分布のフィッティング
ax2 = fig.add_subplot(122)
count, bins, ignored = plt.hist(log_returns, bins=30, density=True, alpha=0.6, color='skyblue', label=r'$\log(S_{i+1}/S_i)$')
x = np.linspace(min(log_returns), max(log_returns), 1000)
pdf = stats.norm.pdf(x, loc=mu_hat, scale=sigma_hat * np.sqrt(delta_t.mean()))
ax2.plot(x, pdf, 'r-', lw=2, label=r'$N(\mu, \sigma^2)$')
ax2.set_title(rf"log returns and normal distribution ($\sigma$={sigma_hat:.5f})")
ax2.set_xlabel("log returns")
ax2.set_ylabel("probability density")
ax2.legend()
ax2.grid(True)
plt.show()
解析結果
以上のコードを実行して得られた結果を図に示します。
| 変数名 | 値 |
|---|---|
| $\hat{\alpha}$ | $2.8\times 10^{-4}$ |
| $\hat{\sigma}$ | $1.2\times 10^{-2}$ |
$\alpha$が一日あたりの平均成長率であることを考えると、一年250営業日であると考えて、年間の平均成長率はおおよそ$2.8\times10^{-4}\times 250 \times 100 = 7$ %と見積もられます ($(1+\hat{\alpha})^{250}\simeq1+250\hat{\alpha}$を用いた)。これがどれほどの値なのか自分にはわからないのですが、投資信託において年利5%で運用するのが割と現実的ラインだと聞いたことある気がするので、オーダーはそこまでずれてはいないと思います。図の左側を見てもそれなりに株価の波形に乗っかっているので、きちんと推定で来ていると言っていいでしょう。ちなみに各銀行の円預金金利を調べてみるとおおよそ0.25%らしいです。結構差があることがわかりますね。
$\sigma$についても図の右側の正規分布をみてみると、尤もらしい値が導かれているのがわかると思います。$\Delta t_i$にばらつきがある(平日の値しか記録されていない)ことから、赤線は$\bar{\Delta t} = \sum_{i=0}^{N-1} \Delta t_i$を代替的に用いた正規分布としましたが、そんなに外れた値でもないのだと思います。中心の指数的成長波形に対しておおよそ1%くらい、すなわち300円程度の変動が見込まれる、という結論になりますね。
まとめ
今回は幾何ブラウン運動を実際の日経平均株価にあてはめ、その平均成長率・ボラティリティを同定しました。簡単なプログラムでうまい具合にフィットしてくれたので、とても面白かったです。引き続き金融工学の勉強を進めていきたいと思います。