「回帰分析から学ぶ計量経済学」をPythonで写経 ～ 第５章「時系列分析」④定常系列と非定常系列
第５章「時系列分析」
書籍の著者 山澤成康 先生
この記事は、書籍「回帰分析から学ぶ計量経済学」第５章「時系列分析」の Python写経活動 を取り扱います。
今回は非定常系列にフォーカスします。
見せかけの回帰、単位根検定、１階階差系列などを実践します！
では書籍を開いて回帰分析の旅に出発です🚀
はじめに
書籍「回帰分析から学ぶ計量経済学」のご紹介
このシリーズは書籍「回帰分析から学ぶ計量経済学: Excelで読み解く経済のしくみ」（オーム社、「テキスト」と呼びます）の Python 写経です。
テキストは、2023年11月に発売され、副題「Excelで読み解く経済のしくみ」のとおり、主に Excel を用いて、計量経済学を平易に学べる素晴らしい書籍です。
テキストの「はじめに」に著者の先生が執筆の動機を書かれています。
社会人の統計リテラシーの向上をテーマの１つとした科研費プロジェクトの最終年度で、広く社会人に向けてわかりやすい経済分析の本を書きたかったのです。
私にとって計量経済学は高嶺の花ですが、このテキストでさまざまな回帰分析のアプローチを知ることができました。
また、書籍の Excel 処理を Python に置き換える「寄り道写経」の実践を通じて、回帰分析のお気持ちに少し近づけた感じがいたします。
回帰分析に慣れ親しむのに丁度良いレベル感と内容ですので、これはぜひともブログにしたい！と思って現在に至ります。
計量経済学の色を薄め、データ分析の色を濃いめに書いてまいります！
引用表記
この記事は、出典に記載の書籍に掲載された文章と配布データを引用し、適宜、掲載文章・配布データを改変して書いています。
【出典】
「回帰分析から学ぶ計量経済学: Excelで読み解く経済のしくみ」
第１版第１刷、著者 山澤成康、オーム社
記事中のイラストは、「かわいいフリー素材集いらすとや」さんのイラストをお借りしています。
ありがとうございます！
第５章 時系列分析
この記事は第５章の以下の節を取り扱います。
5.13 定常系列と非定常系列
5.14 見せかけの回帰の実験
5.15 単位根検定
5.16 消費関数の例
記事に用いるデータは、オーム社の書籍紹介サイトからダウンロードできる Excel ファイル内のデータをもとにしてCSVファイルを作成し、data フォルダに格納しています。
第５章で用いるライブラリをインポートします。
### インポート
# 数値計算
import numpy as np
import pandas as pd
# 統計処理
import scipy.stats as stats
import statsmodels.formula.api as smf # フォーミュラ構文
import statsmodels.tsa.api as tsa # 時系列
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf # コレログラム等
import pmdarima as pm # ARIMA次数自動探索
import lmdiag # 残差プロット
# 描画
import matplotlib.pyplot as plt
import seaborn as sns
import japanize_matplotlib
# ワーニング非表示
import warnings
warnings.simplefilter('ignore')
5.13 定常系列と非定常系列
■ ランダムウォーク p.174～
原系列が非定常で差分系列が定常であるような時系列は単位根過程と呼ばれ、ランダムウォークは単位根過程の一種です。
テキストよりランダムウォークの式をお借りします。
$$
X_t = X_{t-1} + u_t
$$
ランダムウォークデータを生成して眺めてみましょう。
### ランダムウォーク
## 設定と準備
T = 100 # 時間（標本サイズ）
rng = np.random.default_rng(seed=0) # 乱数生成器
## ランダムウォークデータの作成
random_walk = np.cumsum(rng.standard_normal(T))
## ランダムウォークの時系列プロットの描画
plt.figure(figsize=(9, 3))
plt.plot(random_walk)
plt.xlabel('時間')
plt.ylabel('値')
plt.title('ランダムウォーク')
plt.grid(lw=0.5)
# コレログラムと偏自己相関プロットの描画
plot_acf_pacf(random_walk, 30)【実行結果】
式に $${ X_{t-1}}$$ が含まれており、１期間前に影響されて値が生成されています。
偏自己相関係数を見るとラグ１の偏自己相関係数が高くなっています。
ランダムウォークの１階階差をとってみましょう。
１期間前との差分を計算します。
### ランダムウォークの1階階差
## 1階階差の算出
random_walk_I1 = np.diff(random_walk, n=1)
## ランダムウォークの時系列プロットの描画
# 描画領域の設定
plt.figure(figsize=(9, 3))
# ランダムウォークの折れ線グラフの描画
plt.plot(random_walk, label='ランダムウォーク')
# ランダムウォークの1階階差の折れ線グラフの描画
plt.plot(random_walk_I1, color='tab:red', label='1階階差')
# y=0の水平線の描画
plt.axhline(0, color='black', ls='--', lw=1)
# 修飾
plt.xlabel('時間')
plt.ylabel('値')
plt.title('ランダムウォークと1階階差')
plt.grid(lw=0.5)
plt.legend()
plt.show()【実行結果】
青い線がランダムウォーク、赤い線がランダムウォークの１階階差系列です。
１階階差をとると、平均値０に回帰する定常系列に見えますね！
■ ドリフト付きランダムウォーク p.175
ランダムウォークに定数項$${\alpha}$$が付いています。
テキストよりドリフト付きランダムウォークの式をお借りします。
$$
X_t = \alpha + X_{t-1} + u_t
$$
ドリフト付きランダムウォークデータを生成して眺めてみましょう。
### ドリフト付きランダムウォーク
## 設定と準備
T = 100 # 時間（標本サイズ）
alpha = 0.2 # 定数項
rng = np.random.default_rng(seed=0) # 乱数生成器
## ドリフト付きランダムウォークデータの作成
drift_random_walk = np.cumsum(alpha + rng.standard_normal(T))
## ドリフト付きランダムウォークの時系列プロットの描画
# 描画領域の設定
plt.figure(figsize=(9, 3))
# ドリフト付きランダムウォークの折れ線グラフの描画
plt.plot(drift_random_walk, label='ドリフト付きランダムウォーク')
# 通常のランダムウォークの折れ線グラフの描画
plt.plot(random_walk, color='tab:red', ls='--', label='ランダムウォーク')
# 定数項の直線の描画
plt.plot(np.cumsum(np.repeat(alpha, T)), color='tab:orange', ls='--',
label=rf'定数項 $\alpha$={alpha}')
# 修飾
plt.xlabel('時間')
plt.ylabel('値')
plt.title('ドリフト付きランダムウォーク')
plt.grid(lw=0.5)
plt.legend()
plt.show()
# コレログラムと偏自己相関プロットの描画
plot_acf_pacf(random_walk, 30)【実行結果】
毎期間、定数項の同じ値が加算されていますので、定数項は直線的な上方トレンドを示すようです（オレンジ色の点線）。
青色のドリフト付きランダムウォーク全体としても、上方トレンドを持っているようです。
ラグ１の偏自己相関係数が高くなっています。
5.14 見せかけの回帰の実験
非定常系列を使って回帰分析を行うと推定結果が実態以上によくなる問題を、見せかけの回帰と呼んでいるようです。
テキストの日経平均株価データとランダムウォークデータの回帰分析を行ってみましょう。
データを読み込みます。
### データの読み込み
# CSVファイルの読み込み
df6 = pd.read_csv('./data/05_05_nikkei.csv', index_col=0, parse_dates=[0])
# データフレームの表示
print('df6.shape:', df6.shape)
display(df6.head())【実行結果】
３０日分の日次日経平均株価データです。
ドリフト付きランダムウォークデータを生成します。
以下の式で生成します。
$$
\begin{align*}
X_t &= 0.5 + X_{t-1} + u_t \\
u_t &\sim \text{Normal}(0, 1)
\end{align*}
$$
### データの前処理 ドリフト付きランダムウォーク系列 X_t + 0.5 + X_t01 + u_tの作成
## 設定と準備
T = df6.shape[0]
rng = np.random.default_rng(seed=2)
## ドリフト付きランダムウォーク系列の作成
X = np.zeros(T)
X[0] = 10
for t in range(T):
X[t] = 0.5 + X[t-1] + rng.standard_normal()
## データフレームに結合
df6 = pd.concat([df6, pd.DataFrame({'X': X}, index=df6.index)], axis=1)
print('df6.shape: ', df6.shape)
display(df6.head())【実行結果】
$${X}$$がドリフト付きランダムウォークです。
２つのデータを可視化します。
### 可視化
# 描画領域の設定
fig, ax = plt.subplots(figsize=(8, 3))
twinx = ax.twinx()
# 日経平均株価の描画
ax.plot(df6['日経平均'], label='日経平均')
# ランダムウォークの描画
twinx.plot(df6['X'], color='tab:red', label='ランダムウォーク X')
# 凡例処理 axとtwinxの凡例を結合して表示
handles1, labels1 = ax.get_legend_handles_labels()
handles2, labels2 = twinx.get_legend_handles_labels()
ax.legend(handles=handles1+handles2, labels=labels1+labels2)
# 修飾
ax.set(xlabel='日付', ylabel='日経平均株価 [円]',
title='日経平均株価とランダムウォーク')
twinx.set_ylabel('値')
ax.grid(lw=0.5)
plt.show()【実行結果】
よく似た傾向を示しています。
いっけん、２つのデータは相関関係などの関連性を持つ印象を受けますが、しかし、２つのデータは一切関係せずに生成されています。
では回帰分析を行います。
### 回帰分析の実行 見せかけの回帰
result = smf.ols(formula='日経平均 ~ X', data=df6).fit()
display(result.summary())【実行結果】
自由度調整済み決定係数$${0.808}$$は高いですね。当てはまりがいい感じ。
また説明変数$${X}$$の係数の$${p}$$値は$${0.000}$$であり、$${5\%}$$水準で有意です。
ただ、日経平均株価データとドリフト付きランダムウォークデータ$${X}$$は何ら関係のないデータであり、意味のある推定にはなっていません。
見せかけの回帰なんです。
5.15 単位根検定
非定常系列かどうかを単位根検定で調べることができます。
テキストはディッキーフラーテストを紹介しています。
帰無仮説は「単位根を持つ」（非定常系列である）です。
テキストによると単位根検定の基本は以下の式を推計して、$${\rho=0}$$の検定をして有意に０と異ならないときに「単位根がある」とするそうです。
$$
X_t - X_{t-1} = \alpha + \rho X_{t-1} + u_t
$$
$${\rho}$$の係数の$${t}$$値は、非定常系列の場合$${t}$$分布にならないことから、別途、臨界値の表を用いて検定するようです。
単位根検定の臨界値の表をテキストより引用いたします。
$$
\begin{array}{c:c:ccc}
& 有意水準 & & 標本サイズ \\
& & 25 & 50 & 100 \\
\hline
ディッキー・フラーテスト & 1\% & -3.75 & -3.59 & -3.50 \\
& 5\% & -2.99 & -2.93 & -2.90 \\
\hline
通常のt検定(片側検定) & 1\% & -2.50 & -2.41 & -2.37 \\
& 5\% & -1.71 & -1.68 & -1.66 \\
\end{array}
$$
次節で単位根検定を実践します。
5.16 消費関数の例
■ ディッキー・フラーテスト p.179
データを読み込みます。
### データの読み込み
# CSVファイルの読み込み
df7 = pd.read_csv('./data/05_06_consumption.csv', index_col=0)
# 四半期日付のインデックスを設定
df7.index = pd.date_range(start='1994-03-31', periods=df7.shape[0], freq='QE')
df7.index.name = '四半期'
# データフレームの表示
print('df7.shape:', df7.shape)
display(df7.head())【実行結果】
117四半期分の実質GDPと実質民間最終消費支出です。
通常のディッキー・フラーテストを実践します。
まずはデータの前処理をします。
### 通常のディッキー・フラーテスト
## データの準備
# データのコピーを作成
df7dft = df7.copy()
# 実質GDPの1階階差とラグ1の列を追加
df7dft['実質GDP1階階差'] = df7dft['実質GDP'].diff(1)
df7dft['実質GDP1期前'] = df7dft['実質GDP'].shift(1)
print('df7dft.shape: ', df7dft.shape)
display(df7dft)【実行結果】
実質GDPの階差と1期前の列を追加しました。
ディッキー・フラーテストのための回帰分析を実行します。
目的変数：実質GDP階差、説明変数：実質GDP1期前です。
## 回帰分析の実行
result = smf.ols(formula='実質GDP1階階差 ~ 実質GDP1期前', data=df7dft).fit()
display(result.summary())【実行結果】
前節の単位根検定の臨界値の表より、標本サイズ$${100}$$で$${5\%}$$水準の$${t}$$値の臨界値は$${-2.90}$$です。
上記の回帰分析の結果、実質GDP1期前の係数の$${t}$$値は$${-1.730}$$であり、臨界値よりも大きいため、帰無仮説を棄却できず、$${\rho=0}$$の可能性あり⇒単位根を持つ可能性があることを示唆しているそうです。
■ １階階差系列どうしの回帰分析の実行 p.179
実質GDP、実質民間最終消費支出の両方の１階階差系列で回帰分析を行います。
データの前処理を行います。
### 1階階差の回帰分析の実行
## データの前処理
# 実質民間最終消費支出の1階階差を追加
df7dft['実質民間最終消費支出1階階差'] = df7dft['実質民間最終消費支出'].diff(1)
print('df7dft.shape: ', df7dft.shape)
display(df7dft)【実行結果】
回帰分析を実行します。
## 回帰分析の実行
result = smf.ols(formula='実質民間最終消費支出1階階差 ~ 実質GDP1階階差',
data=df7dft).fit()
display(result.summary())【実行結果】
単位根検定の臨界値の表より、実質GDP階差の係数の$${t}$$値$${15.258}$$は有意であり、見せかけの回帰でないと言えそうです。
■ （参考）拡張ディッキー・フラー検定
statsmodels の adfuller 関数を用いることで、簡単に拡張ディッキー・フラー検定を行えます。
帰無仮説は「データは単位根を持つ」です。
実質GDPのADF検定を行います。
### ADF検定 実質GDP
result_adf = tsa.adfuller(df7['実質GDP'], maxlag=10, regression='c')
print('統計量:', result_adf[0])
print('p値 :', result_adf[1])
print('次数 :', result_adf[2])【実行結果】
$${p}$$値は$${0.395}$$であり、$${5\%}$$水準で帰無仮説を棄却できず、単位根が無いとは言えないです。
続いて実質GDPの１階階差のADF検定を行います。
### ADF検定 実質GDP 1階階差
result_adf = tsa.adfuller(df7['実質GDP'].diff(1).dropna(), maxlag=10,
regression='c')
print('統計量:', result_adf[0])
print('p値 :', result_adf[1])
print('次数 :', result_adf[2])【実行結果】
$${p}$$値は$${0.000}$$であり、$${5\%}$$水準で帰無仮説は棄却され、単位根が無い（定常である）と言えそうです。
続いて、実質民間最終消費支出とその１階階差のADF検定を行います。
### ADF検定 実質民間最終消費支出
result_adf = tsa.adfuller(df7['実質民間最終消費支出'], maxlag=10, regression='c')
print('統計量:', result_adf[0])
print('p値 :', result_adf[1])
print('次数 :', result_adf[2])【実行結果】
原系列は単位根を持つようです。
### ADF検定 実質民間最終消費支出 1階階差
result_adf = tsa.adfuller(df7['実質民間最終消費支出'].diff(1).dropna(),
maxlag=10, regression='c')
print('統計量:', result_adf[0])
print('p値 :', result_adf[1])
print('次数 :', result_adf[2])【実行結果】
１階階差系列は単位根を持たないと言えそうです。
今回の写経は以上です。
シリーズの記事
次の記事
前の記事
目次
ブログの紹介
note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！
１．のんびり統計
統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。
２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で
書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！
３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で
書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！
４．楽しい写経 ベイズ・Python等
ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀
５．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で
書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。
６．データサイエンスっぽいことを綴る
統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。
７．Python機械学習プログラミング実践記
書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。
最後までお読みいただきまして、ありがとうございました。
いいなと思ったら応援しよう！
応援ありがとうございます。これからもがんばって記事を作成します！