import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
df = pd.read_excel("./gdp.xlsx", sheet_name="Data")
df["time"] = pd.to_datetime(df["年"].astype(str)+"-"+df["月"].astype(str)) # pandas datetimeを追加
plt.plot(df["time"],df["GDP"])
plt.show()上級マクロ経済学
第11回 時系列データと系列相関
1 時系列データ
時変の数量
時間を通じて変化するあるデータの列を、 (X_1, X_2, \cdots, X_T) だけ観測したとする。 このX_t~(t = 0,1,2,\cdots T)の列を、時系列データと呼ぶ。
また、この列を (X_t)_{t=0}^{T} あるいはtの取る範囲を明示する必要がない場合は (X_t)_t などと書く。
時系列データの例
様々なデータが該当する。
pd.to_datetime
- 時系列データのプロットでは、月次、四半期などのデータを表示することになる。
- プロットの横軸が時間だと、例えば[1996年4月,1996年5月,…]とラベルがならんでいると見にくい。
pandasのDatetime型を用いるとラベルを自動でいい感じに調整してくれる。- 基本は、
year-month-dayを文字列で保存した<list>を作り、pd.to_datetime(<list>)とすれば、datetime型にしてくれる。
df = pd.read_excel("./gdp.xlsx", sheet_name="Data")
df["time"] = pd.to_datetime(df["年"].astype(str)+"-"+df["月"].astype(str)) # pandas datetimeを追加
df| 年 | 月 | GDP | time | |
|---|---|---|---|---|
| 0 | 1994 | 1 | 446237.2 | 1994-01-01 |
| 1 | 1994 | 4 | 443814.0 | 1994-04-01 |
| 2 | 1994 | 7 | 448901.3 | 1994-07-01 |
| 3 | 1994 | 10 | 447176.6 | 1994-10-01 |
| 4 | 1995 | 1 | 452038.7 | 1995-01-01 |
| ... | ... | ... | ... | ... |
| 118 | 2023 | 7 | 557205.7 | 2023-07-01 |
| 119 | 2023 | 10 | 557725.0 | 2023-10-01 |
| 120 | 2024 | 1 | 554281.7 | 2024-01-01 |
| 121 | 2024 | 4 | 557253.4 | 2024-04-01 |
| 122 | 2024 | 7 | 558505.5 | 2024-07-01 |
123 rows × 4 columns
景気循環
時系列データに存在するギザギザが、景気循環の存在を暗示。
- トレンド
- 長期的な水準と解釈される。
- 潜在GDPや自然失業率に当たるものと解釈される。
- サイクル
- いくつかの要因で生じた長期的水準との短期的な乖離。
自然失業率
定常状態において生じる失業率
- 定常状態は長期的な均衡と解釈される
モデルでの導出
- 労働力人口をN、失業者をU、労働者をLとする。
- 失業者は次期にfの割合だけ就職する。
- 労働者は次期にsの割合だけ離職する。
このとき、自然失業率u_nは次のように与えられる。 u_n = \frac{s}{s+f}
導出はこちら
定常状態では、離職者と就職者の数が同じになっている。 \begin{aligned} f U &= s L \\ &= s (N-U) \\ f u_n &= s(1-u_n) \hspace{12pt} \left(u_n \equiv \frac{U}{N}\right) \\ u_n &= \frac{s}{s+f} \end{aligned}
潜在GDP（自然算出水準）
事実上利用可能なすべての生産要素を用いた状態での生産水準を
- 自然算出水準
- 完全雇用GDP
- 潜在GDP
などと呼ぶ。
- 生産関数：Y=F(K,L)。
- 自然失業率の時の労働投入量：\bar{L}。
- 利用可能な資本すべて利用したときの資本投入量：\bar{K}。 \bar{Y} = F(\bar{K}, \bar{L})
フィルタリング
潜在GDPや自然失業率を得るためには、得られたデータからトレンド部分・サイクル部分を抽出すればよい。
いくつかの方法があるが、
- HP（Hodrick–Prescott）フィルター
- 移動平均法
などが代表的。
ここでは移動平均法を概説
移動平均法
四半期データであれば、得られたデータ(X_t)_tについて、移動平均を次のように計算する。
移動平均調整済みのデータ\tilde{X}_tとすると、 \tilde{X}_t = \frac{X_{t} + X_{t-1} + X_{t-2} + X_{t-3}}{4}
移動平均調整済みのデータ\tilde{X}_tとすると、 \tilde{X}_t = \frac{X_{t+2} + 2X_{t+1} + 2X_{t} + 2X_{t-1} + X_{t-2}}{8}
- 月次データなら、12ヶ月の平均をとれるように期間を調整する。
トレンド・サイクル分解
df = pd.read_excel("./gdp.xlsx", sheet_name="Data")
"""
中心移動平均
"""
df["time"] = pd.to_datetime(df["年"].astype(str)+"-"+df["月"].astype(str)) # 時間インデックスを作成
df["tGDP"] = np.nan # トレンドを記録する列をNaNで作成
for i in df.index[2:-2]:
val = 0 # 分子
for j in range(-2,3): # 前後2個ずつのラインを取る。
if np.abs(j) >1:
val = val + df["GDP"].iloc[i+j]
else:
val = val + 2*df["GDP"].iloc[i+j]
df.iloc[i,-1] = val/8
plt.plot(df["time"], df["tGDP"], label="MA") # NaNは無視してプロットしてくれる
plt.plot(df["time"], df["GDP"], label="original", ls ="--")
plt.legend()
plt.show()df = pd.read_excel("./unemp.xlsx", sheet_name="Data")
"""
中心移動平均
"""
df["time"] = pd.to_datetime(df["年"].astype(str)+"-"+df["月"].astype(str)) # 時間インデックスを作成
df["t失業率"] = np.nan # トレンドを記録する列をNaNで作成
for i in df.index[6:-6]:
val = 0 # 分子
for j in range(-6,7): # 前後2個ずつのラインを取る。
if np.abs(j) >5:
val = val + df["失業率"].iloc[i+j]
else:
val = val + 2*df["失業率"].iloc[i+j]
df.iloc[i,-1] = val/24
plt.plot(df["time"], df["t失業率"], label="MA") # NaNは無視してプロットしてくれる
plt.plot(df["time"], df["失業率"], label="original", ls ="--")
plt.legend()
plt.show()GDPギャップ
GDPの実際の値とトレンドの値の差を、GDPギャップと呼ぶ。
- これがサイクル部分。
- 景気拡大期には正、減退期には負
- 内閣府などが計算して公開しているので、普通はそちらを利用する。
- 傾向は似ているが値は絶対値はかなり異なる。
df = pd.read_excel("./gdp.xlsx", sheet_name="Data")
"""
中心移動平均
"""
df["time"] = pd.to_datetime(df["年"].astype(str)+"-"+df["月"].astype(str)) # 時間インデックスを作成
df["tGDP"] = np.nan # トレンドを記録する列をNaNで作成
for i in df.index[2:-2]:
val = 0 # 分子
for j in range(-2,3): # 前後2個ずつのラインを取る。
if np.abs(j) >1:
val = val + df["GDP"].iloc[i+j]
else:
val = val + 2*df["GDP"].iloc[i+j]
df.iloc[i,-1] = val/8
df_gap = pd.DataFrame() # GDPギャップ用のデータフレームを準備
df_gap["time"] = df["time"]
df_gap["GDP_gap"] = df["GDP"].values - df["tGDP"].values # 単に差を取る
df_gap["GDP_gap"] = df_gap["GDP_gap"]/df["tGDP"].values*100 # トレンドで基準化して、%表示
df_gov = pd.read_excel("./gdp_gap.xlsx",sheet_name="Data")
df_gov["time"] = pd.to_datetime(df_gov["年"].astype(str)+"-"+df_gov["月"].astype(str))
plt.plot(df_gap["time"], np.zeros(len(df_gap["time"])), alpha=0.5, color="gray") # 0のラインを引く
plt.plot(df_gap["time"], df_gap["GDP_gap"], label="gap from MA") # NaNは無視してプロットしてくれる
plt.plot(df_gov["time"], df_gov["GDP_gap"], label="goverment estimation")
plt.legend()
plt.show()季節調整について
- 今回のトレンド・サイクル分解は季節調整済みのデータから算出している。
- 原系列と呼ばれる生のデータから季節調整を行う場合も、同じようにフィルタリングを行う。
- ただし、その方法はかなり多岐かつ複雑。
- 最も基本的な考え方は、移動平均。
2 系列相関
データの変動を捉える
- 時系列データの変動を説明することはできるか？
- キーは、自己相関。
GDPギャップを例としよう。
- 今季のGDPギャップが、前期のGDPギャップと相関を持つ場合、次期のGDPギャップも今期のGDPギャップと相関を持つはず。
- すると、相関を元に、次期の値を予測することも正当化される。
これ以降の話は定常性を仮定して続ける。
- 定常性はまた次回に言及。
相関係数
確率変数XとYの相関係数は、次のように与えられる。 \rho = \frac{Cov(X,Y)}{\sqrt{Var(X)}\sqrt{Var(Y)}}
データから計算するならば、(X_i,Y_i)_iの観測した値に対して、 \widehat{Cov}(X,Y)= \frac{\sum_{i} (X_i - \bar{X})(Y_i - \bar{Y})}{n-1} \widehat{Var}(X) = \frac{\sum_{i} (X_i - \bar{X})^2}{n-1} で計算するため、 \hat{\rho} = \frac{\sum_{i} (X_i - \bar{X})(Y_i - \bar{Y})}{\sqrt{\sum_{i} (X_i - \bar{X})^2} \sqrt{\sum_{i} (Y_i - \bar{Y})^2}}
自己相関
現在のy_tに対して、過去の変数y_{t-j}との相関をj次の自己相関という。
確率変数の列(y_t)_tを考えよう。 j次の自己相関を\rho_jとして、次のように定義しよう。 \hat{\rho}_j \equiv \frac{\sum_{t=j+1}^{T} (y_t - \bar{y})(y_{t-j} - \bar{y})}{\sum_{t=１}^{^T} (y_{t} - \bar{y})^２} 分母をT-jで割った \widehat{Cov}(y_t,y_{t-j}) = \frac{\sum_{t=j+1}^{T} (y_t - \bar{y})(y_{t-j} - \bar{y})}{T-j} は、自己共分散と呼ばれる。j=1のケース
自己相関が1次のケースなら、次のようになる。 \hat{\rho}_1 \equiv \frac{(\sum_{t=2}^{T} y_t - \bar{y})(y_{t-1} - \bar{y})}{\sum_{t=１}^{^T} (y_{t} - \bar{y})＾２} 1期ずらして相関をとっているだけ。
コレログラム
(\hat{\rho}_j)_jを並べて図示したものを、コレログラムという。 GDPギャップについて、実際に生成してみよう。
"""
データ読み込みと加工
"""
df = pd.read_excel("./gdp_gap.xlsx", sheet_name="Data") # 内閣府提供データ
df["time"] = pd.to_datetime(df["年"].astype(str) + "-" + df["月"].astype(str)) # 時間のインデックスを作成
print(df.head()) # 中身を確認
## 1~20次自己相関を計算
res = [] # 結果保存のリストを作成
for i in range(1,21):
ybar = np.mean(df["GDP_gap"].values) # GDPギャップのサンプル平均値
cov = np.sum((df["GDP_gap"].iloc[i:].values - ybar) * (df["GDP_gap"].iloc[:-i].values - ybar)) #分子（共分散）
var = np.sum((df["GDP_gap"].values-ybar)**2) #分母
cor = cov/var
res.append(cor) # 結果を保存する。
plt.plot(np.arange(1,21), res, marker = "o") 年 月 GDP_gap time
0 1980 1 -0.2 1980-01-01
1 1980 4 -1.6 1980-04-01
2 1980 7 -0.6 1980-07-01
3 1980 10 0.3 1980-10-01
4 1981 1 0.2 1981-01-01時系列解析に特化したstatsmodels.tsa.apiをimportして、関数acf()を利用
import statsmodels.tsa.api as ts # stsatsmodelsの時系列解析モジュールをインポート
"""
データ読み込みと加工
"""
df = pd.read_excel("./gdp_gap.xlsx", sheet_name="Data") # 内閣府提供データ
df["time"] = pd.to_datetime(df["年"].astype(str) + "-" + df["月"].astype(str)) # 時間のインデックスを作成
print(df.head()) # 中身を確認
## 1~20次自己相関を計算
res = ts.acf(df["GDP_gap"].values, nlags=20) # 0次から20次までの自己相関
plt.plot(np.arange(1,21), res[1:], marker = "o") # 1次以降の自己相関をプロット 年 月 GDP_gap time
0 1980 1 -0.2 1980-01-01
1 1980 4 -1.6 1980-04-01
2 1980 7 -0.6 1980-07-01
3 1980 10 0.3 1980-10-01
4 1981 1 0.2 1981-01-01系列相関のジャッジ
y_tと、少なくともひとつのy_{t-j}~ (t=1,2,\cdots)が自己相関を持っていれば、(y_t)_tには系列相関があるという。
- 1次自己相関係数がゼロであるという帰無仮説のもと、統計検定。
- 帰無仮説が棄却されれば系列相関あり。
- 帰無仮説が棄却されなければ次へ。
- 2次自己相関係数がゼロであるという帰無仮説のもと、統計検定。
- これを繰り返していく。
- 統計検定をどうすればよいかが分かればOK。
- 標準誤差が計算できるので、それを利用。
標準誤差と帰無仮説
- j次以前の自己相関がゼロであるという仮定のもとで、\hat{\rho}_jの標準誤差\widehat{SE}(\hat{\rho}_j)は、次のように計算されることが知られている。 \widehat{SE}(\hat{\rho}_j) = \sqrt{\frac{1}{T}}
さらに、\rho_j=0の帰無仮説のもとで、 \frac{\hat{\rho}_j}{\widehat{SE}(\hat{\rho}_j)} \sim \mathcal{N}(0,1) となることが知られている。
区間推定
帰無仮説に基づくならば、95%の信頼度で -1.96 \le \frac{\hat{\rho}_j}{\widehat{SE}(\hat{\rho}_j)} \le 1.96 となっているはず。 同様に、\hat{\rho}_jは95%の信頼度で、 -1.96\times \widehat{SE}(\hat{\rho}_j) \le \hat{\rho}_j \le 1.96 \times \widehat{SE}(\hat{\rho}_j) を満たしているはず。
- \hat{\rho}_jがこれを満たさない場合、帰無仮説\hat{\rho}_j=0が棄却される。
コレログラムでのチェック
import statsmodels.tsa.api as ts # stsatsmodelsの時系列解析モジュールをインポート
"""
データ読み込みと加工
"""
df = pd.read_excel("./gdp_gap.xlsx", sheet_name="Data") # 内閣府提供データ
df["time"] = pd.to_datetime(df["年"].astype(str) + "-" + df["月"].astype(str)) # 時間のインデックスを作成
## 1~20次自己相関を計算
res = ts.acf(df["GDP_gap"].values, nlags=20) # 0次から20次までの自己相関
## SEの計算
T = len(df) # サンプル期間
SE = (1/T)**(1/2) # SEの計算
## 区間の生成
top = np.ones(20)*SE*1.96 # [1.96, 1.96, ...]を作成
bottom = np.ones(20)*SE*(-1.96) # [-1.96, -1.96, ...]を作成
plt.plot(np.arange(1,21), res[1:], marker = "o", label="correlation") # 1次以降の自己相関をプロット
plt.fill_between(np.arange(1,21), y1=top, y2=bottom, color="gray", alpha = 0.3, label=r"$\pm1.96\times$SE") # シャドー部分を作成
plt.legend() # 凡例を作成
plt.show()ARモデルの導入
1次の自己相関があれば、次のモデルが定式化できる。
y_t = \phi_0 + \phi_1 y_{t-1} + u_t
- このモデルを、自己回帰モデル（AutoregRessive model: AR model）と呼ぶ。
- u_tはt期に発生したショックと呼ばれる。
- E(u_t)=0およびE(u_t^2)=\sigma^2（\sigmaは定数）、さらにu_tは独立かつ同一の分布に従う（iid）と仮定されている。
- このようなショックを、ホワイトノイズと呼ぶ。
ホワイトノイズ
import scipy.stats as stat
x = np.arange(1000)
y0 = stat.uniform.rvs(size=1000, loc=-5, scale=10) # [loc, loc+scale]区間の均一分布からsizeだけランダムに生成
y1 = stat.norm.rvs(size=1000, loc=0, scale=1) # 平均loc、標準偏差scaleの正規分布から、sizeだけランダムに生成
plt.plot(x,y0, ls="-", label="uniform dist.", alpha=0.5)
plt.plot(x,y1,ls="--", label="normal dist.",alpha=0.8)
plt.legend()
plt.show()まとめ
- かなりのデータが時系列データとして解釈可能。
- 時系列データは系列相関を持つことが多い。
- コレログラムでチェック。
- 系列相関があるならば、過去の値から未来の値を予測できるはず。
- 次回からARモデルを概説します。