「回帰分析から学ぶ計量経済学」をPythonで写経 ～ 第５章「時系列分析」①時系列データの加工
第５章「時系列分析」
書籍の著者 山澤成康 先生
この記事は、書籍「回帰分析から学ぶ計量経済学」第５章「時系列分析」の Python写経活動 を取り扱います。
今回より時系列データの分析に移ります。
最初の回は時系列データの加工です。
前期比、移動平均などを実践いたします。
では書籍を開いて回帰分析の旅に出発です🚀
はじめに
書籍「回帰分析から学ぶ計量経済学」のご紹介
このシリーズは書籍「回帰分析から学ぶ計量経済学: Excelで読み解く経済のしくみ」（オーム社、「テキスト」と呼びます）の Python 写経です。
テキストは、2023年11月に発売され、副題「Excelで読み解く経済のしくみ」のとおり、主に Excel を用いて、計量経済学を平易に学べる素晴らしい書籍です。
テキストの「はじめに」に著者の先生が執筆の動機を書かれています。
社会人の統計リテラシーの向上をテーマの１つとした科研費プロジェクトの最終年度で、広く社会人に向けてわかりやすい経済分析の本を書きたかったのです。
私にとって計量経済学は高嶺の花ですが、このテキストでさまざまな回帰分析のアプローチを知ることができました。
また、書籍の Excel 処理を Python に置き換える「寄り道写経」の実践を通じて、回帰分析のお気持ちに少し近づけた感じがいたします。
回帰分析に慣れ親しむのに丁度良いレベル感と内容ですので、これはぜひともブログにしたい！と思って現在に至ります。
計量経済学の色を薄め、データ分析の色を濃いめに書いてまいります！
引用表記
この記事は、出典に記載の書籍に掲載された文章と配布データを引用し、適宜、掲載文章・配布データを改変して書いています。
【出典】
「回帰分析から学ぶ計量経済学: Excelで読み解く経済のしくみ」
第１版第１刷、著者 山澤成康、オーム社
記事中のイラストは、「かわいいフリー素材集いらすとや」さんのイラストをお借りしています。
ありがとうございます！
第５章 時系列分析
この記事は第５章の以下の節を取り扱います。
Talk この章のマホナたちの会話
5.3 時系列データの加工（5.1節、5.2節を含みます）
記事に用いるデータは、オーム社の書籍紹介サイトからダウンロードできる 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')
Talk この章のマホナたちの会話
章冒頭の魔女たちの Talk シーンで描画される２つのグラフに取り組みます。
Talk シーンは「章の取り組みテーマのダイジェスト」のような位置づけです。
■ ランダムウォーク p.148
日経平均株価データと人工的に生成したランダムウォークデータを描画します。
テキストによるとランダムウォーク系列は非定常系列であり、「厄介なのは、ランダムウォークどうしで推定すると、必要以上に結果が良くなってしまうこと」、「特に上昇トレンドを持つドリフト付きランダムウォークの決定係数は高くなる」とのこと。
注意が必要そうです。
後段の節で登場する日経平均株価データを転用いたします。
データを読み込みます。
### CSVファイルの読み込み ※p.157のデータを転用
df1 = pd.read_csv('./data/05_01_nikkei.csv', index_col=0, parse_dates=[0])
print('df1.shape: ', df1.shape)
df1.head()【実行結果】
ランダムウォークデータを人工的に生成して、日経平均株価データと並べて可視化します。
ランダムウォークデータは正規分布乱数の累積和（cumsum）で作成しました。
そうなんです、ランダムウォークデータはランダムに生成しているのです。
### 日経平均とランダムウォークの可視化
## 設定と準備
# 乱数生成器の設定
rng = np.random.default_rng(seed=25)
# ランダムウォークデータの生成
diffs = rng.normal(loc=0, scale=1500, size=len(df1))
rw = diffs.cumsum() + 22000
# データフレームにまとめる
data = pd.concat([df1['日経平均株価'],
pd.Series(rw, index=df1.index, name='人工データ')], axis=1)
## 描画
# 描画領域の設定
data.plot(color=['tab:blue', 'tab:red'], figsize=(10, 3))
# 修飾
plt.legend();【実行結果】
何の関係もない２つのデータはとても良く似た傾向を示しています。
相関係数を確認しましょう。
### 相関係数の算出
data.corr().round(3)【実行結果】
日経平均株価と人工データとの間に強い相関関係が見られます。
日経平均株価を目的変数、人工データを説明変数にして、回帰分析を実行してみます。
### 回帰分析の実行
# 回帰分析の実行
result = smf.ols(formula='日経平均株価 ~ 人工データ', data=data).fit()
# 結果の表示
result.summary()【実行結果】
決定係数はそこそこ高めの$${0.644}$$、各種$${p}$$値も有意な値。
「ランダムウォークどうしで推定すると、必要以上に結果が良くなってしまうこと」を体感できました。
ちょっと気になるので、定常性をADF検定で確認してみます。
ADF検定の帰無仮説は「データに単位根が存在し、データは非定常である」です。
日経平均株価データから。
### 日経平均株価のADF検定 ※帰無仮説「データに単位根が存在し、データは非定常である」
result_adf = tsa.adfuller(data['日経平均株価'])
print('統計量:', result_adf[0])
print('p値 :', result_adf[1])
print('次数 :', result_adf[2])【実行結果】
$${p}$$値は$${5\%}$$水準で有意であり、データは非定常とは言えない、となりました。。。。。あれ？
人工データです。
### 人工データのADF検定 ※帰無仮説「データに単位根が存在し、データは非定常である」
result_adf = tsa.adfuller(data['人工データ'])
print('統計量:', result_adf[0])
print('p値 :', result_adf[1])
print('次数 :', result_adf[2])【実行結果】
$${p}$$値は$${5\%}$$水準で有意ではないため帰無仮説を棄却できず、データは非定常であるを受容します。
人工データはひとまず非定常でした、としておきましょう。
■ ARIMA p.151
以下の式で示されるAR(2)モデルから生成されたデータを可視化します。
１期前・２期前の値が説明変数になっています。
$$
X_t = 30 + 0.5 X_{t-1} - 0.8 X_{t-2} + 標準正規乱数 \\
$$
### データの作成と可視化
## 設定と準備
T = np.arange(1, 22) # 時間
x = np.zeros_like(T) # xの配列の初期化
rng = np.random.default_rng(seed=1234) # 乱数生成器
## xの生成 x_t = 30 + 0.5 * x_t-1 - 0.8 * x_t-2 + e
x[0], x[1] = 20, 20
for t in range(2, len(T)):
x[t] = 30 + 0.5 * x[t-1] - 0.8 * x[t-2] + rng.standard_normal()
## 可視化
# 描画領域の設定
plt.figure(figsize=(8, 3))
# 折れ線グラフの描画
plt.plot(T, x)
# 修飾
plt.xlabel('時間')
plt.ylabel('x')
plt.ylim(16, 27)
plt.grid(lw=0.5)
plt.show()【実行結果】
テキストの乱数生成と異なる結果となり、グラフもテキストとは異なっています。
5.3 時系列データの加工
■ 時系列データと時系列分析
時系列データってなんでしょう。
総務省統計局「なるほど統計学園」の時系列データと時系列分析の概要を確認しましょう。
一つの項目について時間に沿って集めたデータを時系列データといい、時間に沿った変化を分析することができます。時系列データを分析する際は季節変動などに注意する必要があります。
時系列データは上記の他に「一定の時間間隔」「時間順に並ぶ」などの特徴を持っています。
■ 仮想の時系列データの作成
5.3節では８つの加工を実践します。
加工元となる仮想の時系列データを以下の式に基づいて作成します。
$$
\begin{align*}
X_0 &= 300 \\
X_{t} &= 50 + 0.8 X_{t-1} + 乱数（0から49の整数）, \quad t=1, \cdots, 71
\end{align*}
$$
### 仮想データの作成
## 設定と準備
T = 72 # 期間
data1 = np.zeros(T) # データを格納する変数の初期化
rng = np.random.default_rng(seed=2) # 乱数生成器
## データの作成
# t=0の初期値
data1[0] = 300
# t=1以降のデータの生成
for t in range(1, T):
data1[t] = 50 + 0.8 * data1[t-1] + rng.integers(0, 50)
# データフレーム化
data1 = pd.DataFrame({'値': data1})
data1.index = pd.period_range('2018-04-01', periods=T, freq='M')
## 結果の表示
print('data1.shape: ', data1.shape)
display(data1.head())
display(data1.tail())【実行結果】
2018年4月から2024年3月までを想定した72か月の仮想データです。
時系列折れ線グラフで月別の推移を可視化しましょう。
### データの描画
data1.plot()
plt.ylabel('value');【実行結果】
大きな峰が３つあります。
■ 【加工１】前期比 p.154
前期比は１期前からの増減比率です。
テキストの数式をお借りします。
$$
\cfrac{X_t - X_{t-1}}{X_{t-1}} \times 100 = \cfrac{X_t}{X_{t-1}} \times 100 - 100
$$
前期比は pandas の pct_change() メソッドで簡単に計算できます。
### 前期比
data1_pct_change = data1.pct_change() * 100
data1_pct_change.head()【実行結果】
2018年4月の１期前：2018年3月データは存在しないため、１行目が欠損となります。
$${100}$$を乗じており、単位は$${\%}$$です。
可視化します。
### 前期比の可視化
data1_pct_change.plot()
plt.ylabel('前期比 %');【実行結果】
おおよそ$${0\%}$$を堺にして上下動している感じです。
定常性がありそうな感じです。
■ 【加工２】前期比年率 p.154
前期比年率は、前期比（月や四半期ごと）の伸び率が１年間継続すると仮定したときの年率です。
四半期の場合の計算式をテキストからお借りします。
$$
\left(\left(\cfrac{X_t}{X_{t-1}}\right)^4 - 1\right) \times 100
$$
月別の仮想データから四半期合計を算出して四半期データにします。
pandas の resample メソッドを利用して、３月を第四四半期末月（rule='Q-MAR'）にした四半期集計を行います。
### 四半期合計の算出
data1_q = data1.resample(rule='Q-MAR', closed='left').sum()
data1_q.head()【実行結果】
３月の年が適用されるため、例えば 2018年4月～2019年3月の各四半期名称には「2019」が付きます（直感的に１年早い？）。
四半期合計を可視化します。
### 四半期合計の可視化
data1_q.plot()
plt.ylabel('value [四半期]');【実行結果】
月別の推移グラフ（２つ前のグラフ）と比べて、小さな上下動が無くなり、直線的な上下動になりました。
大きな峰が３つある特徴は残っています。
では前期比年率を算出します。
１四半期前（１期前）は pandas の shift() メソッドで引数に$${1}$$を与える（または何も指定しない）ことで算出できます。
### 前期比年率の算出
data1_q_annual_rate = ((data1_q / data1_q.shift(1))**4 - 1) * 100
data1_q_annual_rate.head()【実行結果】
前期比年率を可視化しましょう。
### 前期比年率の可視化
data1_q_annual_rate.plot()
plt.ylabel('前期比年率 %');【実行結果】
月次の前期比伸び率のグラフ（２つ前のグラフ）と比べて、かなり形状（もちろん率も）が変わった感じがします。
■ 【加工３】前年同期比 p.155
１年前の同じ期（同じ月、同じ四半期など）と比べた伸び率です。
月次の前年同期比の数式をテキストからお借りします。
$$
\cfrac{X_t - X_{t-12}}{X_{t-12}} \times 100 = \cfrac{X_t}{X_{t-12}} \times 100 - 100
$$
前年同期比を算出します。
月次の仮想データに対して、前期比を計算できる pct_change() メソッドを利用し、pct_change(periods=12) のように期間を12か月＝１年を設定することで、前年同期比を算出できます。
### 前年同期比の算出
data1_pct_change12 = data1.pct_change(periods=12).dropna() * 100
data1_pct_change12.head()【実行結果】
最初の12ヶ月は前年同期がないため、このデータから除外しています。
では前年同期比を可視化しましょう。
### 前年同期比の可視化
data1_pct_change12.plot()
plt.ylabel('前年同期比 %');【実行結果】
■ 【加工４】年平均成長率 p.155
例えば５年前・後の伸び率を１年間に平均するような感じです。
５年間の場合の数式をテキストよりお借りします。
$$
\left(\left(\cfrac{X_{t}}{X_{t-5}}\right)^{1/5} - 1\right) \times 100
$$
月別の仮想データから年度合計を算出して年度データにします。
pandas の resample メソッドを利用して、３月を年度末（rule='Y-MAR'）にした年集計を行います。
### 年度合計の算出
data1_annual = data1.resample(rule='Y-MAR').sum()
data1_annual【実行結果】
３月の年が適用されるため、例えば 2018年4月～2019年3月は「2019」になります（直感的に１年早い？）。
2019年から2024年の年平均成長率を算出します。
### 2019年から2024年の年平均成長率
((data1_annual.loc['2024'] / data1_annual.loc['2019'])**(1/5) - 1) * 100【実行結果】
５年間の年平均成長率は $$2.88\%$$です。
■ 【加工５】寄与度 p.156
全体の伸び率に対して構成要素ごとの寄与の度合いを示す値です。
テキストのGDPが内需と外需で構成される場合の内需・外需の寄与度の数式をお借りします。
$$
\begin{cases}
内需の寄与度 = \cfrac{D_t - D_{t-1}}{GDP_{t-1}} \times 100 \\
外需の寄与度 = \cfrac{E_t - E_{t-1}}{GDP_{t-1}} \times 100 \\
\end{cases}
$$
寄与度の計算を試してみます。
年度別に集約したデータをGDPに見立てて、GDP一部を内需、残りを外需に割り当てます。
### 仮想的な年度別GDP・内需・外需データの作成
ratios = np.array([0.50, 0.49, 0.51, 0.53, 0.55, 0.48])
data1_annual_contri = data1_annual.copy()
data1_annual_contri.columns = ['GDP']
data1_annual_contri['内需'] = data1_annual_contri['GDP'] * ratios
data1_annual_contri['外需'] = data1_annual_contri['GDP'] * (1 - ratios)
data1_annual_contri【実行結果】
内需の寄与度（単位：％）を算出します。
### 内需の寄与度の算出
(
((data1_annual_contri['内需'] - data1_annual_contri['内需'].shift(1))
/ data1_annual_contri['GDP'].shift(1) * 100)
.rename('内需の寄与度').to_frame()
)【実行結果】
外需の寄与度（単位：％）を算出します。
### 外需の寄与度の算出
(
((data1_annual_contri['外需'] - data1_annual_contri['外需'].shift(1))
/ data1_annual_contri['GDP'].shift(1) * 100)
.rename('外需の寄与度').to_frame()
)【実行結果】
GDPの伸び率を算出します。
### GDP伸び率の算出
(
((data1_annual_contri['GDP'] - data1_annual_contri['GDP'].shift(1))
/ data1_annual_contri['GDP'].shift(1) * 100)
.rename('GDP伸び率').to_frame()
)【実行結果】
内需の寄与度と外需の寄与度の和はGDP伸び率（単位：％）になります。
■ 【加工６】移動平均 p.157
移動平均はある時点の前後のデータを平均する手法です。
テキストの日経平均株価データ（冒頭部でも利用したデータ）を読み込み、「３か月中心移動平均」と「３か月後方移動平均」を算出して可視化します。
データを読み込みます。
### データの読み込み
# CSVファイルの読み込み
df1 = pd.read_csv('./data/05_01_nikkei.csv', index_col=0, parse_dates=[0])
# データフレームの表示
print('df1.shape:', df1.shape)
display(df1.head())【実行結果】
３か月中心移動平均と３か月後方移動平均の計算式をテキストよりお借りします。
$$
３か月中心移動平均 = \cfrac{X_{t-1} + X_t + X_{t+1}}{3} \\
３か月後方移動平均 = \cfrac{X_{t-2} + X_{t-1} + X_t}{3} \\
$$
では各移動平均を算出します。
pandas の rolling() メソッドで期間を「移動」させながら処理を行えます。
rolling() メソッドの引数 windowに期間を設定します。
今回は$${3}$$です。
.mean() とすることで、移動させながら３期間の「平均」を計算できます。
### 3か月中心移動平均と3か月後方移動平均の算出
# 3か月中心移動平均の算出 中心計算center=True, 最小期間min_period
df1['3か月中心移動平均'] = df1['日経平均株価'].rolling(window=3, center=True).mean()
# 3か月後方移動平均の算出 後方計算center=False(デフォルト値)
df1['3か月後方移動平均'] = df1['日経平均株価'].rolling(window=3).mean()
# 結果の表示
print('df1.shape:', df1.shape)
display(df1.head())【実行結果】
３か月中心移動平均は前後１か月のデータがあれば算出できるので、最初の欠損は１期間のみです。
３か月後方移動平均は前２か月のデータが必要なので、最初の欠損は２期間となります。
テキストの図表「日経平均の移動平均」にならって可視化します。
### 可視化
plt.figure(figsize=(10, 4))
sns.lineplot(data=df1)
plt.grid(lw=0.5)【実行結果】
移動平均をとると線が滑らかになっています。
後方移動平均の方が上下動が遅れて反映されるようです。
■ 【加工７】階差・対数階差 p.158
テキストは前期との差（１階階差）を紹介しています。
テキストの数式をお借りします。
$$
\Delta X_t = X_t - X_{t-1} \
$$
１階階差を算出します。
pandas の diff() メソッドで簡単に階差を計算できます。
### １階階差（前期との差）
data1_diff1 = data1.diff(1)
data1_diff1.head()【実行結果】
１階階差系列を可視化します。
### １階階差の可視化
data1_diff1.plot()
plt.ylabel('value [1階階差]');【実行結果】
対数階差をとってみます。
元のデータを対数化し、１階階差をとります。
### 対数階差 ※対数階差は前期比伸び率の近似値になる
data1_log_diff1 = data1.apply(np.log).diff(1)
data1_log_diff1.head()【実行結果】
可視化します。
### 対数階差の可視化
data1_log_diff1.plot()
plt.ylabel('value [対数1階階差]');【実行結果】
テキストによると「対数階差は前期比伸び率の近似値になります」とのこと。
前期比伸び率のグラフと比べてみましょう。
確かに、対数階差のグラフと前期比伸び率のグラフはよく似てます！
■ 【加工８】自己相関係数・偏自己相関係数 p.158
$${X_t}$$と$${X_{t-1}}$$の相関係数、$${X_{t-1}}$$と$${X_{t-2}}$$の相関係数のように、同じ変数の異なる時点間の相関が自己相関係数です（ざっくり）。
$${X_t}$$と$${X_{t-2}}$$の自己相関には$${X_{t-1}}$$を介した相関がありえるので、このような影響を取り除いた偏自己相関係数があります。
月次の仮想データの自己相関係数、偏自己相関係数を可視化しましょう。
仮想データは前期の値$${X_{t-1}}$$に$${0.8}$$を乗じており、ラグ１の自己相関係数$${0.8}$$をイメージして作成しました。
statsmodels の plot_acf （自己相関係数）と plot_pacf （偏自己相関係数）を利用します。
### 自己相関係数と偏自己相関係数の描画
# 描画領域の設定
fig, ax = plt.subplots(1, 2, figsize=(10, 3))
# 自己相関係数（コレログラム）の描画
plot_acf(data1['値'], title='自己相関係数', ax=ax[0])
# 偏自己相関係数の描画
plot_pacf(data1['値'], title='偏自己相関係数', ax=ax[1])
# 修飾
ax[0].set_ylim((-0.6, 1.1))
ax[1].set_ylim((-0.6, 1.1));【実行結果】
右の偏自己相関係数を見ると、ラグ１（左から２本目）に偏自己相関係数$${0.8}$$が見られ、その他のラグには有意な偏自己相関係数が無い（青色が有意で無い領域）ことを示しています。
左の自己相関係数（コレログラム）では、自己相関係数が徐々に減衰する様子が見られます。
今回の写経は以上です。
シリーズの記事
次の記事
前の記事
目次
ブログの紹介
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の教科書です。
よかったらぜひ、お試しくださいませ。
最後までお読みいただきまして、ありがとうございました。