モンテカルロで見える世界

インターネット上から金融関連の価格データを手に入れ、そのデータにトレンドが在るか無いかを見てきました。

トレンドがあれば、確定的トレンドと確率的トレンドに分類され、また、時間との関係では定常時系列と非定常時系列に分類されました。

確定的トレンド、確率的トレンド、定常時系列、または時間トレンド、ランダムウォーク、AR(1)が時系列分析の3種の神器でした。

モンテカルロ・シミュレーションの利用  

日経平均株価の動きの3つのモデルは  

$1.P_{t} = P_{t-1} + σw_{t}$

$2.P_{t} =a + βP_{t-1} + σw_{t}$ 

$3.P_{t} = P_{t-1} + σB_{t}$

です。

$P_{t}$は t 時の価格である。

$P_{t-1}$はtよりも1単位時間前の価格です。

$w_{t}$は平均ゼロ、分散1の正規分布に従う確率変数であり、σは定数です。

$B_{t}$は+1,-1からなるベルヌーイ過程である。シミュレーションではσ=0.00632とする。  

ランダムウォーク  

ランダムウォーク・モデルは  

$P_{t} = P_{t-1} + σw_{t}$

で与えられる。

ここではα= 0としてドリフト無しランダムウォークを検討する。

σは一定である。  

まず1年間の$P_{t}$の動きを人工的に生成してみよう。

np.randomというモジュールの中にある乱数生成器normalを用いて確率変数$w_{t}$を生成する。

株価の初期値 P を1とする。

株式市場の1年の営業日数を250日とし、250個の株価の動きを生成してみよう

In [1]:
%matplotlib inline
import pandas as pd
import numpy as np
sigma=0.1/np.sqrt(250)
print(sigma)
P=[1]
for i in range(1,250):
    w=np.random.normal(0,1)
    P.append(P[i-1]+sigma*w)
0.00632455532034

これでPの時系列が生成できた。価格の差分を取り、その平均、標準偏差、歪度、尖度を確かめてみます。

In [2]:
price=pd.Series(P)
dp=price.diff().dropna()
print("mean %2.5f std %2.5f skew %2.5f kurt %2.5f"\
%(dp.mean(),dp.std(),dp.skew(),dp.kurt()))
mean -0.00006 std 0.00629 skew -0.12995 kurt 0.02848

平均はほぼゼロ、標準偏差はほぼ予定通り、歪度と尖度もほぼゼロに近い。次に結果をグラフで確かめてみます。

In [3]:
import matplotlib.pyplot as plt
plt.figure(figsize=(10,5.2))
plt.subplot(121)
price.plot(color='darkgray')
plt.xlabel('t')
plt.ylabel('$P_t$')
plt.subplot(122)
dp.hist(color='lightgreen')
mx=round(dp.max(),2)
mn=round(dp.min(),2)
plt.xticks([mn,0,mx])
plt.xlabel('$P_t-P_{t-1}$')
plt.ylabel('frequency')
Out[3]:
<matplotlib.text.Text at 0x7f751b1ea4a8>

左側のチャートは1日目の価格、2日目の価格といった具合に250日間の価格の推移を描いていて、右側の頻度図は日々の価格の変化の大きさに対する頻度を表示している。

価格の差分は1回の試行の結果としては、かなりきれいな正規分布である。

プログラムコードの解説

P=[1]

for i in range(250): 
    w=np.random.normal(0,1)

P.append(P[i-1]*(1+sigma*w))

生成された活くを格納する容器の準備を整えている

[1]は初期値を1に設定している。

for文以下は250回繰り返す

平均ゼロ、分散1の確率変数を生成する。

価格の1期間更新している。

P=[] 価格の配列の初期化

high=[0]*250 価格の上限の配列の初期化

low=[1]*250 価格の下限の配列の初期化

for j in range(10000): for文以下を10000回繰り返す

P0=100 価格の初期値を100回に設定

for i in range(250): for文以下を250回繰り返す

w=np.random.normal(0,1) 標準正規文に従う確率変数の生成

P0=P0(1+sigmaw) 価格の更新

if P0>high[i]: 上限の更新

high[i]=P0 if P0<low[i]: 下限の更新

low[i]=P0 P.append(P0) 最終日の価格をPに保存

price=pd.Series(P) リストPを Seriesに変換しpriceとして保存

次に250個の価格データを生成した上述の過程を10000回繰り返すことで、1日目、2日目の価格の最大値、最小値を250日目まで求めることで価格の推移の期間構造を求めよう

In [4]:
P=[]
dP=[]
high=[0]*250
low=[1]*250
for j in range(10000):
    P0=1
    for i in range(250): 
        w=np.random.normal(0,1)
        dp=sigma*w
        P0=P0+dp
        if P0>high[i]:
            high[i]=P0
        if P0<low[i]:
            low[i]=P0
        dP.append(dp)
    P.append(P0)
price=pd.Series(P)
dprice=pd.Series(dP)

価格の時系列はpriceに、価格差の時系列はdpriceにpandasシリーズとして格納した。確認のためにグラフにしてみます。

In [5]:
plt.figure(figsize=(9,4))
plt.subplot(1,2,1)
plt.plot(high,label="high",linestyle='--')
plt.plot(low,label="low",color='darkgray')
plt.title('max$P_t$-min$P_t$')
plt.xlabel('t')
plt.ylabel('$P_t$')
plt.legend(loc='upper left')
plt.subplot(1,2,2)
price.hist(bins=100,color='pink')
plt.xlabel('$P_{250}$')
plt.ylabel('frequency')
print("mean %2.5f std %2.5f skew %2.5f kurt %2.5f"\
%(dprice.mean(),dprice.std(),dprice.skew(),dprice.kurt()))
mean 0.00000 std 0.00633 skew 0.00095 kurt 0.00038

このモンテカルロシミュレーションでは250日間の価格をまず生成し、それを10000回繰り返すことで、時系列の特性を得ようとしました。

左の図は250日間の価格の推移の特性を調べています。

10000回の試行の内の1日目の価格の最大値と最小値をプロットし、それを2日目、3日目にも行い、250日目まで繰り返すと価格の推移の期間構造が得られます。

価格の最大値と最小値の幅がほぼ直線的(緩やかな曲線)に増加しているのが分かります。

右図は250日目の価格の頻度を表している。分布はほぼベル型です。

価格差の平均はゼロ、標準偏差は0.00633と設定値そのものである。歪度と尖度もほぼゼロです。  

このモンテカルロシミュレーションから得られた時系列の特性は実際の価格の動きの特性を説明しているであろうか? 

前回の結果を思い出してほしい

Material 9 変化率の最大値、最小値の期間構造のグラフ参照

最小値(low)のグラフは割とスムーズだが最大値(high)のグラフは凸凹である。

一方、シミュレーションではほぼ曲線的に下落、上昇している。

パラメータが一定のシミュレーションでは現実の世界を説明できないが、そんなに外れてもいない。  

次に得られた時系列から1日間隔の変化率とそれを1日ずらした変化率の散布図を描いてみよう

In [6]:
plt.figure(figsize=(5,5))
mx=round(dprice.max(),2)
mn=round(dprice.min(),2)
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$')
plt.xticks([mn,0,mx])
plt.show()

まん丸の円が描かれた。実際のデータとはだいぶ違う。2つの変化率には相関がないことが見て取れます。

これが分散不均一性もボラティリティ・クラスタリングも無い世界である。  

このように、実際の時系列を分析する際には、モデルの厳密さにこだわるのではなく、まずは大まかな傾向をつかもうとするだけで新しい発見があります。

AR(1)過程の生成  

1次の自己回帰モデルは  

$P_{t} =a + βP_{t-1} + σw_{t}$ 

で与えられる。このβを指定することで、AR(1)モデルがどのような性質をもっているのかを調べてみよう。

プログラムのコードはランダムウォークの時とほとんど変わらないです。

βを指定することによりαを決めなければならないが、これは期待値(μ)から求めることができる   AR(1)の期待値$E(P_{t})$は  

$E(P_{t})= E(α)+E(βP_{t-1}) +E(σw_{t})$

で与えられる。

$E = (P_{t}) = \mu $とすると、

$E(a) = a$、$E(w_{t})$であるから、

$$\mu = \frac{a}{1- \beta } $$

となる。価格の初期値を1と設定すると、  

α=(1-β)である。分散は  

$var(P_{t})=E(P_{t}-μ)^{2} = E(P_{t}^{2})-μ^{2}$

であり、α=0であれば

$$var(P_{t})= \frac{ \sigma ^{2}_{z}}{1- \beta ^{2}} $$

である。期待値も分散も時間とは独立している。

ここがランダムウォークとは違う点である。  

次にdef分を用いて、ar1というm個の価格から成るAR(1)の時系列を標準偏差をsigma、回帰係数をbeta、初期値をp0で指定してn個生成する関数を作ってみよう。

また、単位時間ステップ毎に価格の最大値と最小値も計算しておこう

In [7]:
def ar1(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)
    price=pd.Series(P)
    dprice=pd.Series(dP)
    return price,dprice,high,low

まず最初はbeta=0.9999、sigma=0.00632、n=10000、m=250、p0=1として時系列を生成します。

価格差の平均はゼロ、標準偏差は0.00632、歪度、尖度ともにゼロ近辺という結果になりました。  

β< 1である特徴が僅かながら現れているだろうか

In [8]:
price,dprice,high,low = ar1(0.9999,sigma,10000,250,1)
plt.figure(figsize=(8,4))
plt.subplot(1,2,1)
plt.plot(high,label="high",linestyle='--')
plt.plot(low,label="low",color='darkgray')
plt.title('max$P_t$-min$P_t$')
plt.xlabel('t')
plt.ylabel('$P_t$')
plt.legend(loc='upper left')
plt.subplot(1,2,2)
price.hist(bins=100,color='pink')
plt.xlabel('$P_{250}$')
plt.ylabel('frequency')
print("mean %2.5f std %2.5f skew %2.5f kurt %2.5f"\
%(dprice.mean(),dprice.std(),dprice.skew(),dprice.kurt()))
print("price std %2.5f"%price.std())
mean 0.00000 std 0.00633 skew -0.00001 kurt -0.00054
price std 0.09906

価格差の平均はゼロ、標準偏差は0.00632、歪度、尖度ともにゼロ近辺という結果になった。  

β< 1である特徴が僅かながら現れているだろうか。

結果はランダムウォークとほぼ同等である。区別はできない。次に結果を価格差の散布図にしてみます。1単位時間差を見ています。

In [9]:
plt.figure(figsize=(3,3))
mx=round(dprice.max(),2)
mn=round(dprice.min(),2)
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$')
plt.xticks([mn,0,mx])
Out[9]:
([<matplotlib.axis.XTick at 0x7f75140a7a58>,
  <matplotlib.axis.XTick at 0x7f7513c875c0>,
  <matplotlib.axis.XTick at 0x7f7513c82160>],
 <a list of 3 Text xticklabel objects>)

また、散布図による、価格差の特性は綺麗な円形を示している。  

次にAR(1)の特徴をさらに明確につかむために、βの値を動かしてみよう。

β= 0.99、0.5、1.003、1.005について実験してみよう。まず最初はβ=0.99から始めよう

In [10]:
price,dprice,high,low = ar1(0.99,sigma,10000,250,1)
plt.figure(figsize=(8,4))
plt.subplot(1,2,1)
plt.plot(high,label="high",linestyle='--')
plt.plot(low,label="low",color='darkgray')
plt.title('max$P_t$-min$P_t$')
plt.xlabel('t')
plt.ylabel('$P_t$')
plt.legend(loc='upper left')
plt.subplot(1,2,2)
price.hist(bins=100,color='pink')
plt.xlabel('$P_{250}$')
plt.ylabel('frequency')
print("mean %2.5f std %2.5f skew %2.5f kurt %2.5f"\
%(dprice.mean(),dprice.std(),dprice.skew(),dprice.kurt()))
print("price std %2.5f"%price.std())
mean -0.00000 std 0.00633 skew 0.00122 kurt -0.00668
price std 0.04473

