Pythonで挑む時系列分析：supF検定で見抜く経済構造変化
経済や金融の世界では，時として劇的な変化が起こる．オイルショックやバブル崩壊，金融危機など，これらの出来事は経済構造を根本から変えてしまうことがある．しかし，そのような変化点を特定することは非常に難しい課題である．
ここで登場するのがsupF検定という強力な統計的手法だ．この検定は時系列データの中に潜む構造変化を科学的に検出する方法として評価されている．
本記事では，このsupF検定をPythonで実装する方法を，ステップバイステップで解説する．単なる理論の説明にとどまらず，実際の日本のGDPデータを用いて，バブル崩壊後の経済構造の変化を検出する実践的な例も紹介する．
遅ればせながら，自己紹介を行う．
私の名は穐谷といい，現在，慶應義塾大学商学部で計量経済学を専攻している．
興味のある方はこちらを読んでもらえると嬉しい．
では，本題に戻ろう．
構造変化とは
時系列データを扱う上で，｢構造変化｣という概念は非常に重要である．
時系列データでは，パラメータに構造変化すなわちパラメータ値の変化が生じる可能性がある．たとえばオイルショックやバブル崩壊は変数間の相互関係を変えるため，パラメータ値の変化を引き起こしうる．仮に構造変化があったならば，こうしたパラメータ値のヘカを考慮したうえでモデルを推定しないと，誤ったモデルを推定していることになり，推定結果に歪みが生じる．
構造変化点の求め方
まず構造変化点 $${T_B}$$ が既知の場合の求め方を説明する．
構造変化点である $${T_B}$$ 期において，パラメータに構造変化が生じたかどうかを求めたい．
標本期間は $${1}$$ 期から $${T}$$ 期まであり，これを構造変化点である $${T_B}$$ 期の前後で2分割する．
前半期間が $${1}$$ 期から $${T_B}$$ 期まで，後半期間が $${T_B + 1}$$ 期から $${T}$$ 期までとする．
時点 $${t}$$ が前半期間にあれば，
$$
Y_t = \boldsymbol{\beta}\boldsymbol{X_t} + \epsilon_t
\tag{1}
$$
ここで，$${\boldsymbol{\beta} = (\beta_0, \beta_1, \cdots, \beta_K)^T}$$ ，$${\boldsymbol{X_t} = (1, X_{1t}, X_{2t}, \cdots, X_{Kt})^T}$$ とする．
時点 $${t}$$ が後半期間にあれば，
$$
Y_t = \boldsymbol{\beta'}\boldsymbol{X}_t + \epsilon_t \tag{2}
$$
前半期間と後半期間のパラメータ値が構造変化によって異なっていれば，構造変化があるといえる．したがって，仮説を次のように設定する．
$$
H_0: \boldsymbol{\beta} = \boldsymbol{\beta'} \\
H_1: not \ H_0
$$
帰無仮説はパラメータには構造変化ない事を示し，対立仮説は少なくとも1つのパラメータに構造変化がある事を示す．
F検定を行うために，ダミー変数を定義する．ダミー変数 $${D_t}$$ は時点 $${t}$$ が前半期間ならば 0 をとり，後半期間ならば 1 をとるとする．ダミー変数を用いることで式(1)と(2)の2本の式を一本にまとめることが出来る．
$$
Y_t = \boldsymbol{\beta}\boldsymbol{X}_t + \boldsymbol\theta {D_t}\boldsymbol{X}_t + \epsilon_t \tag{3}
$$
ここで，$${\boldsymbol\theta = \boldsymbol\beta' - \boldsymbol\beta}$$ と定義する．
$${\boldsymbol\theta}$$ を用いると仮説は次のように書き換えられる
$$
H_0: \boldsymbol{\theta} = 0 \\ H_1: not \ H_0
$$
式(3)は対立仮説が正しい．つまり構造変化があるとした制約のないモデルである．この式を推定することで，対立仮説が正しいもとでの残差2乗和 $${{SSR}_1}$$ を計算できる．
帰無仮説が正しい．つまり構造変化がないとした制約のあるモデルは式(1)である．この式を推定することで，残差2乗和 $${{SSR}_0}$$ を推定できる．
これらの情報を用いてF統計量を求めることが出来る．
$$
F = \frac{(SSR_0 - SSR_1) / (K+1)} {SSR_1 / (T - 2(K+1))}.
$$
supF検定の詳細
概要
構造変化点は未知である場合が殆どであり，どの時点で生じても不思議ではない．したがって全ての構造変化点の候補でF検定を行ない，得られたF値の中で最大値を取るものをsupF統計量として用いる．
手順
構造変化点の候補 $${T_B}$$ を定めるために，構造変化点の候補の始期と後期を以下のように定義する．
$$
T_{min} = 0.15 \times T \\ T_{max} = (1 - 0.15) \times T
$$
候補点 $${T_B}$$ は以下のようになる
$$
T_B = T_{min}, T_{min} + 1, \cdots, T_{max} - 1, T_{max}.
$$
$${T = 100}$$ の時， $${T_{min} = 15}$$ ， $${T_{max} = 85}$$ であり，構造変化点の候補は$${15～85}$$ となる．
なお最初と最後の $${15%}$$ のデータを除く理由は，構造変化点の候補の前後で十分なサンプルサイズを確保する必要があるためだ．
全ての $${T_B}$$ で $${F}$$ 検定を行ない， $${F}$$ 値を計算する．
$${sup F}$$ 統計量は
$$
\sup F = \max F(T_{\min}), F(T_{\min} + 1), \cdots, F(T_{\max} -1), F(T_{\max})
$$
(なぜかうまく表示されない．)
得られたF値の中で最大のものをsupF統計量とする．
フローチャートを書くと以下のとおりだ．
[開始]
↓
[データの準備]
↓
[構造変化点候補の設定 ($${T_{min} ~ T_{max}}$$)]
↓
[各候補点でF統計量を計算]
↓
[F統計量の最大値 ($${supF}$$) を求める]
↓
[supF統計量と臨界値を比較]
↓
[構造変化の有無を判断]
↓
[終了]
supF検定の実装方法
データは入門 実践する計量経済学から使用する．
期間 $${1981Q1 - 2007Q4}$$ における日本の実質GDP成長率の四半期データを用いる．
なお，季節性を除くため前年同期比とし，100を掛けることで%表示とした．
式は以下のようになる
$$
{growth}t = \beta_0 + \beta_1{{growth}{t-1}} + \beta_2D_t + \beta_3{D_t}{{growth}_{t-1}} + {\epsilon}_t.
$$
$${growth:}$$ 日本の実質GDP成長率
$${D:}$$ 構造変化後は $${1}$$ を取るダミー変数
ライブラリのインポート
# ライブラリをインポート
# データ操作と数値計算のため
import pandas as pd
import numpy as np
# データ可視化のため
import matplotlib.pyplot as plt
# 統計モデリングと計量経済分析のため
import statsmodels.api as sm
# for ignore warning
import warnings
warnings.filterwarnings('ignore')
データのインポート
# データのインポート
url = "<https://www.fbc.keio.ac.jp/~tyabu/keiryo/growth_data.csv>"
_df = pd.read_csv(url)
# dateオブジェクトを作成
date_range = pd.date_range(
start = '1981-03-01',
end = '2007-12-31',
freq = 'Q' # 四半期
)
_df['date'] = date_range # date列に上書き
_df['growth_lag'] = _df['growth'].shift() # 前期のgrowth
df = _df.dropna()
df.head()
日本のGDP成長率の推移をプロットする
plt.figure(figsize=(12, 6))
plt.plot(
df['date'],
df['growth']
)
plt.xlabel('date')
plt.ylabel('growth')
plt.show()
1988年を過ぎた頃からトレンドが変化しているため，1990年前後に構造変化が起こっていそうだ．
構造変化点の候補を出す
# 構造変化点の候補を計算する．
T = len(df) # サンプルサイズ
# 構造変化点の候補の始期
T_min = round(0.15 * T) # 構造変化点の候補の始期．整数を用いる
T_min_date = _df.sort_values('date').iloc[T_min - 1]['date'] # date列でソートし，T_min番目の値を取得
print(f'構造変化点の候補の始期:{T_min_date}')
# 構造変化点の候補の終期
T_max = round((1 - 0.15) * T) # 構造変化点の候補の終期．整数を用いる
T_max_date = _df.sort_values('date').iloc[T_max - 1]['date'] # date列でソートし，T_max番目の値を取得
print(f'構造変化点の候補の終期:{T_max_date}')
# 構造変化点の候補 T_B
T_B_list = list(range(T_min, T_max + 1)) # 構造変化点の候補のリストを作成．構造変化点の候補の始期: 1984-12-31 00:00:00
構造変化点の候補の終期: 2003-09-30 00:00:00
よって，1984年第4四半期～2003年第3四半期の間に構造変化点があるかどうかを探る．
帰無仮説が正しい元でOLS
# 帰無仮説が正しいという前提で推定
endog = df['growth'] # 被説明変数
exog = df['growth_lag'] # 説明変数
exog = sm.add_constant(exog) # 定数項を追加
H_0 = sm.OLS(endog, exog).fit()
print(H_0.summary()) OLS Regression Results
==============================================================================
Dep. Variable: growth R-squared: 0.691
Model: OLS Adj. R-squared: 0.688
Method: Least Squares F-statistic: 234.5
Date: Sun, 11 Aug 2024 Prob (F-statistic): 1.62e-28
Time: 16:55:18 Log-Likelihood: -172.92
No. Observations: 107 AIC: 349.8
Df Residuals: 105 BIC: 355.2
Df Model: 1
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 0.3965 0.181 2.187 0.031 0.037 0.756
growth_lag 0.8304 0.054 15.314 0.000 0.723 0.938
==============================================================================
Omnibus: 9.098 Durbin-Watson: 2.223
Prob(Omnibus): 0.011 Jarque-Bera (JB): 8.881
Skew: 0.650 Prob(JB): 0.0118
Kurtosis: 3.550 Cond. No. 5.37
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.# 残差２乗和を取得
ssr_0 = np.sum(H_0.resid**2)SSR_0 = 158.7212084003789
構造変化点の候補を1つずつ試す
# F統計量を格納するデータフレームを作成
columns = ['date', 'f_stat']
df_f_stat = pd.DataFrame(columns=columns)
# 構造変化点の候補を一つずつ計算する．
for T_B in T_B_list:
T_B_date = _df.sort_values('date').iloc[T_B - 1]['date'] # date列でソートし，T_B番目の値を取得
_df['d'] = (_df['date'] > T_B_date).astype(int) # 構造変化点の候補T_Bの期以前は0，それ以降は1
_df['d_growth'] = _df['d'] * _df['growth'].shift() # dが1のときだけ前期のgrowthを入れる．
df = _df.dropna()
# 対立仮説が正しいという前提で推定
endog = df['growth'] # 被説明変数
exog = df[['growth_lag', 'd', 'd_growth']] # 説明変数
exog = sm.add_constant(exog) # 定数項
H_1 = sm.OLS(endog, exog).fit()
ssr_1 = np.sum(H_1.resid**2) # 残差２乗和を取得
# F統計量を計算する
K = 1 # 制約なしモデルの説明変数の数
f_stat = (ssr_0 - ssr_1) / (K + 1) / ssr_1 * (T - 2 * (K + 1)) # F統計量を計算する．
# 構造変化点の候補点と得られたF統計量をデータフレームに格納する
new_f_stat = pd.Series(
[T_B_date, f_stat],
index = df_f_stat.columns
)
df_f_stat = pd.concat(
[
df_f_stat,
new_f_stat.to_frame().T
],
ignore_index = True
)supF統計量を算出する
sup_f_stat = df_f_stat['f_stat'].max() # F値の最大値
structural_break_point = df_f_stat.loc[df_f_stat['f_stat'].idxmax(), 'date'] # 構造変化点を取得
print(f'構造変化点: {structural_break_point}')
print(f'sup F 統計量: {sup_f_stat:.02f}')構造変化点: 1991-03-31 00:00:00
sup F 統計量: 9.25
sup F 統計量は9.25となり，構造変化点は1991年第1四半期と判明した．
構造変化点候補別のF統計量をプロットする
F値をプロットする．
plt.figure(figsize=(12, 6))
plt.plot(
df_f_stat['date'],
df_f_stat['f_stat'],
)
plt.xlabel('date')
plt.ylabel('F-stat')
current_ymin, current_ymax = plt.ylim()
plt.ylim(current_ymin, current_ymax * 1.2)
# 最大値の位置にテキストで値を表示
plt.annotate(
f'sup F = {sup_f_stat:.2f}',
xy = (structural_break_point, sup_f_stat),
xytext = (10, 10),
textcoords = 'offset points',
arrowprops = dict(arrowstyle='->', connectionstyle='arc3,rad=0')
)
1991年第1四半期を構造変化点とすると $${F}$$ 統計量が最も大きい値を取っていることがわかる．
表によると，制約数が1の時の有意水準5%の臨界値は $${8.68}$$ である．
supF統計量は $${9.25}$$ なので，5%の有意水準で帰無仮説 $${H_0}$$ を棄却し，対立仮説 $${H_1}$$ を採択した．すなわち，日本のGDP成長率には構造変化があり，その構造変化点は1991年第1四半期だという事がわかった．
構造変化点を1991Q1としてOLS推定
OLS Regression Results
==============================================================================
Dep. Variable: growth R-squared: 0.738
Model: OLS Adj. R-squared: 0.730
Method: Least Squares F-statistic: 96.64
Date: Sun, 11 Aug 2024 Prob (F-statistic): 7.98e-30
Time: 16:55:20 Log-Likelihood: -164.08
No. Observations: 107 AIC: 336.2
Df Residuals: 103 BIC: 346.9
Df Model: 3
Covariance Type: nonrobust
=================================================================================
coef std err t P>|t| [0.025 0.975]
---------------------------------------------------------------------------------
const 2.2968 0.583 3.943 0.000 1.141 3.452
growth_lag 0.5041 0.122 4.144 0.000 0.263 0.745
d_1991 -1.9470 0.611 -3.189 0.002 -3.158 -0.736
d_1991_growth 0.1735 0.151 1.149 0.253 -0.126 0.473
==============================================================================
Omnibus: 7.665 Durbin-Watson: 1.986
Prob(Omnibus): 0.022 Jarque-Bera (JB): 7.269
Skew: 0.560 Prob(JB): 0.0264
Kurtosis: 3.612 Cond. No. 27.4
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.推定結果の解釈
1991年第1四半期に構造変化が起こっている事がわかった．1991年は日本のバブル経済が崩壊しはじめた時期である．不動産価格や株価が大きく下落し始め，景気が悪化しはじめた．これは日本経済の構造に大きな変化をもたらした．
constは定数項であり，これは，他の要因が0の場合の基本的な成長率を示す．
growth_lagは，前期のGDP成長率が今季の成長率に与える影響を示す．構造変化前は，前期の成長率が1%上昇すると，今季の成長率は約0.5%上昇することを意味する．
d_1991は1991年第1四半期以降，成長率の基本水準が約1.95%ポイント低下したことを示す．これはバブル崩壊後の経済成長の鈍化を反映している．
d_1991_growthは構造変化後，前期の成長率が今季の成長率に与える影響が約0.17%に変化したことを示している．ただし，統計的に有意ではない．
まとめ
本記事では，時系列データにおける構造変化を検出するための強力なツールであるsupF検定について，その理論的背景から実践的な実装方法まで詳しく解説した．
主な論点は以下の通りである．
構造変化の重要性
supF検定の理論
Pythonによる実装
実データへの応用
結果の解釈
視覚化
この記事を通じて，時系列分析を学ぶ皆さんの役に立てれば嬉しく思う．
コードは汚く，無駄も多いかもしれないが，来年度からシステムエンジニアとして働くので，より良いコードを書けるよう精進したい．
皆さんのスキとフォローをお待ちしています．
参考文献
※アフィリエイトリンク