トレンドをモデル化しよう

トレンドのモデル化を考えてみよう。 次の一次式の形  

$$Y_{t}=α+β t+u_{t} (1) $$

で表されるトレンドを確定的トレンド、または時間トレンドとよぶ。

時間の経過にともない、その経過時間に比例して価格が変化するモデルである。

αとβは回帰係数である。 ここで t は時間を表し、$u_{t}$は定常なかく乱項または残差項です

線形回帰モデルに息を吹き込む

線形回帰モデル  

2変数からなるモデルを考えよう。

それを  

$$E(Y|X_{i}) =β_{1}+β_{2} X_{i} (2) $$

と表現する。

Yを被説明変数、Xを説明変数という。

$β_{1}$、$β_{2}$は定数です。 これらの定数は説明変数と被説明変数の関係をとらえていて、回帰係数とよばれる。

(2)式は係数と変数に関して線形であるという。それは係数と変数のそれぞれが1次の関数であるからです。

$E()$は条件付き期待値を表している。$E(Y|X_{i})$ は $X_{i}$ の時のY の期待値を表している。

(2)式は $X_{i}$ の時の Y の期待値が必ず $β_{1}+β_{2} X_{i}$ である確定的な関係を表現している。 これを母集団回帰関数とよぶ。 実際には平均からの乖離があるのでそれを表現すると(2)式は  

$$Y_{i} = β_{1}+β_{2} X_{i}+u_{i} = E(Y|X_{i})+u_{i} (3) $$

となる。$E(Y|X_{i})$ は説明変数における被説明変数の平均的な値であり確定的な成分です。 一方、$u_{i}$ は予測不可能であるので確率的な成分です。(3)式を母集団回帰式とよび の数だけ存在する。  

標本と母数  

すべてのデータを母集団というが、 (1)、(2)、(3)式での説明はこの母集団が手に入るという前提に立ち、実際に直面する問題には触れずに来た。

しかし、金融の専門家が母集団を手に入れることは不可能であり、 手にできるのは母集団の部分です。

それを標本という。 母集団の特性を表現する定数(平均、分散など)を母数(parameter)といい、(1)、(2)、(3)式では $β_{1}$ と $β_{2}$ が母数です。

母数は不変だが、未知です。  

標本の特性を表現する定数を統計量(statistics)という。 我々が観測できるのは母集団から抽出された標本であるので、 母数を手に入れることはできない。

しかし、標本から統計量を算出できるので、標本から母数を推定できる

標本回帰式

我々が手にできるのはある固定値 X に対する限られた Y の値です。

従って、観察された実現値、または得られた標本から母集団回帰関数と同等な標本回帰関数を推定しなければならない。

それを

$$\hat{Y}_{i} = \hat{ \beta } _{1} + \hat{\beta } _{2}X_{i}$$

で表す。ここで $\hat{Y}_{i}$、$\hat{ \beta } _{1}$、$\hat{\beta } _{2}X_{i}$ は統計量である。また上式は

$$\hat{Y}_{i} = \hat{ \beta } _{1} + \hat{\beta } _{2}X_{i}+\hat{u}_{i} (3) $$

と書くこともできる。

$\hat{u}_{i}$ は統計量であり、残差項とよぶ。

この式を標本回帰式とよぶ。

$\hat{ \beta } _{1}$、$ \hat{\beta } _{2}$、$\hat{u}_{i}$ は統計量であり、得られた標本からある規則や計算方法により求められる。

この規則や計算方法を推定量(estimator)とよび、そこから得られる数値を推定量(estimate)という。

$\hat{ \beta } _{1}$、$ \hat{\beta } _{2}$、$\hat{u}_{i}$ は推定量です。

最小二乗法

線形回帰関数の統計量は最小二乗法によって求められる。

(3)式を書き直し

$$ \hat{u}_{i} = \hat{Y}_{i} - (\hat{ \beta } _{1} + \hat{\beta } _{2}X_{i})$$

を得る。この両辺を2乗して、残差の2乗和を最小にするような $\hat{ \beta } _{1}$ と$\hat{\beta } _{2}$ を求める。

そうすると

$$\hat{ \beta } _{1} = \frac{ \Sigma X^{2}_{i} \Sigma Y_{i}- \Sigma X_{i} \Sigma X_{i} Y_{i}}{N \Sigma X_{i}^{2}-( \Sigma X_{i})^{2}} $$$$\hat{ \beta } _{2} = \frac{ \Sigma X_{i} Y_{i} - \Sigma X_{i} \Sigma X_{i} Y_{i}}{N \Sigma X^{2}-( \Sigma X_{i})^{2}}$$

が得られる。

