# ARモデルのラグ次数の選択
def get_aic_bic(gap, y, y_hat, param1, param2):
sse = ((y - y_hat) ** 2).sum(axis=0)
aic = math.log(sse / (len(gap) - param1)) + (l + 1) * 2 / (len(gap) - param2)
bic = math.log(sse / (len(gap) - param1)) + (l + 1) * math.log(
(len(gap) - param2)
) / (len(gap) - param2)
return aic, bic
MAX_L = 4
aic_result, bic_result = np.zeros([MAX_L + 1, 3]), np.zeros([MAX_L + 1, 3])
# AR(0)モデルの場合
l = 0
# AIC・BICの計算方法1:推定するすべてのARモデルのサンプルサイズが同じになるよう調整した場合
y = pd.DataFrame(data[MAX_L:])
aic_result[l, 0], bic_result[l, 0] = get_aic_bic(data, y, np.mean(y.to_numpy()),
MAX_L, MAX_L)
# AIC・BICの計算方法2:推定するすべてのARモデルのサンプルサイズが同じになるよう調整した場合（標本分散の自由度を調整）
y = pd.DataFrame(data[MAX_L:])
aic_result[l, 1], bic_result[l, 1] = get_aic_bic(
data, y, np.mean(y.to_numpy()), MAX_L + l + 1, MAX_L
)
# AIC・BICの計算方法3:それぞれのARモデルで利用できる観測値をすべて利用した場合
y = pd.DataFrame(data[l:])
aic_result[l, 2], bic_result[l, 2] = get_aic_bic(data, y, np.mean(y.to_numpy()),
l, l)
# AR(p)モデル（p>0）の場合
for l in range(1, MAX_L + 1):
# AIC・BICの計算方法1:推定するすべてのARモデルのサンプルサイズが同じになるよう調整した場合
x, y = pd.DataFrame(), pd.DataFrame(data[MAX_L:])
for ll in range(1, l + 1):
x_lagged = pd.DataFrame(data[MAX_L - ll : MAX_L - ll + len(y)].to_numpy())
x = pd.concat([x, x_lagged], axis=1)
model = LinearRegression()
model.fit(x, y)
aic_result[l, 0], bic_result[l, 0] = get_aic_bic(data, y, model.predict(x),
MAX_L, MAX_L)
# AIC・BICの計算方法2:推定するすべてのARモデルのサンプルサイズが同じになるよう調整した場合（誤差項の標本分散の自由度を調整）
x, y = pd.DataFrame(), pd.DataFrame(data[MAX_L:])
for ll in range(1, l + 1):
x_lagged = pd.DataFrame(data[MAX_L - ll : MAX_L - ll + len(y)].to_numpy())
x = pd.concat([x, x_lagged], axis=1)
model = LinearRegression()
model.fit(x, y)
aic_result[l, 1], bic_result[l, 1] = get_aic_bic(
data, y, model.predict(x), MAX_L + l + 1, MAX_L
)
# AIC・BICの計算方法3:それぞれのARモデルで利用できる観測値をすべて利用した場合
x, y = pd.DataFrame(), pd.DataFrame(data[l:])
for ll in range(1, l + 1):
x_lagged = pd.DataFrame(data[l - ll : l - ll + len(y)].to_numpy())
x = pd.concat([x, x_lagged], axis=1)
model = LinearRegression()
model.fit(x, y)
aic_result[l, 2], bic_result[l, 2] = get_aic_bic(data, y, model.predict(x),
l, l)
# 結果表の作成
df_aic_result = pd.DataFrame(
aic_result,
columns=[f'方法{i+1}' for i in range(3)],
index=[f'L={i}' for i in range(MAX_L + 1)],
)
df_bic_result = pd.DataFrame(
bic_result,
columns=[f'方法{i+1}' for i in range(3)],
index=[f'L={i}' for i in range(MAX_L + 1)],
)
print('[AIC]\n', df_aic_result)
print('[BIC]\n', df_bic_result)
print('[best AIC lag length]\n', df_aic_result.idxmin())
print('[best BIC lag length]\n', df_bic_result.idxmin())
