金融市場におけるモンテカルロシミュレーション

学習記録
2020.12.08

　こんにちは。システムトレーダーの卵ことKenKenです。今回も「Python3ではじめるシステムトレード ──環境構築と売買戦略」で勉強したことを簡単にまとめたいと思います。今回のテーマは、「モンテカルロシミュレーション」です！

モンテカルロシミュレーション

　モンテカルロシミュレーションとは、独立同一の分布に従う確率変数の和から、人工的に時系列を生成し、その性質を理解することです。今回は、以下の3つのモデルのうち、2つ目まで検証してみます。

  1. 単純ランダムウォーク：\(P_t=R_{t-1}+\sigma w_t\)
  2. AR(1)モデル：\(P_t=\alpha+\beta P_{t-1}+\sigma w_t\)
  3. ベルヌーイ過程：\(P_t=P_{t-1}+\sigma B_t\)

　ここで、\(P_t\)は\(t\)時の価格。\(P_{t-1}\)は\(t\)よりも１単位時間前の価格。\(w_t\)は平均ゼロ、分散１の正規分布に従う確率変数で、\(\sigma\)は定数である。\(B_t\)は+1, -1からなるベルヌーイ過程。今回のシミュレーションでは、\(\sigma=0.1/\sqrt{250}=0.00632\)としている。

1. ランダムウォーク

　まずは、１年間の\(P_t\)の動きを人工的に生成してみる。株価の初期値Pを１として、株式市場の１年間の営業日数を250日として、株価の動きを生成してみる。

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

sigma = 0.1 / np.sqrt(250)
P = [1]
for i in range(1, 250):
    w = np.random.normal(0, 1)
    P.append(P[i-1] + sigma * w)

これで価格の時系列を生成することができた。次に、価格の差分をとり、その平均、標準偏差、歪度、尖度を確かめる。また、それに合わせて、グラフと価格変化のヒストグラムを書いてみる。

price = pd.Series(P)
diff_P = price.diff().dropna() # 価格の変化(差)
# 平均、標準偏差、歪度、尖度を表示
print('mean:{0:.4f}, std:{1:.4f}, skew:{2:.4f}, kurt:{3:.4f}'.format(diff_P.mean(),diff_P.std(), diff_P.skew(), diff_P.kurt()))

# グラフを書く
plt.figure(figsize=(8, 3))
ax1 = plt.subplot(1, 2, 1)
price.plot(color='darkgray', ax=ax1)
plt.xlabel('t')
plt.ylabel('$P_t$')

# ヒストグラムを書く
ax2 = plt.subplot(1, 2, 2)
diff_P.hist(color='lightgreen', ax=ax2)
plt.xlabel('P_t-P_{t-1}')
plt.ylabel('frequency')

　平均はゼロで標準偏差もほぼ想定通りで、歪度と尖度はほぼゼロに近い。ヒストグラムを見てみると、比較的綺麗な正規分布が描かれた。
　次に、250個の価格データを生成した上述の過程を10000回繰り返し、1日目～250日目の最大値、最小値を求め、その推移をグラフに書いてみる。それにより、価格の推移の期間構造を求めてみる。この作業は、前回の記事で日経平均株価に対して行ったものと同じである。詳細は前回の記事をご参照ください。

# 10000回シミュレーションを実施
P_250 = []
P_diff = []
high = [0]*250
low = [1]*250
for j in range(10000):
    P0 = 1
    for i in range(250):
        w = np.random.normal(0, 1)
        P0 = P0 + sigma * w
        if P0 > high[i]:
            high[i] = P0
        if P0 < low[i]:
            low[i] = P0
        P_diff.append(sigma * w)
    P_250.append(P0)

# グラフを書く
plt.figure(figsize=(9, 4))
# 推移
high = pd.Series(high)
low = pd.Series(low)
ax1 = plt.subplot(1, 2, 1)
high.plot(label='high', linestyle='--', ax=ax1)
low.plot(label='low', linestyle='--', ax=ax1)
plt.title('max$P_t$-min$P_t$')
plt.xlabel('x')
plt.ylabel('$P_t$')
plt.legend(loc='upper left')

# 250日目の価格のヒストグラム
ax2 = plt.subplot(1, 2, 2)
P_250 = pd.Series(P_250)
P_250.hist(bins=100)
plt.xlabel('$P_{250}$')
plt.ylabel('frequency')

# 平均、標準偏差、歪度、尖度を表示
print('mean:{0:.4f}, std:{1:.4f}, skew:{2:.4f}, kurt:{3:.4f}'.format(P_250.mean(),P_250.std(), P_250.skew(), P_250.kurt()))

　左図から価格の最大値と最小値はほぼ直線的に伸びていることがわかる。右図から、分布はほぼベル型であり、250日目の価格は設定どおり、平均１、標準偏差ゼロの分布にしたがっている。歪度と尖度もほぼゼロとなっている。
　前回、実際に価格の動きを同様の手順で調べてみた。その時の価格推移は以下であった。

　最小値(low)のグラフは割とスムーズだが、最大値(high)のグラフは凹凸した形をしている。パラメータが一定のシミュレーションでは、現実世界を説明はできないが、似たような動きをしていることは見てとれる。

　次に、得られた時系列から1日間隔の変化率とそれを１日ずらした変化率の散布図を書いてみる。

# 散布図
P_diff = pd.Series(P_diff)
plt.figure(figsize=(6, 6))
plt.scatter(P_diff, P_diff.shift(1), alpha=0.01)
plt.xlabel('$P_{t-1}/P_{t-2}-1$')
plt.ylabel('$P_{t}/P_{t-1}-1$')

実際のデータとは異なり、丸の円が描かれた。つまり、二つの変化率に相関が無いことがわかる。これが分散不均一性もボラティリティクラスタリングもない世界である。

AR(1)モデル

　1次の自己回帰モデルは以下の式で与えられる。

$$P_t=\alpha+\beta P_{t-1}+\sigma w_t$$

この\(\beta\)を指定することで、AR(1)がどのような性質を持っているか調べてみる。\(\beta\)を指定することで\(\alpha\)を決める必要があるが、これは期待値(\(\mu\)から求めることができる。AR(1)の期待値\(E(P_t)\)は以下で与えられる。

$$E(P_t) = E(\alpha) + E(\beta P_{t-1}) + E(\sigma w_t)$$

ここで、\(E(P_t) = \mu \)とすると、\(E(\alpha)=\alpha \)、\(E(w_t)=0\)であるから、
$$ \mu = \alpha + \beta \mu = \frac{\alpha}{1-\beta}$$

となる。価格の初期値を1とすると、\(\alpha = 1 – \beta \)である。分散は、

$$var(P_t) = E(P_t – \mu)^2=E(P_t^2) – \mu^2$$

であり、\(\alpha=0\)であれば、

$$var(P_t) = \frac{\sigma_z^2}{1-\beta^2}$$

である。つまり期待値も分散も時間と独立している。ここがランダムウォークとは異なる点となる。
　以降の分析のために、以下のようにAR(1)の時系列を生成する関数を定義しておく。

# AR(1)時系列生成関数を定義
def AR1(beta, sigma, n, m, p0):
    '''
    AR(1)における時系列データを生成
    beta:回帰係数, sigma:標準偏差, n:生成するデータ数, m:分析期間, p0：初期値 
    '''
    P = []
    dP = []
    high = [0] * m
    low = [p0] * m
    alpha = (1-beta)*p0
    sigma_w = sigma * p0
    for j in range(n):
        P0 = p0
        for i in range(m):
            w = np.random.normal(0, 1)
            P1 = beta * P0 + alpha + sigma_w*w
            dp = P1 - P0
            P0 = P1
            if P0 > high[i]:
                high[i] =P0
            if P0 < low[i]:
                low[i] =P0
            dP.append(dp)
        P.append(P0)
    return pd.Series(price), pd.Series(dprice), pd.Series(high), pd.Series(low)

　まず初めに、beta=0.9999, sigma=0.00632, n=10000, m=250, p0=1として時系列を生成し、最大値・最小値の推移やヒストグラムを表示する。

# 時系列を生成
beta, sigma, n, m, p0 = (0.9999, 0.00632, 10000, 250, 1) # パラメータ設定
price, dprice, high, low = AR1(beta, sigma, n, m, p0)

# グラフを書く
plt.figure(figsize=(8, 4))
# 最大値と最小値の推移
ax1 = plt.subplot(1, 2, 1)
high.plot(label='high', linestyle='--', ax=ax1)
low.plot(label='low', linestyle='--', ax=ax1)
plt.title('max$P_t$-min$P_t$')
plt.xlabel('x')
plt.ylabel('$P_t$')
plt.legend(loc='upper left')

# 250日目の価格のヒストグラム
ax2 = plt.subplot(1, 2, 2)
price.hist(bins=100)
plt.xlabel('$P_{250}$')
plt.ylabel('frequency')

# 平均、標準偏差、歪度、尖度を表示
print('mean:{0:.4f}, std:{1:.4f}, skew:{2:.4f}, kurt:{3:.4f}'.format(dprice.mean(), dprice.std(), dprice.skew(), dprice.kurt()))

　価格差の平均はゼロ、標準偏差は0.0063、歪度、尖度ともにゼロ近辺という結果になった。\(\beta<1\)である特徴がわずかに出ているかもしれないが、ランダムウォークとほぼ同等である。次に、1単位時間差について散布図を書いてみてみる。

# 散布図
plt.figure(figsize=(6, 6))
plt.scatter(dprice, dprice.shift(1), alpha=0.01)
plt.xlabel('$P_{t-1}/P_{t-2}-1$')
plt.ylabel('$P_{t}/P_{t-1}-1$')

価格差は綺麗な円形を示している。つまり、自己相関がほぼない。

　次に、AR(1)の特徴をさらに明確につかむために、\(\beta\)の値を動かしてみる。今回は、\(\beta=0.99, 0.5, 1.003, 1.005, 1.01\)として実験をしてみる。

\(\beta=0.99\)
# 時系列を生成
beta, sigma, n, m, p0 = (0.99, 0.00632, 10000, 250, 1) # パラメータ設定
price, dprice, high, low = AR1(beta, sigma, n, m, p0
# 以下省略 #

結果は、価格差の平均はゼロ、標準偏差は0.0063、歪度と尖度はほぼゼロとなった。\(\beta=0.99\)では、ある日数を経過と、価格の最大値と最小値の幅が一定値にとどまる様子が見られる。これが、中心回帰、定常確率過程の特徴らしい。それでも定常性を得るために、50日程度かかっていることがわかる。

\(\beta=0.5\)

　次に、\(\beta=0.5\)にしてみる。

# 時系列を生成
beta, sigma, n, m, p0 = (0.5, 0.00632, 10000, 250, 1) # パラメータ設定
price, dprice, high, low = AR1(beta, sigma, n, m, p0)
# 以下省略 #

結果は、価格差の平均はゼロ、標準偏差は0.0063、歪度と尖度はほぼゼロとなった。\(\beta=0.5\)ではほぼ主運管的に中心回帰の定常確率過程の特徴が見えてくる。

\(\beta=1.003\)

次に、AR(1)が発散する場合(\(\beta>1\)として\(\beta=1.003\)を見てみる。

# 時系列を生成
beta, sigma, n, m, p0 = (1.003, 0.00632, 10000, 250, 1) # パラメータ設定
price, dprice, high, low = AR1(beta, sigma, n, m, p0)
# 以下省略 #

　結果は、価格差の平均はゼロ、標準偏差は0.0063、歪度と尖度はほぼゼロとなった。\(\beta=1.003\)ではまだ発散に至っておらず、ランダムウォークとの区別が難しい。

\(\beta=1.005\)

　\(\beta=1.005\)にしてみる。

# 時系列を生成
beta, sigma, n, m, p0 = (1.005, 0.00632, 10000, 250, 1) # パラメータ設定
price, dprice, high, low = AR1(beta, sigma, n, m, p0)
# 以下省略 #

結果は、価格差の平均はゼロ、標準偏差は0.0063、歪度と尖度はほぼゼロとなった。\(\beta=1.005\)では発散、の一歩手前である様子が見られる。

\(\beta=1.01\)

最後に\(\beta=1.01\)にしてみる。

# 時系列を生成
beta, sigma, n, m, p0 = (1.01, 0.00632, 10000, 250, 1) # パラメータ設定
price, dprice, high, low = AR1(beta, sigma, n, m, p0)
# 以下省略 #S

結果は、価格差の平均はゼロ、標準偏差は0.0063、歪度と尖度はほぼゼロとなった。\(\beta=1.01\)では発散一歩の様子を見ることができた。

AR(1)まとめ

　AR(1)モデルについてここまで(\beta\)を何通りか動かしシミュレートしてみた。\(\beta>1\)であってもその程度が小さければ、価格過程は破綻することが無いことがわかった。また、それは観測期間にも左右される様子が見られた。期間が短ければ、\(\beta\)が1よりも大きくても破綻せずに済む確率は上がる。また、\(\beta<1\)近辺では、ランダムウォークとの区別が難しい結果となった。

　自己回帰モデルの最大の特徴は、価格の中心回帰の動きである。ランダムウォークの価格の動きは時間の経過と共に分散が大きくなるのに対して、自己回帰モデルではある点に収束する。その為、将来予測を行うためには便利である。つまり、価格の変化幅が一定に収まるためってことなのかな。。。

ベルヌーイ過程

　最後にベルヌーイ過程についてシミュレートしてみる。結果が二つしかない実験をベルヌーイ試行と呼び、とびとびの値からなる独立な離散時間の確率変数列をベルヌーイ過程と呼ぶ。今回は、とりうる二つの値を{-5, +5}として1万個の確率変数を発生させ、その和を一日分の取引とする。
　まずは、以下の通り、生成器を関数で定義し、試しに価格の推移をグラフにしてみた。

#ベルヌーイ過程に基づく価格を生成
def bernoulli(N, p0):
    #N：個数、p0：価格の初期値
    s = (np.random.binomial(1, 0.5, N)-0.5) * 10 # -0.5 or 0.5の乱数をn個生成後、10を書けることで{-5, +5}にしている
    P = [p0]
    for i in range(N):
        P.append(P[i] + s[i])
    return pd.Series(P)

P = bernoulli(10000, 0)
# 価格推移を書く
plt.figure(figsize=(8, 4))
P.plot(color='darkgray')
plt.xlabel('step $t$')
plt.ylabel('$\sum_{1}^{t} B_t \cdot 5$')

　実際の価格の動きに似たものを再現できている。これはランダムウォークであり、一定方向に価格が上昇したり下落しているのは確率的トレンドである。

　次にこの実験を10000回繰り返し、それぞれの試行で生成した確率変数の和をヒストグラムにしてみる。

N = 10000
M = 10000
Q = []
for j in range(M):
    P =bernoulli(N, 0)
    Q.append(P[N])
plt.figure()
plt.hist(Q, density=True, bins=25)
plt.xlabel('$\sum_{1}^{t} B_t \cdot 5$')
plt.ylabel('probility deinsity function')

結果は、正規分布に近いベル型の形となった。次に1日の取引回数を3000回に減らしてみる。

N = 3000
M = 10000
Q = []
for j in range(M):
    P =bernoulli(N, 0)
    Q.append(P[N])
plt.figure()
plt.hist(Q, density=True, bins=25)
plt.xlabel('$\sum_{1}^{t} B_t \cdot 5$')
plt.ylabel('probility deinsity function')

結果として、分布の幅が狭くなっていることがわかる。つまり、一日の価格の動きが取引回数に依存していることがわかる。また、先ほどに比べ分布の滑らかさが維持されていないことにも気づく。つまり、ボラティリティは取引数と関係があることがシミュレーションを通してもわかる。

まとめ

　今回は、単純ランダムウォーク、ARモデル、ベルヌーイ過程の３つのモデルについてモンテカルロシミュレーションを試してみました。単純ランダムウォークのシミュレーションでは、実際の日経平均株価と比較し、現実の価格推移と似たような動きをすることが確認できました。ARモデルでは、\(\beta\)の値をいろいろ動かしてみてその特性についてみてみました。回帰係数からそのモデルの特徴を掴むのは難しいですが、シミュレーションを通して視覚的にARモデルについて理解を深めることができたのではないかと思います。この点がモンテカルロシミュレーションを行ってみるメリットの一つかもしれないですね。ベルヌーイ過程では、シミュレーションを通してボラティリティと取引数の関係性についてみることができました。

　仮定が成立しない現実世界ではシミュレーション通りになることは滅多にないと思います。しかし、仮説が成立した場合にモデルがどのような特性を持つか知っておくことは非常に大切なことだと思います。モデルの特性を理解しておくことで分析結果に対してより深い考察ができるようになるのではないかと思います。今後のシストレ開発でも、単にリアルデータに対して検証するに留まらず、モンテカルロシミュレーションを取り入れていきたいと思います。

以上

タイトルとURLをコピーしました