ここで N は標本数です。

最小二乗法の仮定

統計モデルは幾つかの仮定のもとに成り立つ。

古典的線形回帰モデルの仮定をここで列挙しておきます。  

*回帰関数は線形でなければならない。  

* $X_{i}$ は確率変数であってはならない。  

* $u_{i}$ の平均はゼロである。  

* $u_{i}$ の分散は一定である。  

* $u_{i}$ と $u_{i+j}$ の相関はゼロである。j≠1

* $u_{i}$​ と $X_{i}$​ の共分散はゼロである。

ここで注目してほしいのは残差項に関する仮定が多いことです。  

推定の信頼性  

単純な時間トレンドの標本線形回帰式には $\hat{ a }$ , $\hat{ \beta }$ , $\hat{ u }$ の3つの推定値があります。  

標準誤差

標準誤差(Standard Error,Std Err)は母数の推定値(統計量)と未知の母数との差です。

これは統計量の正確さの測度であり

$$se(\hat{ \beta } _{1}) = \sqrt{\frac{ \Sigma X_{i}^{2} }{N \Sigma( X_{i}- \overline{X})^{2}} \sigma} $$$$se(\hat{ \beta } _{2}) = \frac{\sigma}{\sqrt{\Sigma( X_{i}- \overline{X})^{2}}} $$

から求めることができる。

$\overline{X}$ の標本平均です。

ここで σ は次式で与えられる $u_{i}$ の分散で推定の標準誤差とよばれる。

$$\sigma^{2} = E(\hat{{\sigma}}^{2}) = \frac{E(\Sigma u_{i}^{2})}{N-2} $$

この母数は未知であるので、その推定量は、

$$\hat{{\sigma}}^{2} = \Sigma\frac{\hat{u}_{i}^{2}}{N-2} = \frac{\Sigma[Y_{i}-E(Y|X_{i})]^{2}}{{N-2}} $$

である。標準誤差は、統計量のバラツキ具合、つまり精度の測度であり、

推定の標準誤差は標本線形回帰線とデータとの適合度(goodness of fit)の目安です。

※決定係数 $R^{2}$ (R-squared)

$R^{2} $ (R-squared) は標本線形回帰線とデータとの間の適合度(goodness of fit)を表す測度として知られ、 決定係数(coefficient of determination)という。

次式で与えられる。

$$r^{2} = \frac{\Sigma(\hat{Y}_{i}-\overline{Y})^{2}}{\Sigma({Y}_{i}-\overline{Y})^{2}} $$

ここで$\overline{Y}$はYの標本平均です。

決定係数が1に近いほど、相対的なバラツキは少ない。

自動調整済み決定係数(Adj R-squared)は説明変数の数の効果を考慮した係数です。

説明変数の数が多くなると決定係数は良くなる傾向があるので、その分を調整している。  

2乗平均平方根誤差

実測値と予測値の間の残差の目安が2乗平均平方根誤差 (RMSE:root-mean-square error)です。

$$RMSE^{2} = MSE = E(\hat{Y}-Y)^{2} = var(\hat{Y})+E[E(\hat{Y})-Y]^{2} $$

で与えられる。ここで$E(\hat{Y}-Y)^{2}$はバイアスです。

$var(\hat{Y})$は$\hat{Y}$のの分散であって、実測値と予測値の間の残差の2乗の平均値を平均2乗誤差といい、2乗平均平方根誤差はその平方根です。

小さいほうが良い。$var(\hat{Y})$がゼロであれば、

$RMSE^{2} = E(\hat{Y}-Y)^{2}$です。

※信頼区間

母平均μが確率1-xで区間 $μ^{*}$と$μ^{**}$の間にいるとき、1-x を信頼係数(the confidence coefficient)、または信頼水準(the significance level)といい、$μ^{**}$<μ< $μ^{*}$を信頼区間という。

$μ^{*}$、$μ^{**}$をそれぞれ信頼上限、信頼下限という。xを有意水準(the level of significance)という。  

※p-値

信頼区間を標準化変換を用いて、

$$ \theta _{0}- \delta < \hat{ \theta } <\theta _{0}+\delta$$

と書き換えることができる。

$\theta _{0}$は帰無仮説の母数の値です。

$\delta$は平均がゼロ、分散が1の標準正規分布に従う

$$Z = ( \widehat{ \theta } - \theta _{0})/se(\widehat{ \theta })$$

から得ることができる。

母標準誤差が未知であるとき、 母標準誤差を$se(\widehat{ \theta })$に置き換え、ステューデントのt分布推定量(回帰係数)の信頼区間が得られる。

この関係を利用すると推定量$\widehat{ \theta}$(回帰係数)の信頼区間が得られる。

$\theta _{0}$をゼロとし、標本平均を$\theta _{0}$とすると確率1-p の信頼区間は

$$0 - t_{p,n}se( \widehat{ \theta }) < \widehat{ \theta } < 0 + t_{p,n}se( \widehat{ \theta })$$

で与えられる。

ここで$t_{p,n}$は確率1-p の信頼区間を与える自由度n-1のt分布の値で臨界値とよぶ。

ここで与えられる区間を採択域(the region of acceptance)といい、この外側の領域を棄却域(the region(s)of rejection)という。

これを危険域(the critical region(s))と呼ぶことがある。

そうすると大きな|t|値は帰無仮説の棄却域にいることになる。

$$t_{p,n} = \frac{ \widehat{ \theta } }{se( \widehat{ \theta })} $$

であるから、$t_{p,n}$からp-値(p-value)を計算できる。

p-値臨界値が与える棄却域の確率であると定義される。

一般に次の表が棄却、棄却が難しいの目安になります。

関係(pはp-値を表す) 解釈
p<0.01 帰無仮説を棄却する
0.01<p<0.1 帰無仮説を棄却するに足る
0.1<p 帰無仮説を棄却するのは難しい

「帰無仮説(否定)を棄却する」とは0.01以下の確率でしか起こらないことが起こった、ということです。

「帰無仮説(否定)の棄却は難しい」は棄却するに十分な証拠がないということです。

統計学の目的は極力誤った判断を減らすことにあります。

                                                                                                                                                                                                   ※:分析の判断目安

日経平均株価の確定的トレンド

時間の経過とともに価格が線形に上昇、または下落する傾向をもつ時、その時系列は次のようにモデル化される。  

$$Y_{t}=α+β・t+u_{t} $$

ここでYtは日経平均株価、t は経過時間である。$Y_{t}$はトレンドの傾きである。$u_{t}$は残差項である。

このようなトレンドを確定的トレンドという。$Y_{t}$が一時的にトレンドから大きく乖離したとしても、その長期的なトレンドに変化はなく、

日経平均株価には確定的トレンドに復帰する力があると考えます。

このような確定的トレンドが成り立つ理由として、 株式市場の動向が景気循環と関係があり、景気の悪化により株価が下落しても、

その景気の悪化が景気循環の一次的な現象であり、 時間の経過と共に回復し、それにともない株式市場も回復すると考えるからです。  

静的分析  

日経平均株価はバブル崩壊までは強い上昇トレンドを経験し、

その後は幾つかの上昇トレンドと下落トレンドを繰り返しながら、いまだに最高値を更新できないでいます。

今回の分析は対数価格を用いる。

興味のある場合は価格で分析を行ってほしいと思います。

その違いが判るはずです。

In [70]:
import pandas as pd
import pandas_datareader.data as web
import statsmodels.api as sm
import numpy as np
end='2016/9/30'
n225 = web.DataReader("NIKKEI225", 'fred',"1949/5/16",end).dropna()
lnn225=np.log(n225.dropna())
lnn225.columns=['Close']
y=lnn225
x=range(len(lnn225))
x=sm.add_constant(x)
model=sm.OLS(y,x)
results=model.fit()

print(results.summary())
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                  Close   R-squared:                       0.756
Model:                            OLS   Adj. R-squared:                  0.756
Method:                 Least Squares   F-statistic:                 5.196e+04
Date:                Wed, 02 Aug 2017   Prob (F-statistic):               0.00
Time:                        12:47:55   Log-Likelihood:                -18601.
No. Observations:               16769   AIC:                         3.721e+04
Df Residuals:                   16767   BIC:                         3.722e+04
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          6.1967      0.011    546.859      0.000       6.174       6.219
x1             0.0003   1.17e-06    227.943      0.000       0.000       0.000
==============================================================================
Omnibus:                      560.110   Durbin-Watson:                   0.000
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              284.479
Skew:                          -0.119   Prob(JB):                     1.68e-62
Kurtosis:                       2.408   Cond. No.                     1.94e+04
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.94e+04. This might indicate that there are
strong multicollinearity or other numerical problems.

結果の解釈  

標本数(Number of Ovservations)は16769である。  

決定係数(R-sqaured)は0.756である。

これはモデルがどれ程うまく実現値を説明しているかの尺度であり、最大値は1である

回帰係数αとβの帰無仮説はそれぞれα=0とβ=0である。

それぞれの標準誤差(Standard Error,Std Err)は0.011と0.0000である。αの帰無仮説が棄却されないと$Y_{t}=β・t+u_{t} $である。

αのp-値がゼロであることで帰無仮説は棄却されα= 6.1967である。

βの帰無仮説が棄却されないとYt=α+utでである。βのp-値がゼロであることで帰無仮説は棄却されβ= 0.0003である。

従って1949年以降の日経平均株価は上昇トレンドをもっている。  

$$lnY_{t}=6.1967+0.0003・t$$

とモデル化された。

[95.0% Conf.int.]は信頼区間である 以上のレポートだけでは結果の判断はできない。

最小二乗法の仮定を思い出してほしい。まず、モデルの期待値と実際の日経平均株価をプロットして、その差を目視で確かめよう

Statsmobelの解説 : statmodelの線形回帰分析メゾット

sm.OLS(y,x)

xは説明変数、yは被説明変数である。ここではyを日経平均株価、xは経過時間とする。

プログラムコードの解説

n225 = web.DataReader("NIKKEI225", 'fred',"1949/5/16",end).dropna()

fredからのデータの取得

lnn225=np.log(n225.dropna())

対数価格の計算、n/aとなるデータdropna()で削除

lnn225.columns=['Close']

数値の列にCloseと名前をつけている。

列の名前をCloseに設定

y=lnn225

変数yの作成

x=range(len(lnn225))

x=sm.add_constant(x)

変数xの作成

lnn225と同じ長さに設定。Rangeによって1,2,...lnn225のデータ数までの整数を設定

model=sm.OLS(y,x)

results=model.fit()

線形回帰の切片のために列(要素)を作成

線形回帰分析を設定をmodelとして保存

modelをfitを用いて最適化

In [39]:
%matplotlib inline
import matplotlib.pyplot as plt
plt.plot(y,label='Close',color="darkgray")
results.fittedvalues.plot(label='prediction',style='--')
plt.ylabel('log(n225 index)')
plt.legend(loc='upper left')
Out[39]:
<matplotlib.legend.Legend at 0x7f935c7dbac8>

fittedvaluesはt時の予測値、または期待値を与えてくれる。 次に残差項だけを取り出し、時系列としてチャートを描いてみよう。residは残差項を与えてくれる。

In [40]:
results.resid.plot(color="seagreen")
plt.ylabel('residual')
Out[40]:
<matplotlib.text.Text at 0x7f935bef25c0>

明らかに残差の時系列の中にトレンドが在りそうだ。残差のヒストグラムを描いてみるとどうだろうか?

In [41]:
results.resid.hist(bins=100,color="lightgray")
plt.xlabel('residual')
plt.ylabel('frequency')
Out[41]:
<matplotlib.text.Text at 0x7f935bf335f8>

残差はいくつかの分布で構成されていて、残差の平均と分散はいくつか存在するように見える。残差が古典的線形回帰モデルの条件を満たしているようには見えない。

結果の解釈

時間トレンドが算出する日経平均株価の期待値は、実際の終値の推移を直線でうまく説明しているように思える。

しかし残差項のチャートではさらに期間の短いトレンドが存在するのではないかと疑わざるを得ない。

最初の乖離がマイナス方向に大きく、時間の経過と共にそれが小さくなり、次に残差はプラスに転換し、その後は時間の経過と共にプラス方法に乖離幅が大きくなりバブルのピークの時に最大となっている。

そして、その後再びゼロに近づき、マイナス方向に拡大していくことが分かる。このようなトレンドの存在は残差のヒストグラムに幾つかの分布が混合しているのではないかと思われることからも確認できる。残差の平均と分散が一定しているとは思えない。

これらの結果は目による確認なので残差に関してはさらなる統計的な分析が必要である.

景気循環と日経平均株価の確定的トレンド

1949年以降直近までの日経平均株価の時間トレンドの分析では明確な答えは得られていない。

分析結果から直感的にはそれぞれの景気循環期に対して時間トレンドを求めることでより明確な判断ができるのではないかと考えた。

そこで内閣府の定めた景気循環期に順じて分析してみよう。また、それに加えてバブル経済期の始まりから日経平均株価のピークまでとピークから暴落終焉までの期間を設定した

戦後復興期(recover)

戦後復興期は2つの景気循環期で構成した。通称、特需景気と投資景気からなる。朝鮮戦争の特需から始まり、その休戦協定により終結した。

この休戦協定は1951年7月10日から始まり1953年7月27日に締結されている。 R-squared=0.76と1949年以降のデータをすべて用いた場合よりも改善している。

p-値から回帰係数と切片の帰無仮説を棄却することはできないので、この2つの推定値には意味がある。

それぞれの景気循環期と時間トレンドの間には何か関係を見出せるかもしれないという期待を抱かせる。

モデルが算出する日経平均株価の期待値と終値を比べて時間トレンドで指数の動きが説明できるかどうかを見てみよう。

In [44]:
y=lnn225.ix[:'1954/11/30'].dropna()
x=range(len(y))
x=sm.add_constant(x)
model=sm.OLS(y,x)
results=model.fit()
print(results.summary())
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                  Close   R-squared:                       0.762
Model:                            OLS   Adj. R-squared:                  0.761
Method:                 Least Squares   F-statistic:                     4438.
Date:                Wed, 02 Aug 2017   Prob (F-statistic):               0.00
Time:                        12:15:38   Log-Likelihood:                -65.389
No. Observations:                1391   AIC:                             134.8
Df Residuals:                    1389   BIC:                             145.3
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          4.5236      0.014    332.550      0.000       4.497       4.550
x1             0.0011   1.69e-05     66.618      0.000       0.001       0.001
==============================================================================
Omnibus:                      160.558   Durbin-Watson:                   0.003
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              121.287
Skew:                           0.624   Prob(JB):                     4.60e-27
Kurtosis:                       2.268   Cond. No.                     1.60e+03
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.6e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
In [45]:
results.resid.std()
Out[45]:
0.2537083247963737
In [46]:
plt.plot(y,label='Close',color='darkgray')
results.fittedvalues.plot(label='prediction',style='--')
plt.legend(loc='upper left')
plt.ylabel('log(n225 index)')
Out[46]:
<matplotlib.text.Text at 0x7f935bdce320>

この時期の日経平均株価は必ずしも景気の山と谷に一致していない。また、景気循環期をさらに別の基準で分割する必要がありそうです。

高度経済成長期 (growth)

神武景気から始まりニクソン不況で終わる4つの景気循環期を高度経済成長期は含んでいる。

神武景気から岩戸景気の山までは日経平均株価は強いが、その後の証券不況で長い低迷期に入る。

オリンピックの終焉と共に不況に入った日本経済は低迷し、山一證券の取り付け騒ぎに至った。

大蔵大臣田中角栄による無制限、無担保の日銀特融により、事態は沈静化したが、日本経済が危機的状況に陥れば政府がラストリゾートになるという考え方が生まれた。

いざなぎ景気、ニクソンショックの期間において日経平均株価は堅調に推移した

In [47]:
y=lnn225.ix['1954/12/1':'1971/12/31'].dropna()
x=range(len(y))
x=sm.add_constant(x)
model=sm.OLS(y,x)
results=model.fit()
print(results.summary())
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                  Close   R-squared:                       0.824
Model:                            OLS   Adj. R-squared:                  0.824
Method:                 Least Squares   F-statistic:                 1.995e+04
Date:                Wed, 02 Aug 2017   Prob (F-statistic):               0.00
Time:                        12:15:53   Log-Likelihood:                 257.20
No. Observations:                4272   AIC:                            -510.4
Df Residuals:                    4270   BIC:                            -497.7
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          6.1675      0.007    884.618      0.000       6.154       6.181
x1             0.0004   2.83e-06    141.244      0.000       0.000       0.000
==============================================================================
Omnibus:                      396.700   Durbin-Watson:                   0.002
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              516.452
Skew:                           0.851   Prob(JB):                    7.14e-113
Kurtosis:                       3.082   Cond. No.                     4.93e+03
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 4.93e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
In [48]:
plt.plot(y,label='Close',color='seagreen')
results.fittedvalues.plot(label='prediction',style='--')
plt.legend(loc='upper left')
plt.ylabel('log(n225 index)')
Out[48]:
<matplotlib.text.Text at 0x7f935ba5c0f0>

証券不況の間しばらく日経平均株価が低迷したことが分かる。証券不況まで、証券不況の間、そして証券不況からの回復とさらに期間を分けることで時間トレンドで日経平均株価の動きを説明できそうです。

安定期 (stable)

列島改造景気から2度目の円高不況までの4つの経済循環期が日経平均株価が安定上昇した時期です。

この時期には第一次石油危機、第2次石油危機、アメリカとの貿易摩擦、プラザ合意などが含まれている。

町工場の倒産も続出した。鉄鋼、造船、石油産業は構造不況業種とよばれていた。

それでも第3次産業の情報処理産業が日本経済を牽引し始めた時期ですね。

In [49]:
y=lnn225.ix['1972/1/1':'1986/11/30'].dropna()
x=range(len(y))
x=sm.add_constant(x)
model=sm.OLS(y,x)
results=model.fit()
print(results.summary())
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                  Close   R-squared:                       0.910
Model:                            OLS   Adj. R-squared:                  0.910
Method:                 Least Squares   F-statistic:                 3.791e+04
Date:                Wed, 02 Aug 2017   Prob (F-statistic):               0.00
Time:                        12:15:57   Log-Likelihood:                 2451.9
No. Observations:                3768   AIC:                            -4900.
Df Residuals:                    3766   BIC:                            -4887.
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          8.1001      0.004   1969.386      0.000       8.092       8.108
x1             0.0004   1.89e-06    194.702      0.000       0.000       0.000
==============================================================================
Omnibus:                      548.042   Durbin-Watson:                   0.004
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              813.397
Skew:                           1.077   Prob(JB):                    2.36e-177
Kurtosis:                       3.736   Cond. No.                     4.35e+03
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 4.35e+03. This might indicate that there are
strong multicollinearity or other numerical problems.

日本経済にとって大変に難しい時期であるにも関わらず日経平均株価は順調に推移した。R-squared=0.91と前の2つの景気循環期よりも結果は良好である。チャートを用いて確認してみよう

In [50]:
plt.plot(y,label='Close',color='hotpink')
results.fittedvalues.plot(label='prediction',style='--')
plt.legend(loc='upper left')
plt.ylabel('log(n225 index)')
Out[50]:
<matplotlib.text.Text at 0x7f935b972320>

非常に安定した上昇トレンドが見て取れる。日本経済の高度成長時代が終わり安定期に入ると思われた。

人々の心は熱狂的なブームから一気にパニックに陥り、不安に駆られた時期です。

設備投資と個人消費が見込めず、急激な円高で貿易の拡大が見込めない中、財政出動により構造変化を成し遂げた。証券市場は経済の構造的変化に対応してその姿を変え、マネーフローに大きな変化があらわれました。

貸付市場から証券市場への相対的な移行でもある時期でした。

このような金融の構造変化の中で不況下の株高がもたらされた

バブル成長期 (bubble)

プラザ合意をきっかけにした円急騰が交易条件の改善をもたらし、それが実質所得の増加と結びつき内需主導の景気拡大が自律的に発生した。

上昇を続けていた株式市場はブラックマンデーにより大幅下落するが、そのダメージからいち早く回復したため、バブルであるという認識が強くなった。

日本発の世界株式の暴落を避けたい大蔵省は決算弾力化方針を打ち出しました。

また、銀行の貸し出し過多、証券不祥事に代表される損失補てんなどのモラルを逸脱した金融機関の業務拡大もバブルの要因でした。

また、大企業のエクイティーファイナンスと外債発行による新しい資金調達手段もその要因の1つに挙げられる

In [51]:
y=lnn225.ix['1986/12/1':'1993/10/31'].dropna()
x=range(len(y))
x=sm.add_constant(x)
model=sm.OLS(y,x)
results=model.fit()
print(results.summary())
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                  Close   R-squared:                       0.215
Model:                            OLS   Adj. R-squared:                  0.215
Method:                 Least Squares   F-statistic:                     467.5
Date:                Wed, 02 Aug 2017   Prob (F-statistic):           7.99e-92
Time:                        12:16:02   Log-Likelihood:                 305.16
No. Observations:                1707   AIC:                            -606.3
Df Residuals:                    1705   BIC:                            -595.4
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         10.2897      0.010   1050.279      0.000      10.271      10.309
x1            -0.0002   9.95e-06    -21.622      0.000      -0.000      -0.000
==============================================================================
Omnibus:                       40.208   Durbin-Watson:                   0.005
Prob(Omnibus):                  0.000   Jarque-Bera (JB):               22.069
Skew:                           0.077   Prob(JB):                     1.61e-05
Kurtosis:                       2.465   Cond. No.                     1.97e+03
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.97e+03. This might indicate that there are
strong multicollinearity or other numerical problems.

バブルによる上昇相場とバブル終焉による暴落相場を含んでいるためにR-squared=0.22と好ましい数値ではない。チャートにより確認してみよう

In [52]:
plt.plot(y,label='Close',color='seagreen')
results.fittedvalues.plot(label='prediction',style='--')
plt.legend()
plt.ylabel('log(n225 index)')
Out[52]:
<matplotlib.text.Text at 0x7f935b804b38>

バブルのピークの前と後では明らかにトレンドの形成の仕方が異なる。バブルの前では明確な強い上昇トレンド、ピークの後は明確な下落トレンドを形成している

バブル成長期(bubble)(日経平均株価ピークまで)  

特金・ファントラが貯蓄超過経済を象徴した。

また、画一的な機関投資家の運用姿勢も株価の上下動を大きくしたと考えられた。

NTT株売り出しによる株フィーバーは株式ブームの象徴であった。

ブラックマンデー以降も上昇を続ける株価はPER70倍前後を維持した。

In [53]:
y=lnn225.ix['1986/12/1':'1989/12/31'].dropna()
x=range(len(y))
x=sm.add_constant(x)
model=sm.OLS(y,x)
results=model.fit()
print(results.summary())
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                  Close   R-squared:                       0.913
Model:                            OLS   Adj. R-squared:                  0.913
Method:                 Least Squares   F-statistic:                     8011.
Date:                Wed, 02 Aug 2017   Prob (F-statistic):               0.00
Time:                        12:16:06   Log-Likelihood:                 1142.0
No. Observations:                 765   AIC:                            -2280.
Df Residuals:                     763   BIC:                            -2271.
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          9.9122      0.004   2520.009      0.000       9.904       9.920
x1             0.0008   8.91e-06     89.506      0.000       0.001       0.001
==============================================================================
Omnibus:                        1.148   Durbin-Watson:                   0.043
Prob(Omnibus):                  0.563   Jarque-Bera (JB):                1.193
Skew:                           0.092   Prob(JB):                        0.551
Kurtosis:                       2.943   Cond. No.                         882.
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

また回帰係数は今までで最も大きな値である。これだけでは強いトレンドが在るとは判断できない。そこで目視でモデルの期待値と実際の日経平均株価の推移を比べてみよう

In [54]:
print("return ",np.exp(y.Close).pct_change().mean()*250)
print("volatility ",y.Close.diff().std()*np.sqrt(250))
print("std of residual",results.resid.std())
plt.plot(y,label='Close',color='darkgray')
results.fittedvalues.plot(label='prediction',style='--')
plt.legend(loc='upper left')
plt.ylabel('log(n225 index)')
return  0.2624257002375478
volatility  0.177502142933
std of residual 0.054413688241352745
Out[54]:
<matplotlib.text.Text at 0x7f935bba9438>

米国での1987年11月のブラックマンデーの影響による暴落で上昇トレンドは一旦途切れたように見えるが順調に回復している。さらに残差を調べてみよう

In [55]:
results.resid.hist(bins=10,color="mistyrose")
plt.xlabel('residual')
plt.ylabel('frequency')
Out[55]:
<matplotlib.text.Text at 0x7f935b622e10>

ブラックマンデー後の大きな上昇トレンドからの乖離を含んでいるにも関わらす、目視では残差は1つの分布にしたがい平均、分散は一定しているように見える

バブル暴落時(日経平均株価ピークから) (bubble)

1990年の年初から1992年の夏まで若干の上下動はあるものの急激な下落である。この相場をバブルの上昇期のように時間トレンドでとらえることができるだろうか?

In [56]:
y=lnn225.ix['1990/1/1':'1992/8/31'].dropna()
x=range(len(y))
x=sm.add_constant(x)
model=sm.OLS(y,x)
results=model.fit()
print(results.summary())
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                  Close   R-squared:                       0.816
Model:                            OLS   Adj. R-squared:                  0.815
Method:                 Least Squares   F-statistic:                     2891.
Date:                Wed, 02 Aug 2017   Prob (F-statistic):          3.24e-242
Time:                        12:16:16   Log-Likelihood:                 622.33
No. Observations:                 656   AIC:                            -1241.
Df Residuals:                     654   BIC:                            -1232.
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         10.4342      0.007   1425.476      0.000      10.420      10.449
x1            -0.0010   1.93e-05    -53.767      0.000      -0.001      -0.001
==============================================================================
Omnibus:                       49.072   Durbin-Watson:                   0.038
Prob(Omnibus):                  0.000   Jarque-Bera (JB):               39.778
Skew:                          -0.516   Prob(JB):                     2.30e-09
Kurtosis:                       2.375   Cond. No.                         756.
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

結果はR-squaredも0.816と決して悪くはない。期待値と実測値をチャートに描いてみよう

In [57]:
print("return ",np.exp(y.Close).pct_change().mean()*250)
print("volatility ",y.Close.diff().std()*np.sqrt(250))
print("std of residual",results.resid.std())
plt.plot(y,label='Close',color='seagreen')
results.fittedvalues.plot(label='prediction',style='--')
plt.legend()
plt.ylabel('log(n225 index)')
return  -0.2494591088293264
volatility  0.287514843973
std of residual 0.0937749947107439
Out[57]:
<matplotlib.text.Text at 0x7f935b4fb780>

こちらも見た目うまく回帰できているように見える。次に残差のヒストグラムを見てみよう

In [58]:
results.resid.hist(bins=10,color="mistyrose")
plt.xlabel('residual')
plt.ylabel('frequency')
Out[58]:
<matplotlib.text.Text at 0x7f935b445be0>

この結果からは、直感的に残差が1つの分布にしたがうのではなく、よって平均がゼロで分散が一定しているというふうには感じられない

経済変革期(reform) 

カンフル景気から欧州危機までとそれ以降を経済変革期としている。

バブル崩壊後も大ブームの酔いから醒めることはできずに、政府が景気後退を認めるのは1992年1月であった。

ヨーロッパのEMS(欧州通貨システム)の混乱を狙った通貨投機の後、1993年から1995年まで急激な円高が進み、日本経済は再び混乱した。

経済変革期には4つの気循環期が含まれている。また、この時期は、拡張期よりも後退期に強い印象がある。

平成第2、3次不況、世界同時不況、そして欧州危機と続く。

In [59]:
y=lnn225.ix['1993/11/1':].dropna()
x=range(len(y))
x=sm.add_constant(x)
model=sm.OLS(y,x)
results=model.fit()
print(results.summary())
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                  Close   R-squared:                       0.138
Model:                            OLS   Adj. R-squared:                  0.138
Method:                 Least Squares   F-statistic:                     902.6
Date:                Wed, 02 Aug 2017   Prob (F-statistic):          4.76e-184
Time:                        12:16:26   Log-Likelihood:                -528.67
No. Observations:                5631   AIC:                             1061.
Df Residuals:                    5629   BIC:                             1075.
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          9.7139      0.007   1371.196      0.000       9.700       9.728
x1         -6.547e-05   2.18e-06    -30.043      0.000   -6.97e-05   -6.12e-05
==============================================================================
Omnibus:                      865.725   Durbin-Watson:                   0.003
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              213.450
Skew:                          -0.124   Prob(JB):                     4.47e-47
Kurtosis:                       2.079   Cond. No.                     6.50e+03
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 6.5e+03. This might indicate that there are
strong multicollinearity or other numerical problems.

結果はいままでで最悪である。R-squared=0.138と共に最悪の結果である。チャートを用いて確認してみよう

In [60]:
print("return ",np.exp(y.Close).pct_change().mean()*250)
print("volatility ",y.Close.diff().std()*np.sqrt(250))
print("std of residual",results.resid.std())
plt.plot(y,label='Close',color='hotpink')
results.fittedvalues.plot(label='prediction',style='--')
plt.legend()
plt.ylabel('log(n225 index)')
return  0.021601615793050276
volatility  0.241061913359
std of residual 0.26581255278187066
Out[60]:
<matplotlib.text.Text at 0x7f935b3968d0>

チャートからも回帰の期間の取り方が適切でない様子が見て取れる。残差のヒストグラムでも確認してみよう

In [61]:
results.resid.hist(color='lightgray')
plt.xlabel('residual')
plt.ylabel('frequency')
Out[61]:
<matplotlib.text.Text at 0x7f935b2368d0>

結果は予想通りに、残差の分布は幾つかの分布から構成されていると考えられる。この期間をさらに短く分割することで、時間トレンドで日経平均株価の動きを説明できる可能性をこのヒストグラムは示唆している

景気 期間始点   終点 ror Vol 係数 切片   残差STD
戦後復興期(recover) 1949/5/16 1954/11/30 14% 23%  0.011 4.52 0.25
高度経済成長期(growth) 1954/12/1 1971/12/31 13% 14% 0.004 6.15 0.23
安定期(stable) 1972/1/1 1986/11/30 13% 13% 0.004 8.10 0.13
バブル経済期(bubble) 1986/12/0 1993/10/31 3% 23% -0.0003 10.28 0.20
バブルのピークまで   1986/12/1   1989/12/31   26% 18% 0.0008 9.9 0.05
バブルのピークから谷 1990/1/1 1992/8/31   -25% 29% -0.0010 10.43 0.09
経済変革期(reform) 1993/11/1 現在   2% 24% -0.0001 9.7 0.26

結果の解釈 

時間の経過に比例して価格が上昇、または下落するという一見単純すぎる時間トレンドモデルのもつ意味を理解いただけましたか?  

資料の例から時間トレンド、または確定的トレンドの有無の判断が統計的に如何に難しいかが理解できます。

決定係数(R-squared)などが良好であっても、残差のチェックが必要である。

また、バブルが崩壊するまでは日経平均株価の上昇率は一定水準を確保し、下落の後に常に回復してきていた。

しかしバブル経済の終焉と共にその傾向は無くなった。

背後にある経済の状態の理解こそがこのモデルの採用の要因であることには間違いがない。

確定的トレンドは人々の経済活動の結果であり、継続がもたらす成果でもあります。

In [ ]: