"2020年"の日経225をGeometric Brownian MotionでMonte Carlo Simulationする

一言で

2019年のデータを用いて2020年の日経平均株価のシミュレーションをしたい。

概要

確率過程 $S$ が次のStochastic Differential Equation (SDE)を満たすとする。

$dS(t) = \mu S(t) dt + \sigma S(t)dW(t).$

$W(t):$ Brownian Motion, $\mu :$ Drift coef., $\sigma :$ Diffusion coef..

このSDEには次のsolutionがある。

$S(t) = S\left(0\right)e^{\left(\mu-\frac{\sigma^2}{2}\right)t+\sigma W(t)}.$

シミュレーションの際には、資産価格 $S$ がGBMに従うと仮定した上で、次の離散化した形を用いる。

$S\left(t_{i+1}\right) = S\left(t_{i}\right)e^{\left(\mu-\frac{\sigma^2}{2}\right)\left(t_{i+1}-t_i \right)+\sigma Z\sqrt{t_{i+1}-t_i}}.$
ただし、 $Z_i \sim N\left(0,1\right);i=0,1,...,T-1.$

やってみよう

ライブラリなど

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

ASSET = '^N225'
START = '2019-01-17'
UNTIL = '2020-12-30'

df = yf.download(ASSET, start=START, end=UNTIL)
returns = np.log(df.Close / df.Close.shift(1)).dropna()

準備

TRAIN_START = '2019-01-17'
TRAIN_UNTIL = '2019-12-30'
TEST_START = '2020-01-06'
TEST_UNTIL = '2020-12-30'

train = returns[TRAIN_START:TRAIN_UNTIL]
test = returns[TEST_START:TEST_UNTIL]

INDEX =[_iter.date() for _iter in df['Adj Close'][TRAIN_UNTIL:TEST_UNTIL].index]

N = len(test)
T = len(test)
S_0 = df['Adj Close'][TRAIN_UNTIL]
mu = train.mean()
sigma = train.std()

GBMの関数を用意

Brownian Increment $dW$ の（ $シミュレーション回数 \times N$ ）の行列を用意する。行列は各行が一つのPathに該当するので、Brownian Path $W$ を累積和で計算。

def geometricBrownianMotionSimulation(_n_simulations):

    dt = T / N
    dW = np.random.normal(scale=np.sqrt(dt),
                          size=(_n_simulations, T + 1))
    def _expCumSum(x):
        return np.exp(np.cumsum((mu - 0.5 * sigma ** 2) * dt + sigma * x, axis=1))

    S_t = S_0 * _expCumSum(dW)
    S_t[:, 0] = S_0
    S_t = np.transpose(S_t)

    return S_t

プロットする関数を用意

def _plot(_simulated_df):
    ax = _simulated_df.plot(alpha=0.1, legend=False, figsize=(20,8))
    line_1, = ax.plot(INDEX, _simulated_df.mean(axis=1), color='red')
    line_2, = ax.plot(INDEX, df['Adj Close'][TRAIN_UNTIL:TEST_UNTIL], color='blue')
    ax.set_title(f'Simulation of Nikkei 225 w/ GBM; # of Simulations : {n_simulations}')
    ax.legend((line_1, line_2), ('"Counterfactual"', 'Factual'))

結果

n_simulations = 10
simulated_df = pd.DataFrame(geometricBrownianMotionSimulation(n_simulations),
                            index=INDEX)
_plot(simulated_df)

n_simulations = 100
simulated_df = pd.DataFrame(geometricBrownianMotionSimulation(n_simulations),
                            index=INDEX)
_plot(simulated_df)

n_simulations = 1000
simulated_df = pd.DataFrame(geometricBrownianMotionSimulation(n_simulations),
                            index=INDEX)
_plot(simulated_df)

n_simulations = 10000
simulated_df = pd.DataFrame(geometricBrownianMotionSimulation(n_simulations),
                            index=INDEX)
_plot(simulated_df)

感想/展望

３月以降が面白い
髪の毛みたい。