【Python】株価をランダムウォークでシミュレートしてみる

前回はPythonを使ってランダムウォークを再現してみたので、今回は私にとっては本題である株価を再現してみたいと思います。

※python2で書いています。

はじめに

(単純)ランダムウォークだとダメ？

前回使ったランダムウォークを使えばできそうな気がしますが、無作為な数字を加えて行くため途中でマイナスの値になることがある。株価がマイナスってありえないので不都合ですよね。

例として画像で示します。以下。

最初 40 からスタートして単純ランダムウォーク(1000歩)をした結果です。-1 か +1 の確率 1/2 なのにこんなにも下がります。

上記程度の偏りだったら「-1 か +1」 を 「-1% か +1%」 と考えて計算すればマイナスの値にはなりません。

しかし、-100% を超えるとこれも駄目です。他にも当然ですが価格が大きいと変動が凄まじいことになります。

最初 22000 からスタートして -1% か +1% でランダムウォークした結果です。ありえなくはないですが、求めているものとはちょっと違いますよね。

そこで、金融市場などのモデルでよく利用されているらしい「幾何ブラウン運動」というものを使います。

幾何ブラウン運動

期待収益率やボラティリティを指定でき、ある程度望みの株価変動を再現してくれる数式。
詳しいことは ウィキペディア などで調べてください。

幾何ブラウン運動の数式

dS_t: 価格の変化量
dB_t: ブラウン運動の変化量
μ: 期待収益率
σ: ボラティリティ

上記の数式を 伊藤の公式 を使ったりしていくと S_t(価格) が求まるらしい。

解

S_t: 価格
S₀: 最初の価格
e: 自然対数の底(ネイピア数)
μ: 期待収益率
σ: ボラティリティ
B_t: ブラウン運動の値

今回はこちらの式を使用して株価をシミュレートしてみます。

株価をシミュレートする

ここから上記の幾何ブラウン運動とPythonを使って株価をシミュレートしてみたいと思います。
サードパーティライブラリは numpy と pandas を使います。

ブラウン運動

幾何ブラウン運動をする前に、ブラウン運動が使用されているので「ブラウン運動(B_tの部分)」をPythonで書いてみます。

bmotion.py

#! /usr/bin/env python
# coding: utf-8

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt


def bm(T, n, seed=0):
    ''' ブラウン運動をする関数 Brownian Motion '''
    np.random.seed(seed)
    dt = T / n
    brownian = np.random.standard_normal(n) * np.sqrt(dt)
    return brownian.cumsum()


if __name__ == '__main__':
    T = 1.0  # 期間
    n = 250  # 期間を分割する数
    plt.show(pd.Series(bm(T, n)).plot())

実行結果

実行結果を見た感じ問題なくブラウン運動ができていそうです。

ざっくり解説。

def bm(T, n, seed=0):

ブラウン運動をするために定義した関数。
引数T が期間、引数n が期間を分割する数、引数seed は疑似乱数を作るための元。

np.random.seed(seed)

引数seed で受け取った値を numpy.random.seed() に渡して疑似乱数を設定しています。
記事を書く都合上、実行ごとに結果が異なると不都合なため書いた。

dt = T / n

変数dt に T / n を代入。
やっていることは期間(T)を細かく分けているだけ。
今回は株価をシミュレートすることを想定しているため、期間(T) = 1.0年間 に 営業日が約250日あるので、分割する数(n) = 250 とした。なので、変数dt の中身は 0.004 になっている。

ちなみに ⊿t(変化量(差)) と dt(微小差) は厳密には違う用途で使われますが、私はあまり区別していません。

brownian = np.random.standard_normal(n) * np.sqrt(dt)

変数brownian にブラウン運動の増分(⊿B)を配列で代入してます。
numpy.random.standard_normal(n) で標準正規分布(平均0, 分散1)に従う無作為な数字を n個 生成しています。
numpy.sqrt(dt) を掛けているのは dt換算 しています。ヒストリカルボラティリティのときに年換算したのと同じ流れですね。
分散を n等分 していると考えると分かりやすい？かもしれない。

return brownian.cumsum()

ブラウン運動の増分(⊿B) を一つずつ足してブラウン運動にして返しています。

幾何ブラウン運動

ブラウン運動の関数ができたので、次に幾何ブラウン運動を作ります。
別ファイルを作って先程のスクリプトをインポートして利用する形で作ります。

gbmotion.py

#! /usr/bin/env python
# coding: utf-8

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import bmotion  # 上記で書いた bmotion.py スクリプト


def gbm(S0, mu, sigma, B, T, n):
    ''' 幾何ブラウン運動をする関数 Geometric Brownian Motion '''
    t = np.linspace(0, T, n)
    trend = (mu - 0.5 * sigma**2) * t
    nonsense = sigma * B
    S = S0 * np.exp(trend + nonsense)
    return S


if __name__ == '__main__':
    S0 = 100
    mu = 0.05
    sigma = 0.16
    T = 1.0
    n = 250
    B = bmotion.bm(T, n)
    result = gbm(S0, mu, sigma, B, T, n)
    plt.show(pd.Series(result).plot())

実行結果

多分問題なく幾何ブラウン運動になっていると思われます。価格の変動具合も株価のそれと似たような感じになっていますね。
ブラウン運動と折れ線の具合がほぼ同じ形状になっているが、mu や sigma の値を変えれば形状も若干変化します。

ざっくり解説。

def gbm(S0, mu, sigma, B, T, n):

幾何ブラウン運動をするために定義した関数。
引数S0 が最初の価格、mu が期待収益率、sigma がボラティリティ、B がブラウン運動、T が期間、n が分割数。

t = np.linspace(0, T, n)

変数t に 0~T を n 等分した等差数列を配列で代入。今回は 0.0~1.0 を 250 等分した等差数列なので [0.0, 0.00401606, 0.00803213, ... , 1.0] となる。

trend = (mu - 0.5 * sigma**2) * t
nonsense = sigma * B
return S0 * np.exp(trend + nonsense)

幾何ブラウン運動の解の式に当てはめて返している。
変数trend がトレンドを表すドリフト項 (μ-σ²/2)t
変数nonsense が予測不可能な部分 σB_t

if __name__ == '__main__':
    S0 = 100
    mu = 0.05
    sigma = 0.16
    T = 1.0
    n = 250
    B = bmotion.bm(T, n)
    result = gbm(S0, mu, sigma, B, T, n)
    plt.show(pd.Series(result).plot())

最初の価格、期待収益率など幾何ブラウン運動の計算に必要な値を用意している。
そして、変数result に gbm()関数 の結果を代入し、matplotlib でグラフ表示している。

おわりに

株価をシミュレートとはいったが、μ や σ を算出して過去の相場をシミュレートしても同じグラフが出てくることはまずないです。あくまで相場のでっち上げなので、自分のトレードルールで損小利大になるかの検証に使うくらいが懸命なのかなと思います。
最近はPCもハイスペックになってきているので、多量にシミュレートして過去の相場に近似しているものから乱数をある程度絞り込めるかもしれませんが、個人でやるには現実的ではなさそうです。もちろん趣味なら全然ありです。

何かの参考になれば幸いです。