季節性とマクロ変数を加える?

今まで見てきた時間トレンド、ランダムウォークモデル、ARモデルなどの線形のモデルを統一的に扱うことができるのだろうか? そしてそれに季節性、マクロ変数も加えることができるのだろうか? もちろん可能である。 そして、例えば、  

時系列データ = 時間トレンド+ランダムウォーク+ AR(1)+ 季節性 +マクロ変数

と書き、システム的な方法を取ることでより豊かな表現が可能になる。

特に季節性とマクロ変数との関係は金融市場の効率性とも密接な関係にあることから長い研究の歴史がある。

ただし、同時に変数の数を増やすことで、問題も生じてくる。本章では季節性とマクロ変数の分析を簡単に紹介します。

その前に

ニュートンも解けない3体問題を知る  

$$P_{t} = T_{t} + W_{t} + Z_{t} + \xi _{t}$$

$P_{t}$は観測値 $\xi _{t}$はかく乱項 そして、それぞれのシステムは  

$Y_{t} = a + \beta _{t} + u{t}$ 

$W{t} = a + W{t-1} + w{t}$ 

$Z{t} = a + \beta Z{t-1} + z{t}$

と表現できます。

上から、時間トレンド、ランダムウォーク、AR(1)モデルである。

それぞれの方程式がかく乱項をもつことで、観測値の複雑な挙動を説明できる柔軟性を確保している。

そしてこれらのかく乱項をさらに別の方程式で説明することで、さらなる豊かさが確保される。

季節性とマクロ変数の関係をシステムモデルとして加えてもいいし、残差項を季節性とマクロ変数などで説明しても構わない。

また、それぞれのシステムモデルを動的・静的に扱う選択も可能です。

一般に

1.データに多くのノイズが含まれていないか?  

2.価格の測定単位としての通貨の価値が安定であるか?  

3.モデルは間違っていないか?  

4.残差項が多すぎないか?  

5.現象そのものがランダムではないか?  

などに十分な注意を払うべきです。  

実は高度な数学、統計学を用いても、正確に金融市場の価格の動きを説明できる可能性は少ないです。

現在の金融のモデル化の手法は価格の動きを運動として方程式を立て、それを積分して解を導くという手法を取っている。

しかし、このような手法はニュートン力学の延長線上にあります。

2つの天体の間の運動方程式を積分により解くことで安定周期を得るニュートン力学では天体の数が3つ以上になるとその相互作用はかなり複雑になり解析的に解を得ることはできなくなる。

それを3体問題とよび、変数が増えると解析的には解けないという問題が生じる

季節性の分析  

金融市場の季節性は、学術的に積極的に検証されているものから格言的なものまで多数あります。ここでそれらを整理をしておきます。

*月次効果  

a.ヘッジファンドの決算売り(5月)  

b.米国ファンドの決算売り(6月)  

c.夏枯れ相場(7、8月)  

d.米国の節税売り(10月)  

e.ヘッジファンドの決算売り(11月)  

f.12月の節税売り(12月)

*曜日効果  

a.月曜効果  

b.金曜効果  

c.週末効果

*月末効果  

a.節分天井彼岸底 : 3月末の本決算を前に機関投資家、年金などが利益確定売りに出るため。  

b.掉尾の一振 年末にかけて株価が上昇すること。  

c.サンタクロースラリー クリスマスから年末にかけて株価が上がること。  

d.もちつき相場 : 年末のボラティリティーが上がること  

月末の期間をどの程度の長さに取るかによって、月末効果の評価は大きく変わる。また、月末効果は標本の数が少ないので、結果の解釈には注意が必要である。

*週次効果  

a.第1金曜:米国雇用統計  

b.第2金曜:先物SQ

*祝日相場  

a.ゴールデンウイーク相場:昭和の日の前日からゴールデンウイーク明けにかけて相場が上昇すること。  

b.盆休み相場 :8月10日前後にポジションを調整する動き。  

c.ラマダン相場 :ラマダン期間中は相場は横ばいか下げやすい

平均値の検定  

得られた実現値の平均値がその母集団の平均値と等しいか否かを検定する方法を考えてみよう。

母標準偏差の値は未知である。しかし実現値から算出した標本標準偏差と母集団の標準偏差が等しいとすることで、仮説検定が季節性の判別に使えますね。

例えば、バブル崩壊前の1989年1月の価格から計算した終値の変化率の母平均について判定してみます。  

仮説検定では2つの仮説を立てる。ひとつは証明したい仮説の反対を示すもので帰無仮説(Null Hypothesis)といい$H_{0}$と書きます。

帰無仮説は棄却を前提に立てられる。 もう一方は対立仮説で帰無仮説に反する仮説で$H_{1}$と書き

例えば、1月効果ではその効果の棄却を前提に帰無仮説を立てます。

母集団の平均値$\mu $が想定される値$\mu _{0} $に等しいとします。

そして、対立仮説では母平均$\mu $は$\mu _{0} = 0$より大きいと一般には3つの基本形があります。

1. 母平均はプラス  

$a. H_{0} : \mu = 0$ 

$b. H_{1} : \mu > 0 $

2.母平均はマイナス  

$a. H_{0} : \mu = 0$ 

$b. H_{1} : \mu < 0 $

3.母平均はゼロ  

$a. H_{0} : \mu = 0$ 

$b. H_{1} : \mu \neq 0 $

 これらの仮説が実際の観察値と整合性があるかどうかの判定は、統計的仮説の検定(The Test of Significant Appoarch)を用いて行います。

統計量Tは

$$T = \frac{ \overline{ \mu } }{ \frac{8}{ \sqrt{n} } } $$

で定義されます。

そしてその実現値をtで表します。  

統計的仮説の検定では実現値tが採択域にいれば帰無仮説を棄却することなく、また棄却域にいれば帰無仮説を棄却すると判断できます。

表にまとめておきますね。

平均値 仮説のタイプ 帰無仮説 対立仮説 決定規則
プラス 右片側 μ= 0 u>0 t>tx,df
マイナス 左片側 μ= 0 u<0 t<-tx,df
ゼロ 両側 μ≠ 0 μ≠ 0 t>tx/2,df

表で示した決定規則が成り立てば、帰無仮説を棄却する。1989年の1月の変化率の標本平均はμ=0.00242、標本標準偏差はs =0.00514、標本数はn=18である。

$H_{0} : \mu = 0、H_{1} : \mu > 0$とすると、帰無仮説が棄却されれば1月の平均値がプラスである可能性があります。

$t \geq t{x,n-1}$ならば$H_{0}$を棄却し、 $t < t_{x,n-1}$ならば$H_{0}$を採択する。

ここでxは有意水準である。計算してみよう。t=1.949。 $t_{0.05,17} = 1.74$であるから、$H_{0}$は棄却され、1989年の1月の平均値がプラスの可能性が否定できないです。

季節性の具体例  

標本から推定した母平均μがゼロに等しいか否かで季節性を検定してみます。

その際には、標本標準偏差を用いるので、分散不均一性の問題が生じる。この問題を避けるために1年間毎に分けて実現値を分析し、分析対象期間でどれくらいの割合で、帰無仮説が棄却されるかを見てみます。  

月次の季節性  

それぞれの月の間に起こる値動きに明確な特徴があるかどうかを過去にさかのぼり調べてみました。

例えば1月効果であれば、1949年5月16日から2017年8月3日までの各年の1月の日々の変化率からその平均値がプラスであるか否か有意水準10%の統計的仮説の検定を用いて判断し、

そして、帰無仮説が棄却された年の数で、季節性の判定をしてみました。

In [6]:
from scipy.stats import t
import pandas as pd
import pandas_datareader.data as web
import numpy as np
end='2017/8/3'
n225 = web.DataReader("NIKKEI225", 'fred',"1949/5/16",end).dropna()
develop=n225.loc[:'1989/12/31']
reform=n225.loc['1989/12/31':]
year=n225.loc['1989']
years=[x+1950 for x in range(66)]
m=lambda x:x.month
count=[0]*12
alpha=0.1
for i in range(len(years)):
    year=n225.loc[str(years[i])]
    r=year.pct_change().groupby([m])
    tv=r.mean()/r.std()*np.sqrt(r.count())
    t0=t.ppf(1-alpha,len(r)-1)
    for j in range(12):
        if float(tv.iloc[j])>t0:# and years[i]>=1990:
            count[j]+=1
print(count)
[24, 15, 14, 9, 14, 10, 9, 14, 16, 9, 14, 15]
In [8]:
from scipy.stats import t
import pandas as pd
import pandas_datareader.data as web
import numpy as np
m=lambda x:x.month
count=[0]*12
alpha=0.1
for i in range(len(years)):
    year=n225.ix[str(years[i])]
    r=year.pct_change().groupby([m])
    tv=r.mean()/r.std()*np.sqrt(r.count())
    t0=t.ppf(1-alpha,len(r)-1)
    for j in range(12):
        if float(tv.iloc[j])>t0 and years[i]>=1990:
            count[j]+=1
print(count)
print(t0)
[2, 3, 2, 2, 4, 2, 2, 1, 5, 1, 3, 4]
1.36343031802

[24, 15, 14, 9, 14, 10, 9, 14, 16, 9, 14, 15] 全期間の結果

[2, 3, 2, 2, 4, 2, 2, 1, 5, 1, 3, 4] 1990年以降の結果

結果は左から順に1月、2月、…、12月となる。平均値がプラスであると判定された月は1950年から2017年までの67年間で1月では24回あります。

その内1990年以降に起きたのは2回であり、1月効果はバブル崩壊前の現象であった可能性が高い。39年間の内に22回起きているので、1月にプラスになる確率は高い。

バブル崩壊以降は、どの月においても明確な方向性はなく、多くの月で平均値ゼロと判断される。

プログラムコードの解説:無名関数ラムダ 月次、還次のデータの抽出には無名関数lambdaを用いる。

月の最初の日の変化率はその前の月の最終日の価格があれば、 それを用いて計算する。

データが月の初めから始まれば、最初の日の変化率はNaNが返される。

y=lambda x:x,month 無名関数による月の設定 r=year.pct_change().groupby([y]) グループ化

プログラムコードの解説:年次データの取得 for i in range(len(year)): year=n225.ix[str(year[i]) 1年間のデータの抽出

プログラムコードの解説:月次の推定量からt値の算出 from scipy.stats import t ステューデントの1分布のインポート

tv=r.mean()/r.std()*np.sqrt(r.count()) t値の算出 t0=t.ppf(0.90,len(r)-1) t分布の臨界値 for j in range(12): 月の選択 If float(tv.iloc[j])>t0: 判定 count+=1 H0 棄却数

また、同様の分析をマイナスの変化率に対しても行ってみました。その結果は

In [9]:
from scipy.stats import t
import pandas as pd
import pandas_datareader.data as pdr
import numpy as np
m=lambda x:x.month
count=[0]*12
for i in range(len(years)):
    year=n225.ix[str(years[i])]
    r=year.pct_change().groupby([m])
    tv=r.mean()/r.std()*np.sqrt(r.count())
    t0=-t.ppf(1-alpha,len(r)-1)
    for j in range(12):
        if float(tv.iloc[j])<t0:# and years[i]>=1990:
            count[j]+=1
print(count)
[2, 5, 5, 2, 9, 4, 6, 6, 5, 4, 7, 3]
In [10]:
from scipy.stats import t
import pandas as pd
import pandas_datareader.data as pdr
import numpy as np
m=lambda x:x.month
count=[0]*12
for i in range(len(years)):
    year=n225.ix[str(years[i])]
    r=year.pct_change().groupby([m])
    tv=r.mean()/r.std()*np.sqrt(r.count())
    t0=-t.ppf(1-alpha,len(r)-1)
    for j in range(12):
        if float(tv.iloc[j])<t0 and years[i]>=1990:
            count[j]+=1
print(count)
[1, 2, 2, 0, 4, 2, 2, 2, 3, 0, 3, 1]

[2, 5, 5, 2, 9, 4, 6, 6, 5, 4, 7, 3]全期間の結果

[1, 2, 2, 0, 4, 2, 2, 2, 3, 0, 3, 1]1990年以降の結果

です。変化率がプラスの検定の場合と同様に、左端が1月であり右端が12月である。1月効果のような現象は見られない。

12月の節税売りも見られない。あえていえば5月のヘッジファンド売りがあるかもしれないが5月は収益率がプラスになる可能性も高い

曜日の季節性  

次に曜日について分析してみました。月曜効果、金曜効果、週末効果が注目されるが、本分析では金曜の終値から月曜の終値までの分析となるので、月曜効果と週末効果が複合して分析されている

In [11]:
w=lambda x:x.week
count=[0]*5
for i in range(len(years)):
    year=n225.ix[str(years[i])]
    r=year.pct_change().groupby([w])
    tv=r.mean()/r.std()*np.sqrt(r.count())
    t0=t.ppf(1-alpha,len(r)-1)
    for j in range(5):
        if float(tv.iloc[j])>t0:
            count[j]+=1
print(count)
[15, 17, 14, 13, 10]

65年間で月曜が15回、火曜が17回、水曜が14回、木曜が13回、金曜が10回の平均値が有意にプラスであると判断された。

月次効果と同じように、ほとんどの有意な曜日をもつ年はバブル崩壊前に起きていて、バブル崩壊後には曜日には全く特徴が見られない。

プログラムコードの解説:lambda関数、曜日の設定

python y=lambda x:x,weekday 無形関数による曜日の設定

変化率のマイナス側の結果は次の通りである

In [12]:
w=lambda x:x.week
count=[0]*5
for i in range(len(years)):
    year=n225.ix[str(years[i])]
    r=year.pct_change().groupby([w])
    tv=r.mean()/r.std()*np.sqrt(r.count())
    t0=-t.ppf(1-alpha,len(r)-1)
    for j in range(5):
        if float(tv.iloc[j])<t0:
            count[j]+=1
print(count)
[7, 10, 5, 3, 7]
In [13]:
w=lambda x:x.week
count=[0]*5
for i in range(len(years)):
    year=n225.ix[str(years[i])]
    r=year.pct_change().groupby([w])
    tv=r.mean()/r.std()*np.sqrt(r.count())
    t0=t.ppf(1-alpha,len(r)-1)
    for j in range(5):
        if float(tv.iloc[j])>t0 and years[j]>1990:
            count[j]+=1
print(count)
[0, 0, 0, 0, 0]

[7, 10, 5, 3, 7]全期間の結果

[0, 0, 0, 0, 0]1990年以降の結果

特に注目する結果は見られない

マクロ変数との関係(単回帰と多変量解析)  

今まで、内閣府の景気循環期をもとに日経平均株価を分析してきた。

ここでマクロ変数と日経平均株価の関係を考えてみよう。

誰にとってもなじみの深い経済変数は、外国為替レート、国内総生産、金利、物価水準、労働人口、マネーサプライなどではないだろうか? 

ここではFred(セントルイス連銀サイト)からできるだけ長い期間のデータでダウンロードできる変数を3つ選びました。

1つは労働人口、国内総生産、そしてドル円の為替レートです。

国内総生産(出所:OECD)の指数はドル建てであるために、円建てに直し

ドル円の為替レートはドル建ての指数を円建てに変換するために用いた。日経平均株価を労働人口と国内総生産の国内要因のみで説明するモデルを考えてみました。

説明変数と被説明変数の変数すべてに対数をとるので、モデルは対数線形モデルとよばれます。

この場合にはある説明変数の傾きの係数はその説明変数に対する被説明変数の弾力性を表しています。

単回帰分析  

日経平均を分析する前に手始めとして、労働人口と国内総生産について調べてみました。

労働人口は経済の原動力であり、国内総生産はその結果です

In [12]:
import pandas as pd
import pandas_datareader.data as web
import numpy as np
start='1971/12/1'
end='2017/8/3'
workpop = web.DataReader('LFWA64TTJPM647S',"fred",start,end).dropna()
gdp = web.DataReader('MKTGDPJPA646NWDB',"fred",start,end).dropna()
gdp=gdp.resample('A',loffset='-1d').last().dropna()
fx = web.DataReader('DEXJPUS',"fred",start,end).dropna()
fx=fx.resample('A',loffset='-1d').last().dropna()
workpop=workpop['1972':].resample('A',loffset='-1d').last().dropna()
gdpjpy=gdp.MKTGDPJPA646NWDB*fx.DEXJPUS
gdpjpy=np.log(gdpjpy).dropna()
workpop=np.log(workpop).dropna()
In [13]:
import statsmodels.api as sm
x=sm.add_constant(gdpjpy)
model=sm.OLS(gdpjpy,x)
results=model.fit()
print(results.summary())
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       1.000
Model:                            OLS   Adj. R-squared:                  1.000
Method:                 Least Squares   F-statistic:                 3.379e+27
Date:                Fri, 04 Aug 2017   Prob (F-statistic):               0.00
Time:                        01:53:25   Log-Likelihood:                 1311.3
No. Observations:                  45   AIC:                            -2619.
Df Residuals:                      43   BIC:                            -2615.
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const       3.411e-13   5.77e-13      0.591      0.558   -8.23e-13     1.5e-12
0              1.0000   1.72e-14   5.81e+13      0.000       1.000       1.000
==============================================================================
Omnibus:                       14.676   Durbin-Watson:                   0.001
Prob(Omnibus):                  0.001   Jarque-Bera (JB):               16.013
Skew:                          -1.378   Prob(JB):                     0.000333
Kurtosis:                       3.971   Cond. No.                     2.37e+03
==============================================================================

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

ドル建てのGDPは円換算して対数をとっています。

これを被説明変数、労働人口を説明変数として回帰を行い、GDPのデータは年次で、他のデータは日次、月次であるために、GDPのデータに期間を合わせています。

切片、回帰係数ともにp-値がゼロとなっています。

$R^{2}$も0.586と高い数値となりました。

つぎに期待値(fitted)と実現値をグラフで比べてみます。

In [8]:
%matplotlib inline 
import matplotlib.pyplot as plt
f,ax = plt.subplots()#2軸のグラフの準備
ax.plot(gdpjpy,label='gdp',linestyle="--")
ax2=ax.twinx()#2軸目をax2として設定
ax2.plot((workpop),label='workpop')#2軸目にプロット
results.fittedvalues.plot(label='fitted',style=':',ax=ax)
ax.set_ylabel('log GDP')#1軸目にラベルを設定
ax2.set_ylabel('workshop')#2軸目にラベルを設定
ax.legend(loc='lower right')
ax2.legend(loc='upper left')
Out[8]:
<matplotlib.legend.Legend at 0x7fcf7c916f60>

期待値は点線で示さている。労働人口(workpop)はグラフでの比較が容易になるように2軸目に表示してあることに注意してほしいです。

2000年以降、期待値は実現値をうまく説明できていない。労働生産性のような別の要素が影響をしているのかもしれない

ランダムウォークの悪魔  

つぎに残差についてみてみよう。一般に複数の非定常な時系列を最小二乗法を用いて回帰する際には、関係の無い2つの変数の間の回帰係数が有意な推定値をもつと判断されてしまう可能性がある。

これを、見せかけの回帰(spurious regression)とよび、単位根問題の1つとして長い研究の歴史があります。

日経平均株価、国内総生産、労働人口はランダムウォークであるので、見せかけの回帰を常に頭に置いておく必要がありますね

In [9]:
import matplotlib.pyplot as plt
results.resid.hist(label='residual',color='seagreen')
plt.xlabel('residual:gdp vs work population')
plt.ylabel('frequency')
plt.legend(loc='upper right')
Out[9]:
<matplotlib.legend.Legend at 0x7fcf78414a20>

直感的には、残差は正規分布とみなせると思うが、標本数は多くないので確認のために、Jarque-Bera(JB)の正規性の検定を紹介しますね。

Jarque-Bara検定  

JBでは標本数がnの残差からその歪度(S)と尖度(K)を求め、次の検定統計量を算出し

$$JB = n[ \frac{S^2}{6} + \frac{(K-3)^2}{24}] $$

残差の正規性を検定しますね。統計量は自由度2の$\chi ^{2}$乗分布に従います。帰無仮説は残差の正規性であり、得られたp-値(両側検定)が十分に小さければ帰無仮説を棄却し、そうでなければ採択する。  

Statsmodel.OLSのsummary()の中にはJarque-Bera(JB)が含まれている。そのp-値は0.00297と非常に小さい。したがって残差は正規分布に従っているとはいえない

プログラムコードの解説:データの取得と準備

workpop = pdr.DataReader('LFWA64TTJPM647S',"fred",start,end).dropna()

労働人口の取得

gdp = pdr.DataReader('MKTGDPJPA646NWDB',"fred",start,end).dropna()

ドル建て国内総生産の取得

gdp=gdp.resample('A',loffset='-1d').last().dropna()

データの日付が1月1日であるので、これを12月30日に変更

fx = pdr.DataReader('DEXJPUS',"fred",start,end).dropna()

ドル/円データの取得

fx=fx.resample('A',loffset='-1d').last().dropna()

ドル/円のデータを年次に変換し、12月30日に一番近い日のデータを用いる。

workpop=workpop['1972':].resample('A',loffset='-1d').last().dropna()

労働人口のデータを年次に変換し12月30日に近いデータを使うためにlast()を用した。

'A'は年次にデータを変換するか日付は年初になってしまう。それでoffset='-ld'を用いて

日付を1日前にずらしている。

x=sm.add_constant(workpop) 切片のために列の追加

model=sm.OLS(gdpjpy,x) 回帰モデルの設定

results=model.fit() モデルの最適化と結果の格納

print(results.summary()) 結果の要約を出力

多変量解析  

つぎに日経平均株価を被説明変数、労働人口と国内総生産を説明変数として多回帰分析を試みてみよう。日経平均株価も対数をとってあることを忘れないでほしいです。

In [36]:
import pandas as pd
lnn225 = np.log(pdr.DataReader("NIKKEI225", 'fred',start,end).dropna())
lnn225=lnn225.resample('A',loffset='-1d').last().dropna()
port=pd.concat([lnn225,x,gdpjpy],axis=1).dropna()
port.columns=["n225","const","workpop","gdpjpy"]
model=sm.OLS(port.n225,port.iloc[0:,1:])
results=model.fit()
print(results.summary())
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                   n225   R-squared:                       0.552
Model:                            OLS   Adj. R-squared:                  0.542
Method:                 Least Squares   F-statistic:                     53.09
Date:                Fri, 04 Aug 2017   Prob (F-statistic):           4.97e-09
Time:                        02:32:52   Log-Likelihood:                -18.712
No. Observations:                  45   AIC:                             41.42
Df Residuals:                      43   BIC:                             45.04
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const        -19.4310      3.956     -4.912      0.000     -27.409     -11.453
workpop        0.4296      0.059      7.286      0.000       0.311       0.549
gdpjpy         0.4296      0.059      7.286      0.000       0.311       0.549
==============================================================================
Omnibus:                        6.742   Durbin-Watson:                   0.335
Prob(Omnibus):                  0.034   Jarque-Bera (JB):                5.702
Skew:                           0.837   Prob(JB):                       0.0578
Kurtosis:                       3.489   Cond. No.                     7.57e+17
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The smallest eigenvalue is 1.77e-31. This might indicate that there are
strong multicollinearity problems or that the design matrix is singular.

$R^{2}$は0.552です。また、それぞれの回帰係数のp-値もどれも1%以下です。

また、JBのp-値は1%以下。結果を目で見で確認してみます。

日経平均株価の対数とその期待値をプロットして、モデルの期待値がどれほど実際の指数の動きを説明できるかみてみよう

In [22]:
#多変量解析:折れ線グラフ
f,ax = plt.subplots()#2軸のグラフの準備
(port.gdpjpy-24).plot(label='gdp',linestyle="--",ax=ax)
port.n225.plot(label='n225',ax=ax)
ax2=ax.twinx()#2軸目をax2として設定
(port.workpop).plot(label='workpop',ax=ax2,style='o')
results.fittedvalues.plot(label='fitted',style=':',ax=ax)
plt.legend(loc='upper left')
ax.set_ylabel('log Nikkei225 index')#1軸目にラベルを設定
ax2.set_ylabel('workshop')#2軸目にラベルを設定
ax.legend(loc='lower right')
ax2.legend(loc='upper left')
Out[22]:
<matplotlib.legend.Legend at 0x7fa85b4e22e8>

バブル成長期に日経平均株価がオーバーシュートし、バブル崩壊により期待値に接近している。

また、バブル崩壊後には米国のインターネットバブル崩壊、日本の第3次平成不況期に期待値を下方にオーバーシュートし、その後回復している。

また、リーマンショックで期待値を下回る様子が見てとれる。  

JBは0.203で正規分布ではないとはいいきれない状態にある。一応、目視で残差を確認しておきます。

In [23]:
#多変量解析:ヒストグラム
results.resid.hist(label='residual',color='mistyrose')
plt.xlabel('residual:gdp vs work population')
plt.ylabel('frequency')
plt.legend(loc='upper right')
Out[23]:
<matplotlib.legend.Legend at 0x7fa860cd2a90>

経済の構造変化

日経平均株価、国内総生産、労働人口がバブルの崩壊を機にそれ以前の上昇から、下落のトレンドに転換したようにみえるので、それを確かめてみよう。分析はバブル崩壊前と後に2つ時期に分けて行った

In [25]:
#バブル崩壊前のグラフ
f,ax = plt.subplots()#2軸のグラフの準備
(port[:'1990/1/1'].gdpjpy-24).plot(label='gdp',linestyle="--",ax=ax)
port[:'1990/1/1'].n225.plot(label='n225',ax=ax)
ax2=ax.twinx()#2軸目をax2として設定
(port[:'1990/1/1'].workpop).plot(label='workpop',style='o',ax=ax2)
#書籍のグラフはworkpopから8.5引いてしまっているので、こちらが正しいグラフです。
results_b.fittedvalues.plot(label='fitted',style=':',ax=ax)
ax.set_ylabel('log Nikkei225 index')#1軸目にラベルを設定
ax2.set_ylabel('workshop')#2軸目にラベルを設定
ax.legend(loc='lower right')
ax2.legend(loc='upper left')
Out[25]:
<matplotlib.legend.Legend at 0x7fa85b292f60>

JBで0.612となった残差の状態をグラフで確認しよう

In [21]:
#バブル崩壊前:ヒストグラム
results_b.resid.hist(label='residual',color='seagreen')
plt.xlabel('residual:gdp vs work population')
plt.ylabel('frequency')
plt.legend(loc='upper right')
Out[21]:
<matplotlib.legend.Legend at 0x7fcf782f03c8>

バブル崩壊後  

つぎにバブル崩壊後について崩壊前と同じ分析を行ってみよう。結果はつぎのとおりです。

In [22]:
#バブル崩壊後
port_a=port['1990/1/1':]
results_a=(sm.OLS(port_a.n225,port_a.iloc[0:,1:])).fit()
print(results_a.summary())
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                   n225   R-squared:                       0.019
Model:                            OLS   Adj. R-squared:                 -0.020
Method:                 Least Squares   F-statistic:                    0.4820
Date:                Fri, 04 Aug 2017   Prob (F-statistic):              0.494
Time:                        02:02:57   Log-Likelihood:                -5.9628
No. Observations:                  27   AIC:                             15.93
Df Residuals:                      25   BIC:                             18.52
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         -6.3701     22.986     -0.277      0.784     -53.711      40.971
workpop        0.2356      0.339      0.694      0.494      -0.463       0.935
gdpjpy         0.2356      0.339      0.694      0.494      -0.463       0.935
==============================================================================
Omnibus:                        1.826   Durbin-Watson:                   0.511
Prob(Omnibus):                  0.401   Jarque-Bera (JB):                1.108
Skew:                          -0.132   Prob(JB):                        0.575
Kurtosis:                       2.044   Cond. No.                     1.22e+19
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The smallest eigenvalue is 4.17e-34. This might indicate that there are
strong multicollinearity problems or that the design matrix is singular.

$R^{2}$、AIC、BICともに悪化している。回帰係数のp-値もすべて20%以上である。これはどの説明変数も被説明変数の動きをとらえていない可能性を示唆しているのである。

このように経済はある時点を境に大きく変化する可能性がある。  

そこで説明変数としてドル円の為替レートを加えてみよう。特にバブル崩壊後は強い円高に日本経済は苦しめられてきた

In [26]:
import pandas as pd
lnn225 = np.log(pdr.DataReader("NIKKEI225", 'fred',start,end).dropna())
lnn225=lnn225.resample('A',loffset='-1d').last().dropna()
lnfx=np.log(fx)
port1=pd.concat([lnn225,x,gdpjpy,lnfx],axis=1).dropna()
port1.columns=["n225","const","workpop","gdpjpy","fx"]
model1=sm.OLS(port1.n225,port1.iloc[0:,1:])
results1=model1.fit()
print(results1.summary())
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                   n225   R-squared:                       0.555
Model:                            OLS   Adj. R-squared:                  0.534
Method:                 Least Squares   F-statistic:                     26.22
Date:                Fri, 04 Aug 2017   Prob (F-statistic):           4.07e-08
Time:                        02:11:48   Log-Likelihood:                -18.571
No. Observations:                  45   AIC:                             43.14
Df Residuals:                      42   BIC:                             48.56
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const        -24.0115      9.763     -2.459      0.018     -43.715      -4.308
workpop        0.4865      0.126      3.874      0.000       0.233       0.740
gdpjpy         0.4865      0.126      3.874      0.000       0.233       0.740
fx             0.1545      0.300      0.514      0.610      -0.452       0.761
==============================================================================
Omnibus:                        7.804   Durbin-Watson:                   0.358
Prob(Omnibus):                  0.020   Jarque-Bera (JB):                6.884
Skew:                           0.921   Prob(JB):                       0.0320
Kurtosis:                       3.530   Cond. No.                     7.72e+17
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The smallest eigenvalue is 1.72e-31. This might indicate that there are
strong multicollinearity problems or that the design matrix is singular.
In [27]:
#バブル崩壊後:要素にドル円の為替レートを追加
port1_a=port1['1990/1/1':]
results1_a=(sm.OLS(port1_a.n225,port1_a.iloc[0:,1:])).fit()
print(results1_a.summary())
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                   n225   R-squared:                       0.324
Model:                            OLS   Adj. R-squared:                  0.268
Method:                 Least Squares   F-statistic:                     5.753
Date:                Fri, 04 Aug 2017   Prob (F-statistic):            0.00910
Time:                        02:12:06   Log-Likelihood:               -0.93333
No. Observations:                  27   AIC:                             7.867
Df Residuals:                      24   BIC:                             11.75
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          4.7868     19.766      0.242      0.811     -36.008      45.581
workpop       -0.0144      0.297     -0.048      0.962      -0.628       0.599
gdpjpy        -0.0144      0.297     -0.048      0.962      -0.628       0.599
fx             1.2308      0.374      3.292      0.003       0.459       2.002
==============================================================================
Omnibus:                        2.948   Durbin-Watson:                   0.792
Prob(Omnibus):                  0.229   Jarque-Bera (JB):                1.802
Skew:                          -0.619   Prob(JB):                        0.406
Kurtosis:                       3.260   Cond. No.                     1.10e+19
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The smallest eigenvalue is 5.2e-34. This might indicate that there are
strong multicollinearity problems or that the design matrix is singular.

ドル円の為替レートを加えたことで大きな変化が現れた。それはそれぞれの回帰係数のp-値である。

切片、国内総生産、そして労働人口の回帰係数が大幅に上昇し、一方でドル円の為替レートの回帰係数は5%以下に収まっている。

また、$R^{2}$、AIC、BICともに前のモデルよりかは改善している。しかしJBに関しては若干の悪化がみられる。

モデルが示している期待値と実際の日経株価平均、ドル円の為替レートをプロットしてみよう

In [26]:
##バブル崩壊後のグラフ:ドル円の為替レート
(port1['1990/1/1':].fx+5).plot(label='fx',linestyle="--")
port1['1990/1/1':].n225.plot(label='n225')
results1_a.fittedvalues.plot(label='fitted',style=':')
plt.ylabel('log Nikkei225 Index')
plt.legend(loc='lower left')
Out[26]:
<matplotlib.legend.Legend at 0x7fe2f1e272b0>

グラフからはリーマンショック以降に期待値は実際の動きをよく説明している。

経済は生き物  

バブル崩壊後の経済をもう少し細かくみてみよう。これらは過去を振り返った分析であるので注意が必要である。

期間は1990/1/1から2000/1/1、2000/1/1から2008/1/1、2008/1/1以降の3つに分割した。

結果をつぎの表にまとめました。

In [27]:
#バブル崩壊後:細分化
def report(port):
    results1_a=(sm.OLS(port1_a.n225,port1_a.ix[0:,['const','gdpjpy','workpop','fx']]))\
    .fit()
    print("R-squared: ",results1_a.rsquared," F-pvalue: ",results1_a.f_pvalue," AIC: "\
          ,results1_a.aic," BIC: ",results1_a.bic)
    print("pvalues: ")
    print(results1_a.pvalues)
    from statsmodels.compat import lzip
    import statsmodels.stats.api as sms
    test=sms.jarque_bera(results1_a.resid)
    print("jbpv: ",test[1])
port1_a=port1['1990/1/1':'2000/1/1']
report(port1_a)
R-squared:  0.637234115352  F-pvalue:  0.0889983480357  AIC:  -10.4131864571  BIC:  -9.20284608516
pvalues: 
const      0.023541
gdpjpy     0.254732
workpop    0.026042
fx         0.095029
dtype: float64
jbpv:  0.651686690925
In [28]:
#バブル崩壊後:細分化2
port1_a=port1['2000/1/1':'2008/1/1']
report(port1_a)
R-squared:  0.926156045386  F-pvalue:  0.00996898531012  AIC:  -13.7258843227  BIC:  -13.408118156
pvalues: 
const      0.061178
gdpjpy     0.006358
workpop    0.007690
fx         0.057891
dtype: float64
jbpv:  0.731005909439
In [29]:
#バブル崩壊後:細分化3
port1_a=port1['2008/1/1':]
report(port1_a)
R-squared:  0.949635785933  F-pvalue:  0.00467542454257  AIC:  -12.4292010344  BIC:  -12.1114348677
pvalues: 
const      0.182073
gdpjpy     0.726831
workpop    0.122300
fx         0.050233
dtype: float64
jbpv:  0.600603318664
1990以前 1990以降 1990ー2000 2000ー2008 2008以降
R2 0.94 0.32 0.69 0.93 0.95
Fp-値 0.00 0.03 0.05 0.01 0.00
AIC -7.45 10.18 -12.18 -13.76 -12.37
BIC -4.78 15.22 -10.92 -13.45 -12.05
回帰係数 p-値
切片 0.00 0.99 0.01 0.05 0.17
国内総生産(gdp) 0.02 0.99 0.13 0.01 0.76
労働人口(workpop) 0.00 0.87 0.01 0.01 0.11
ドル/円(fx) - 0.01 0.06 0.06 0.05
JB p-値 0.61 0.47 0.84 0.73 0.60
                                                         ※fitted=期待値 n225=日経平均株価

表からわかるようにそれぞれの期間で説明変数のp-値は大きく異なり、被説明変数に与える影響が異なるという結果になっている。

1990-2000の期間では労働人口とドル円の為替レートが日経平均株価の有力な説明変数になっている一方、2000-2008の期間ではそれに国内総生産が加わる。

2008以降ではドル円の為替レートだけが有力な説明変数である。これは期間の取り方により微妙に変化する

プログラムコードの解説:report関数

print("R-squared: ",results1_a.rsquared," F-pvalue: ",results1_a.f_pvalue," AIC: "\ ,results1_a.aic," BIC: ",results1_a.bic)

演算結果の受け渡し

R-squared: self.rsqaured

F-p-値: self.f_pvalue

AIC: self.aic

BIC: self.bic

print("pvalues: ")

print(results1_a.pvalues)

回帰係数の出力:self.pvalues

from statsmodels.compat import lzip

import statsmodels.stats.api as sms

test=sms.jarque_bera(results1_a.resid)

JBに関してはOLS.fit()で演算されているにもかかわらず、

取り出すことができないので別途計算した。演算結果をtestに格納

print("jbpv: ",test[1])

[0]は統計量、[1]はp-値である。

多変量解析の原理  

ここで用いた多変量解析のモデルは、実は自己回帰移動平均(ARMA)モデルで使用したものと同じ原理です。

Material8での編自己相関を思い出してほしい。

自己回帰モデルAR(n)とは実現値となる被説明変数がn個の自己の過去の変数の和のモデルであった。

このように複数の要因(確率変数)があり、2つの特定の組となる変数の間の相関を求めたいときに、他の確率変数の動きを固定して、編自己相関なるものを

求めた。これは条件付確率とよばれる方法を用いている。同じ方法が多変量解析にも用いられているのです。

全般を通して 

長期投資の収益源である確定的トレンドの存在と短期の収益源である季節性の存在が、統計的には必ずしも有意だとは限らない事実をみてみよう。  

確定的トレンドや季節性、マクロ経済指標との関係などの判断は統計的手法、数学的モデルだけではなく、その現象の背後にある事実関係を知る必要があります。

統計的にモデルの正当性が今説明できなくとも、現実を見て人々が経済の発展(回復)をもたらす強い政策と意思と努力があるのだと確信できればそのモデルは正しくなる。

これらすべての問題点を認識した時点が出発点です。

In [ ]: