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

# # モンテカルロで見える世界
# 
# インターネット上から金融関連の価格データを手に入れ、そのデータにトレンドが在るか無いかを見てきました。
# 
# トレンドがあれば、確定的トレンドと確率的トレンドに分類され、また、時間との関係では定常時系列と非定常時系列に分類されました。
# 
# 確定的トレンド、確率的トレンド、定常時系列、または時間トレンド、ランダムウォーク、AR（１）が時系列分析の３種の神器でした。
# 
# ## モンテカルロ・シミュレーションの利用 　
#  
#  日経平均株価の動きの３つのモデルは 　
#  
#   $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よりも１単位時間前の価格です。
#  
#  $w_{t}$は平均ゼロ、分散１の正規分布に従う確率変数であり、σは定数です。 
#  
# $B_{t}$は＋1,－1からなるベルヌーイ過程である。シミュレーションではσ＝0.00632とする。 　
#  
#  ## ランダムウォーク 　
#  
#  ランダムウォーク・モデルは 　
#  
#  $P_{t} = P_{t-1} + σw_{t}$ 
#  
#  で与えられる。
#  
#  ここではα= 0としてドリフト無しランダムウォークを検討する。
#  
#  σは一定である。 　
#  
#  まず１年間の$P_{t}$の動きを人工的に生成してみよう。
#  
#  np.randomというモジュールの中にある乱数生成器normalを用いて確率変数$w_{t}$を生成する。 
#  
#  株価の初期値 P を１とする。 
#  
#  株式市場の１年の営業日数を250日とし、250個の株価の動きを生成してみよう

# In[1]:


get_ipython().run_line_magic('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)


# これで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()))


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

# 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')


# 左側のチャートは１日目の価格、２日目の価格といった具合に250日間の価格の推移を描いていて、右側の頻度図は日々の価格の変化の大きさに対する頻度を表示している。 
# 
# 価格の差分は１回の試行の結果としては、かなりきれいな正規分布である。
# 
# #### プログラムコードの解説
# 
# 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+sigma*w)            価格の更新
# 
# 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回繰り返すことで、１日目、２日目の価格の最大値、最小値を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()))


# このモンテカルロシミュレーションでは250日間の価格をまず生成し、それを10000回繰り返すことで、時系列の特性を得ようとしました。 
# 
# 左の図は250日間の価格の推移の特性を調べています。
# 
# 10000回の試行の内の１日目の価格の最大値と最小値をプロットし、それを２日目、３日目にも行い、250日目まで繰り返すと価格の推移の期間構造が得られます。
# 
# 価格の最大値と最小値の幅がほぼ直線的（緩やかな曲線）に増加しているのが分かります。
# 
# 右図は250日目の価格の頻度を表している。分布はほぼベル型です。
# 
# 価格差の平均はゼロ、標準偏差は0.00633と設定値そのものである。歪度と尖度もほぼゼロです。 　
# 
# このモンテカルロシミュレーションから得られた時系列の特性は実際の価格の動きの特性を説明しているであろうか？　
# 
# 前回の結果を思い出してほしい

# Material 9 変化率の最大値、最小値の期間構造のグラフ参照
# 
# 最小値（low）のグラフは割とスムーズだが最大値（high）のグラフは凸凹である。
# 
# 一方、シミュレーションではほぼ曲線的に下落、上昇している。
# 
# パラメータが一定のシミュレーションでは現実の世界を説明できないが、そんなに外れてもいない。 　
# 
# 次に得られた時系列から１日間隔の変化率とそれを１日ずらした変化率の散布図を描いてみよう

# 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()


# まん丸の円が描かれた。実際のデータとはだいぶ違う。２つの変化率には相関がないことが見て取れます。
# 
# これが分散不均一性もボラティリティ・クラスタリングも無い世界である。 　
# 
# このように、実際の時系列を分析する際には、モデルの厳密さにこだわるのではなく、まずは大まかな傾向をつかもうとするだけで新しい発見があります。

# ## AR（1）過程の生成 　
# 
# １次の自己回帰モデルは 　
# 
#  $P_{t} =a + βP_{t-1} + σw_{t}$　 
# 
# で与えられる。このβを指定することで、AR（１）モデルがどのような性質をもっているのかを調べてみよう。
# 
# プログラムのコードはランダムウォークの時とほとんど変わらないです。
# 
# βを指定することによりαを決めなければならないが、これは期待値（μ）から求めることができる
# 　
#  AR（１）の期待値$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－β)である。分散は 　
#  
#  $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（１）の時系列を標準偏差を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())


# 価格差の平均はゼロ、標準偏差は0.00632、歪度、尖度ともにゼロ近辺という結果になった。 　
# 
# β< 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])


# また、散布図による、価格差の特性は綺麗な円形を示している。 　
# 
# 次にAR（１）の特徴をさらに明確につかむために、βの値を動かしてみよう。
# 
# β= 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())


# 結果は価格差の平均ゼロ、標準偏差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())


# 結果は価格差の平均ゼロ、標準偏差000730、歪度と尖度はほぼゼロである。 　
# 
# β= 0.5ではほぼ瞬間的に中心回帰の定常確率過程の特徴が見えてくる。
# 
# 次にAR（１） が発散するβが１を超える場合を見てみよう。最初はβ=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())


# 結果は価格差の平均ゼロ、標準偏差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())


# 結果は価格差の平均ゼロ、標準偏差000636、歪度と尖度はほぼゼロである。β= 1.005では発散の一歩手前である。 　
# 
# βの違いは、AR（１）の動きの特性に大きな影響を与えたであろうか？　
# 
# β> 1であってもその程度が小さければ、価格過程は破たんすることが無いことが分かる。
# 
# またそれは観測期間にも左右されることも分かる。
# 
# 期間が短ければ、βが１よりも大きくても破たんせずにすむ確率は上がる。また、β< 1で１近辺の場合はどうであろうか？　
# 
# ランダムウォークとの区別はうまく行くだろうか？　これもまた難しいという結果になった。 　
# 
# 自己回帰モデルの最大の特徴は価格の中心回帰の動きである。
# 
# ランダムウォークの価格の動きは時間の経過とともに分散が大きくなるのに対して、自己回帰モデルではある点に収束する。
# 
# そのために将来の予測を行うためには便利です。 

# ## ベルヌーイ試行 　
#  
#  金融資産の価格の動きの大きさがx に定められていると、上昇は＋ x 、下落は－ x となる。
#  
#  このように結果が２つしかない実験をベルヌーイ試行とよび、とびとびな２つの値から成る独立な離散時間の確率変数列をベルヌーイ過程とよぶ。
#  
#  このような２つのとびとびな確率変数からなる確率過程の性質を乱数生成器を用いて調べてみます。
#  
#  ここではx を５とした時の１万個の確率変数を発生させその和を求め、それを１日分の取引としました。

# 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$')


# 実際の価格の動きに似た動きが再現できていると思う。これはランダムウォークであり、一定の方向に価格が上昇したり下落しているのは確率的トレンドである。 　
# 
# 
# つぎにこの実験を10000回繰り返し、１回１回の実験で生成した確率変数の和をヒストグラムとして描いてみます。

# 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()))


# 結果はかなり正規分布に近いベル型の形をしている。次に１日の取引の回数を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()))


# 分布の形状は分布の幅が狭くなっている。これは１日の価格の動きが取引の回数に依存することを示しています。
# 
# ２番目の特徴は、ベル型の形は維持しているが、分布の滑らかさは維持されていない。
# 
# これが最も重要な点で、これこそが本題です。 　
# 
#  ボラティリティは取引数と関係があるということです。

# In[ ]: