import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.tsa.api as ts
df = pd.read_excel("gdp_gap.xlsx", sheet_name="Data") # データの読み込み
df["time"] = pd.to_datetime(df["年"].astype(str) + "-" + df["月"].astype(str)) # 時間列の作成
df0 = df[df["time"] < "2016-01-01"] # 2015年中までのデータが既知、2016年以降のデータが未知だとして予測する。
res = ts.AutoReg(df0["GDP_gap"], lags=3).fit(cov_type="HAC", cov_kwds={"maxlags":4}) # AR(3)で、4期ラグのHAC標準誤差を用いる
## 推計パラメータ
phi0 = res.params[0]
phi1 = res.params[1]
phi2 = res.params[2]
phi3 = res.params[3]
sigmahat = (np.mean(res.resid**2))**(1/2) # 誤差の推定値
## GDPギャップデータをappendで扱えるリストに変換
y = df0["GDP_gap"].values.tolist()
N = 100 # シミュレーションの数
h = 34 # maxの予測先
np.random.seed(1234) # シードの固定
# 点推定量
yhat_point = y.copy() # 点推定の結果を入れるリストとして、yのコピーを作成
for j in range(h):
yhat = phi0 + phi1*y[-1] + phi2*y[-2] + phi3*y[-3] # yhatを計算
yhat_point.append(yhat) # 一番後ろに足す。
## 区間予測のためのシミュレーションを行う
sim_res = [] # シミュレーション結果を入れるリストを作成
for i in range(N):
yhat_int = y.copy() # 結果を入れるリストとして、yのコピーを作成
u = np.random.normal(size=h,loc=0, scale=sigmahat) # 平均ゼロ、標準偏差sigmahatの正規分布からランダムに100個の変数をピック
for j in range(h):
yhat = phi0 + phi1*y[-1] + phi2*y[-2] + phi3*y[-3] + u[j] # 確率変動を含むyhatを計算する
yhat_int.append(yhat)
sim_res.append(yhat_int) #シミュレーションした系列を保存
## シミュレーション結果から区間を構築する。
## sim_res[i][j]成分がi回目のトライアルのj期目の値。
sim_res = np.array(sim_res) # j列を取り出しやすいように、ndarrayにして処理
lb = [] # 5%水準の値を保存するリスト
ub = [] # 95%水準の値を保存するリスト
ind_b = int(np.round(N*0.05,0)) # 5%水準のインデックス
ind_u = int(np.round(N*0.95,0)) # 95%水準のインデックス
for j in range(len(y)+h):
vals = sim_res[:,j] # j期目の値のシミュレーション結果を取り出す
lb.append(np.sort(vals)[ind_b]) # j期の下限をリストに加える
ub.append(np.sort(vals)[ind_u]) # j期の上限をリストに加える
## プロットする
plt.rcParams["font.size"]=10
plt.plot(df["time"],np.zeros(len(df)), color="gray", alpha=0.3)
plt.plot(df["time"], df["GDP_gap"], label="data", alpha=0.5) # データのプロット
plt.plot(df["time"].iloc[len(y)-1:len(y)+h], yhat_point[len(y)-1:len(y)+h],ls="--", label="estimate", color="black") # 予測値のプロット
plt.fill_between(df["time"].iloc[len(y)-1:len(y)+h], y1=ub[len(y)-1:len(y)+h], y2=lb[len(y)-1:len(y)+h], color="blue", alpha=0.2, label="90 confidence interval")
plt.legend()
plt.show()