§ Introduction

せっかく年号が変わって令和になったので、これまでの年号が続いた年数について時系列分析して、年号の続きやすさについて傾向を調べてみることにしました。

§1. 準備

ひとまず年号が続いた年数のデータを準備してきました。以下からcsvファイルでダウンロードできます。

・nengo : 年号
・start : 年号が始まった年
・years : 年号が続いた年数

室町時代に北朝と南朝に分かれていた時期があったので、今回は北朝に統一してデータを取得することにしました。データを読み込んで、記述統計を行うことにします。

# libraryの読み込み
library(readr)

# 年号データ
dat_nengo <- read_csv("nengo.csv")
dat_nengo

§2 記述統計

原系列と差分系列をプロットして系列の傾向をチェックします。

n_years <- as.ts(dat_nengo$years)
plot(n_years)    # 時系列のプロット
plot(diff(n_years))    # 差分系列のプロット

今回は差分系列に注目しました。この系列を見るといくつかの変化点があり、変化点で分割した各期間ごとでいえば定常な時系列としてモデリングできそうに見えます。

§2. 変化点を調べる

差分系列の変化点を把握することにしました。この程度であれば目で見てもよかったのですが、今回はlibrary(changepoint)を用いました。

# 必要なライブラリの読み込み
library(changepoint)

# 差分系列のオブジェクトを作る。
diff_nyears <- diff(n_years)

# 分散が変化すると仮定して、変化点を検出する。
res_cpts <- cpt.var(diff_nyears, method = "PELT")
param.est(res_cpts)    # 推定されたパラメータの確認
plot(res_cpts)    # 変化点の可視化
cpts(res_cpts)    # 変化点

変化点は、年号で言うと延喜/延長の間・康応 / 明徳の間・元治 / 慶応の間で起こっているようです。ちなみに、

・延喜 / 延長の間 : 辛酉の年に改元する習わしを決めた時
・康応 / 明徳の間 : 室町時代の南北朝の統一
・元治 / 慶応の間 : 大政奉還

が起こった頃で、何か関係があるのかなと考えているところです。

§3. 時系列モデリング

変化点ごとに期間を切って、年号の続いた年数を調べてみることにしました。今回は特に、辛酉の年に改元する習わしが決まった延喜の直後になる延長から室町時代に南北朝統一が行われた明徳までの期間でモデリングしてみましょう。（この期間を選んだのは、時点数が143個と一番多かったからです。）

あとで、モデルを評価するために最後の9時点分をテストデータとして分割し、残りの134時点でモデルを推定します。

# train test split
n_years_term <- as.ts(dat_nengo[30:172, ]$years)
train <- window(n_years_term, 1, 134)
test <- window(n_years_term, 135, 143)

まず、この期間の時系列を改めて確認しておきましょう。

# 原系列のプロットと差分系列のプロット
plot(train)    # 原系列のプロット
plot(diff(train))    # 差分系列のプロット

次に、差分系列のコレログラムを描いて自己相関をチェックします。

# 延長から明徳までの差分系列のオブジェクト
train_diff <- diff(train)

# コレログラム
acf(train_diff)

差分系列はラグ3までの自己相関を持っていそうに見えます。これから、恐らく原系列はARIMA(0, 1, 3)モデルであろうことが分かります。実際、AICが最も小さくなるようなARIMAモデルを探してみましょう。

＃ 必要なライブラリの読み込み
library(forecast)

# ARIMAモデルを探索する。
res <- auto.arima(train, ic = "aic",
                  stepwise = FALSE, approximation = FALSE,
                  max.p = 10, max.q = 10, max.order = 20, 
                  parallel = TRUE, num.cores = 4)
res

結果、やはりARIMA(0, 1, 3)が得られ、以下のような推定結果が得られました。

・ラグ1のMAの係数 : -0.8051, s.e. = 0.0863
・ラグ2のMAの係数 : +0.1249, s.e. = 0.1043
・ラグ3のMAの係数 : -0.2513, s.e. = 0.0834
・誤差の分散 : 4.305

原系列の数の大きさに比べて誤差分散の大きさがやや大きい点が気になるところです。年号ごとのビッグイベントの有無を表す変数を外生変数として追加することで、モデルをrefineすることができるかもしれません。

§4. 残差分析

推定されたモデルの妥当性を検討するために、残差分析を行いましょう。

# 残差分析
tsdiag(res)

最初の残差系列のプロットから誤差の等分散性は確認できると判断してよさそうです。また残差系列の自己相関がラグ0でしか見られず、Ljung-Box検定でも任意のラグで有意でないことから、残差系列の白色性も確認できると考えてよさそうです。

§5. 予測精度の評価

さらに、モデルの予測精度も確認しておきましょう。

# ARIMA(0, 1, 3)モデルを推定して、h時点後までを予測する関数を作る。
pred_func <- function(x, h){
  res <- Arima(x, order = c(0,1,3))
  return(forecast(res, h = h))
}

# time series cross validation
error_all <- tsCV(n_years_term, pred_func, h = 1)
error_test <- window(error_all, 135, 143)
var(error_test)    # 残差の分散
c(mean(abs(error_test)), mean(abs(error_test/test)))    # MAE, MAPE

time series cross validationの結果、1時点後の予測に対するMAEは1.68, MAPEは41.6%でした。推定の結果得られた時系列モデルの予測精度はそこまで高いものではないと言えそうです。一方で残差の分散は4.86とモデル推定時の誤差の分散と近く、過剰適合なモデルではないと考えられます。

§ まとめ

思いついてから急いでやってみたという記事になってしまいました。もう少し適切に分析プロセスを回せる部分がある気がするので、後でしっかり書き換えたいと思います。

謝辞 : この記事はピースオブケイクさんでの統計勉強会で時系列分析を教える際のイントロダクションとして書いたものです。勉強会参加のみなさまへ、このような機会を頂いていることに心から感謝申し上げます。また、@hanaoriのおかげでスクリプトをrefineする切っ掛けを得ました。ありがとうございます。