結果は価格差の平均ゼロ、標準偏差000634、歪度と尖度はほぼゼロである。  

β= 0.99では、ある日数を経過する価格の変化率の最大値と最小値の幅が一定値に留まるようになる。

これが中心回帰、定常確率過程の特徴である。それでも定常性を得るために50日程度かかっているのが分かる。  

次にβ= 0.5にしてみよう

In [11]:
price,dprice,high,low = ar1(0.5,sigma,10000,250,1)
plt.figure(figsize=(8,4))
plt.subplot(1,2,1)
plt.plot(high,label="high",linestyle='--')
plt.plot(low,label="low",color='darkgray')
plt.title('max$P_t$-min$P_t$')
plt.xlabel('t')
plt.ylabel('$P_t$')
plt.legend(loc='upper left')
plt.subplot(1,2,2)
price.hist(bins=100,color='pink')
plt.xlabel('$P_{250}$')
plt.ylabel('frequency')
print("mean %2.5f std %2.5f skew %2.5f kurt %2.5f"\
%(dprice.mean(),dprice.std(),dprice.skew(),dprice.kurt()))
print("price std %2.5f"%price.std())
mean 0.00000 std 0.00730 skew 0.00070 kurt -0.00035
price std 0.00729

結果は価格差の平均ゼロ、標準偏差000730、歪度と尖度はほぼゼロである。  

β= 0.5ではほぼ瞬間的に中心回帰の定常確率過程の特徴が見えてくる。

次にAR(1) が発散するβが1を超える場合を見てみよう。最初はβ=1.003にしてみよう

In [12]:
price,dprice,high,low = ar1(1.003,sigma,10000,250,1)
plt.figure(figsize=(8,4))
plt.subplot(1,2,1)
plt.plot(high,label="high",linestyle='--')
plt.plot(low,label="low",color='darkgray')
plt.title('max$P_t$-min$P_t$')
plt.xlabel('t')
plt.ylabel('$P_t$')
plt.legend(loc='upper left')
plt.subplot(1,2,2)
price.hist(bins=100,color='pink')
plt.xlabel('$P_{250}$')
plt.ylabel('frequency')
print("mean %2.5f std %2.5f skew %2.5f kurt %2.5f"\
%(dprice.mean(),dprice.std(),dprice.skew(),dprice.kurt()))
print("price std %2.5f"%price.std())
mean 0.00000 std 0.00633 skew 0.00057 kurt 0.00336
price std 0.15265

結果は価格差の平均ゼロ、標準偏差0.00633、歪度と尖度はほぼゼロである。  β= 1.003ではまだ発散には至っていない。

むしろランダムウォークとの区別が難しい。β= 1.005ではどうだろうか

In [13]:
price,dprice,high,low = ar1(1.005,sigma,10000,250,1)
plt.figure(figsize=(8,4))
plt.subplot(1,2,1)
plt.plot(high,label="high",linestyle='--')
plt.plot(low,label="low",color='darkgray')
plt.title('max$P_t$-min$P_t$')
plt.xlabel('t')
plt.ylabel('$P_t$')
plt.legend(loc='upper left')
plt.subplot(1,2,2)
price.hist(bins=100,color='pink')
plt.xlabel('$P_{250}$')
plt.ylabel('frequency')
print("mean %2.5f std %2.5f skew %2.5f kurt %2.5f"\
%(dprice.mean(),dprice.std(),dprice.skew(),dprice.kurt()))
print("price std %2.5f"%price.std())
mean -0.00000 std 0.00636 skew -0.00215 kurt 0.00121
price std 0.20891

結果は価格差の平均ゼロ、標準偏差000636、歪度と尖度はほぼゼロである。β= 1.005では発散の一歩手前である。  

βの違いは、AR(1)の動きの特性に大きな影響を与えたであろうか? 

β> 1であってもその程度が小さければ、価格過程は破たんすることが無いことが分かる。

またそれは観測期間にも左右されることも分かる。

期間が短ければ、βが1よりも大きくても破たんせずにすむ確率は上がる。また、β< 1で1近辺の場合はどうであろうか? 

ランダムウォークとの区別はうまく行くだろうか? これもまた難しいという結果になった。  

自己回帰モデルの最大の特徴は価格の中心回帰の動きである。

ランダムウォークの価格の動きは時間の経過とともに分散が大きくなるのに対して、自己回帰モデルではある点に収束する。

そのために将来の予測を行うためには便利です。

ベルヌーイ試行  

金融資産の価格の動きの大きさがx に定められていると、上昇は+ x 、下落は- x となる。

このように結果が2つしかない実験をベルヌーイ試行とよび、とびとびな2つの値から成る独立な離散時間の確率変数列をベルヌーイ過程とよぶ。

このような2つのとびとびな確率変数からなる確率過程の性質を乱数生成器を用いて調べてみます。

ここではx を5とした時の1万個の確率変数を発生させその和を求め、それを1日分の取引としました。

In [14]:
def bernoulli(p,x,N,p0):
    s = (np.random.binomial(1, p, N)-0.5)*x
    P0=p0
    P=[]
    for i in range(N):
        P0=P0+s[i]
        P.append(P0)
    return P
P=bernoulli(0.5,5,10000,0)
plt.figure(figsize=(10,6))
plt.plot(P,color='darkgray')
plt.xlabel('steps $t$')
plt.ylabel('$\sum_1^t B_t \cdot 5$')
Out[14]:
<matplotlib.text.Text at 0x7f7517677240>

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

つぎにこの実験を10000回繰り返し、1回1回の実験で生成した確率変数の和をヒストグラムとして描いてみます。

In [15]:
N=10000
M=10000
Q=[]
for j in range(N):
    P=bernoulli(0.5,5,M,0)
    Q.append(P[M-1])
plt.figure(figsize=(8,5))
plt.figure.left=-0.1
plt.hist(Q,normed=True,histtype='stepfilled',color='pink',bins=25)
plt.xlabel('$\sum_1^t B_t \cdot 5$')
plt.ylabel('probability density function')
price=pd.Series(Q)
dprice=price.diff()
print("mean %2.5f std %2.5f skew %2.5f kurt %2.5f"\
%(dprice.mean(),dprice.std(),dprice.skew(),dprice.kurt()))
mean -0.01500 std 352.97407 skew -0.02056 kurt 0.03577

結果はかなり正規分布に近いベル型の形をしている。次に1日の取引の回数を3000回に減らして、その実験の回数は10000回とします。

In [16]:
N=10000
M=3000
Q=[]
for j in range(N):
    P=bernoulli(0.5,5,M,0)
    Q.append(P[M-1])
plt.figure(figsize=(8,5))
plt.figure.left=-0.1
plt.hist(Q,normed=True,histtype='stepfilled',color='pink',bins=25)
plt.xlabel('$\sum_1^t B_t \cdot 5$')
plt.ylabel('probability density function')
price=pd.Series(Q)
dprice=price.diff()
print("mean %2.5f std %2.5f skew %2.5f kurt %2.5f"\
%(dprice.mean(),dprice.std(),dprice.skew(),dprice.kurt()))
mean 0.00200 std 194.97902 skew 0.02862 kurt 0.06520

分布の形状は分布の幅が狭くなっている。これは1日の価格の動きが取引の回数に依存することを示しています。

2番目の特徴は、ベル型の形は維持しているが、分布の滑らかさは維持されていない。

これが最も重要な点で、これこそが本題です。  

ボラティリティは取引数と関係があるということです。

In [ ]: