日経平均株価とマクロ変数との関係について分析してみた

学習記録

　こんにちは。システムトレーダーの卵ことKenKenです。今回も「Python3ではじめるシステムトレード ──環境構築と売買戦略」で勉強したことを簡単にまとめたいと思います。今回のテーマは「マクロ変数」です！これまでは日経平均株価にのみ着目して分析してきました。今回は、日経平均株価とマクロ変数の関係を見てみたいと思います。使用するマクロ変数は、労働人口、国内総生産にします。データ取得は、pandas-datareaderからFredより取得します。国内総生産は、ドル建てになっているため、円換算にするためにドル円の為替レートも併せて使用します。

マクロ変数の関係性について

　まずは、今回使用する労働人口と国内総生産について調べてみる。労働人口は経済の原動力であり、国内総生産はその結果である。２変数について回帰分析してみた結果は以下の通り。

import pandas_datareader as pdr
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm

# データ取得
start = '1971/12/31'
end = '2016/8/31'
WorkPop = pdr.DataReader('LFWA64TTJPM647S', 'fred', start, end).dropna() # 労働人口(月次データ:1972/1/1, 1972/2/1,...)
GDP = pdr.DataReader('MKTGDPJPA646NWDB', 'fred', start, end).dropna() # 国内総生産（年次データ:1972/1/1, 1973/1/1）
fx = pdr.DataReader('DEXJPUS', 'fred', start, end).dropna() # ドル円為替レート(日次データ:1971/12/31, 1972/1/3, ...)
# GDPに期間を合わせる(年次にする、日付はすべて年末にする)
GDP = GDP.resample('A').last().dropna() # 1972/1/1⇒1972/12/31
WorkPop = WorkPop.resample('A').last().dropna() # 1972/12/1⇒1972/12/31
fx = fx.resample('A').last() # 1972/12/29⇒1972/12/31
# GDPを円換算
GDP_J = GDP.MKTGDPJPA646NWDB * fx.DEXJPUS
# 対数変換
GDP_J = np.log(GDP_J).dropna()
WorkPop = np.log(WorkPop).dropna()

# 回帰分析
x = sm.add_constant(WorkPop)
model = sm.OLS(GDP_J, x)
result = model.fit()
print(result.summary())

結果として、切片、回帰係数ともにp値がゼロとなっているため意味のある変数となった。また、決定係数は0.529とそこそこ高い値となった。グラフを書いてみると、以下の通り。

# グラフを書く
f, ax = plt.subplots() # 2軸グラフ
# 1軸目
ax.plot(GDP_J, label='GDP', linestyle='--')
result.fittedvalues.plot(label='fitted', style=':', ax=ax)
ax.set_ylabel('log GDP')
ax.legend(loc='lower right')
# 2軸目
ax2 = ax.twinx()
ax2.plot(WorkPop, label='WorkPop')
ax2.set_ylabel('log WorkPop')
ax2.legend(loc='upper left')

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

　次に残差についてみてみる。一般的に複数の非定常な時系列を最小二乗法を用いて回帰する際には、関係のない２つの変数間の回帰係数が有意な推定値をもつと判断されてしまう可能性がある。これを見せかけの回帰と呼ぶ。日経平均株価、国内総生産、労働人口はランダムウォークであるので見せかけの回帰については常に頭に置いておく必要がある。

# 残差のヒストグラムを描く
result.resid.hist(label='residual')
plt.xlabel('residual: GDP vs Work Population')
plt.ylabel('frequency')
plt.legend(loc='upper right')

直感的には、残差は正規分布とみなせるような気もするが、標本数が多く無い為、Jarque-Bera(JB)の正規性の検定を行ってみる。

Jarque-Bara検定

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

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

統計量は自由度2の\(\chi^2\)分布に従う。帰無仮説は残差の正規性であり、p値（両側検定)が十分に小さければ帰無仮説は棄却される。

ありがたいことに、Pythonのstatsmodel.OLSのsummary()の中に、Jarque-Bera(JB)が含まれている。上で載せたサマリーを見てみると、p値は0.0282と小さい。したがって、残差は正規分布に従っているとは言えない。

マクロ変数による多変量解析

　次に、日経平均株価を被説明変数、労働人口と国内総生産を説明変数として多変量解析を行ってみる。コードは、以下の通り。

# データ取得
N225 = pdr.DataReader('NIKKEI225', 'fred', start, end).dropna()
ln225 = np.log(N225)
# 日次⇒年次
ln225 = ln225.resample('A').last().dropna()

# OLS
x = pd.concat([WorkPop, GDP_J], axis=1)
x.columns = ['WorkPop', 'GDP']
x = sm.add_constant(x)
y = ln225
model = sm.OLS(y, x)
result = model.fit()
print(result.summary())

\(R^2\)は0.610であり、それぞれの回帰係数のp値はどれも1％未満である。また、JBのp値は0.250となった。次に、日経平均株価の対数と期待値をプロットして、モデルの期待値がどれほど実際の指数を説明できているかを見てみる。

# グラフを書く
f, ax = plt.subplots() # 2軸グラフ
# 1軸目
ax.plot(ln225, label='N225')
result.fittedvalues.plot(label='fitted', style='o', ax=ax)
ax.plot(x['GDP']-24, label='GDP', linestyle='--') # 調整
ax.set_ylabel('log Nikkei225 Index')
ax.legend(loc='lower right')
# 2軸目
ax2 = ax.twinx()
ax2.plot(x['WorkPop'], label='WorkPop', linestyle=':')
ax2.legend(loc='upper left')

　バブル成長期に日経平均株価はオーバーシュートし、バブル崩壊により期待値に接近している様子が見られる。また、バブル崩壊後には米国のインターネットバブル崩壊、日本の第三次平成不況期に期待値を下方にオーバーシュートし、その後回復している。また、リーマンショックで期待値を下回る様子が見て取れる。
　JBのp値は0.250と正規分布でないとは言い切れない状態にある。そこで、目視により残差を確認するため、以下を実行する。

# ヒストグラム：残差
result.resid.hist(label='residual')
plt.xlabel('residual')
plt.ylabel('frequency')
plt.legend(loc='upper right')

　ここまで見てきた通り、日経平均株価、国内総生産、労働人口がバブル崩壊を機に上昇から下落に転じた様子が見られた。そこで、それを確かめてみる。分析は、バブル崩壊前と崩壊後の２つの時期に分ける。

バブル崩壊前
# バブル崩壊前
x_before_bubble = x[:'1990/1/1']
y_before_bubble = ln225[:'1990/1/1']
model = sm.OLS(y_before_bubble, x_before_bubble)
result = model.fit()
print(result.summary())

# グラフを書く
f, ax = plt.subplots() # 2軸グラフ
# 1軸目
ax.plot(y_before_bubble, label='N225')
result.fittedvalues.plot(label='fitted', style='o', ax=ax)
ax.plot(x_before_bubble['GDP']-24, label='GDP', linestyle='--') # 調整
ax.set_ylabel('log Nikkei225 Index')
ax.legend(loc='lower right')
# 2軸目
ax2 = ax.twinx()
ax2.plot(x_before_bubble['WorkPop'], label='WorkPop', linestyle=':')
ax2.legend(loc='upper left')

# ヒストグラム：残差
result.resid.hist(label='residual')
plt.xlabel('residual')
plt.ylabel('frequency')
plt.legend(loc='upper right')

　決定係数は0.94とすべての分析期間から大きく改善した。またAIC、BIC、JBともに改善していることもわかる。回帰係数のp値もすべて5％以下に収まっている。

バブル崩壊後

# バブル崩壊後
x_after_bubble = x['1990/1/1':]
y_after_bubble = ln225['1990/1/1':]
model = sm.OLS(y_after_bubble, x_after_bubble)
result = model.fit()
print(result.summary())

\(R^2\)、AIC、BICともに悪化している。回帰係数のp値はすべて30％以上である。これはどの説明変数も被説明変数の動きをとらえていない可能性を示唆している。バブル崩壊前後で経済が大きく変化する可能性があることがわかる。

　そこで説明変数にドル円の為替レートを追加してみる。

# ドル円を追加
x = pd.concat([x, fx], axis=1).dropna()

# バブル崩壊後 + ドル円
x_after_bubble_fx = x['1990/1/1':]
y_after_bubble = ln225['1990/1/1':]
model = sm.OLS(y_after_bubble, x_after_bubble_fx)
result = model.fit()
print(result.summary())

# グラフを書く
f, ax = plt.subplots() # 2軸グラフ
# 1軸目
ax.plot(y_after_bubble, label='N225')
result.fittedvalues.plot(label='fitted', style='x', ax=ax)
ax.plot(x_after_bubble['GDP']-24, label='GDP', linestyle='--') # 調整
ax.set_ylabel('log Nikkei225 Index')
ax.legend(loc='lower right')
# 2軸目
ax2 = ax.twinx()
ax2.plot(x_after_bubble['WorkPop'], label='WorkPop', linestyle=':')
ax2.legend(loc='upper left')

　回帰係数のp値について着目してみると、切片、労働人口、国内総生産は大幅に上昇し、ドル円の為替レートは5％以下に収まっている。決定係数、AIC、BIC、JBは前のモデルより改善している。グラフを見てみると、リーマンショック以降の期待値は実際の動きをよく説明できているように思われる。

バブル崩壊後をより詳しく。。。

　バブル崩壊後の経済についてより詳しく見てみる。期間は、以下のように分けた。

  • 1990/1/1 ～ 2000/1/1
  • 2000/1/1 ～ 2008/1/1
  • 2008/1/1 ～ 2016/8/31
import statsmodels.stats.api as sms
sms.jarque_bera(result.resid)

terms = ['1990/1/1', '2000/1/1', '2008/1/1', end]
columns = [] # 表の列名作成用
R2 = []
AIC = []
BIC = []
p_values = []
JB = []
for t, term in enumerate(terms[:-1]) :
    columns.append(terms[t]+ '～' + terms[t+1])
    x_t = x[terms[t]:terms[t+1]]
    y_t = ln225[terms[t]:terms[t+1]]
    model = sm.OLS(y_t, x_t)
    result = model.fit()
    R2.append(result.rsquared)
    AIC.append(result.aic)
    BIC.append(result.bic)
    p_values.append(result.pvalues.values)
    JB.append(sms.jarque_bera(result.resid)[1])
    
# dfにまとめる
df_0 = pd.DataFrame([R2, AIC, BIC, JB], index=["R2", "AIC", "BIC", "JB"])
df_1 = pd.DataFrame(p_values).T
df_1.index = ['切片', '労働人口', '国内総生産', 'ドル円']
df = pd.concat([df_0, df_1])
df.columns = columns

pd.options.display.float_format = '{:.2f}'.format # 桁数を調整
df

表から各期間で説明変数のp値が大きく異なり被説明変数に与える影響が異なることがわかる。1990~2000では、労働人口とドル円が日経平均株価の有力な説明変数となっている。2000~2008では、国内総生産がそれに加わっている。2008~2016では、為替レートだけが有力な説明変数となった。つまり、期間ごとに適切な説明変数を選択していく必要があるように思われる。

まとめ

　今回は、日経平均株価とマクロ変数との関係について分析してみました。分析を通して、期間ごとで有効なマクロ変数があることがあることがわかりました。マクロ変数を追加することで、モデルを改善できることもわかりました。ただし、今回使用したマクロ変数は３つ程度なので、他のマクロ変数を加えてより詳細に分析をしてみるなど、色々試す余地は残っています。

　分析では、簡単な残差分析、多変量解析を行いました。残差の分析では、Jarque-Bara検定を用いました。また、多変量解析では、OLSにより説明位変数の有効性確認や統計量からモデルの良さを判定しました。更にグラフを書くことでモデルの予測値(期待値)や残差を目視で確認しました。ここまでの内容で最低限の時系列分析はできるようになるのではないかと思います。今後は、これまでの知識を生かした売買ルールの作成など行ってみたいと思います。

以上

タイトルとURLをコピーしました