時系列データとは
日あるいは月、時・分・秒など一定の間隔で取られた一連のデータ。
時系列データでないデータを、トランザクションデータと呼ぶ。
時系列データの取り扱い
- 「一日のデータ」は「一日に一回しか手に入らない」。
- 2000年1月1日の気温データは1つしか存在しないが、推測統計の考え方を適用し、2000年1月1日という日が無数にあったと仮定したとき、「無数に存在する2000年1月1日の気温」が母集団となる。手元にある一つの2000年1月1日のデータから、母平均を推定する、つまり「無数に存在する2000年1月1日の気温の平均」を求めるようなことをしたい。
データ生成過程（Data Generation Process）
時間に従って変化する確率分布の事。確率過程、あるいは単に過程とも呼ばれる。
自己相関
時系列データの特徴は、データに前後の関係があること。
自己相関とは、過去と未来の相関を取ったもののこと。
- 正の自己相関：昨日の気温が高ければ今日も高い
- 負の自己相関：昨日の気温が高ければ今日は低い
自己相関の情報を未来予測に利用できる。
ホワイトノイズ
未来を予測する情報がほとんど含まれていない純粋な雑音。
「期待値0、分散が一定、自己相関が0」であることが要件。
# ホワイトノイズの例
white_noise <- rnorm(n = 400)
plot(white_noise, type = 'l')
コレログラム
時系列データにおいて、どのくらい離れた時期とどの程度相関があるのかをグラフにしたもの。
Rだとacf関数でコレログラムを描くことができる。
acf(white_noise)
ランダムウォーク
データの変化がランダムである過程。ホワイトノイズのような乱数の累積和のようなもの。
ホワイトノイズの累積和を取ると簡単にランダムウォークが作れる。
random_walk <- cumsum(white_noise)
plot(random_walk, type = 'l', main = 'ランダムウォーク')
時系列データの構造
時系列データは以下の要素に分解できる。
- 短期の自己相関
- 周期性変動
- トレンド
- 外因性
- ホワイトノイズ
その他用語
| 用語 | 英語 | 意味 |
|---|---|---|
| 定常過程 | Stationary Process | 時系列の性質が変化しない過程 |
| 非定常過程 | Non-Stationary Process | 時系列の性質が変化する過程 |
| ノイズ | Noise | 他の成分で説明できないデータのばらつき、誤差 |
| 単位根検定 | Unit root test | 時系列データが定常過程であることの検定 |
| 差分 | Difference | ある時点のデータとn個前のデータとの差、変化量 |
| 偏自己相関 | Partial autorecorrelation | 間のデータの影響を除去した自己相関 |
| 残差 | Redidual | 実データとモデルの推定値との差 |
Box-Jenkins法（時系列分析のフレームワーク）
| ステップ | 内容 | ポイント |
|---|---|---|
| 1 | データを分析しやすいように変換する（差分・対数・対数差分など） | 定常過程への変換 |
| 2 | データにARIMAモデルやそれに準ずつモデルを適用する | ・ARIMA,SARIMAなどの仕組み。 ・モデル同定の手順 |
| 3 | 推定されたモデルを評価する | 評価の方法（残差の自己相関、残差の正規性の検定） |
| 4 | 推定されたモデルを用いて予測する | 評価の方法（予測精度の評価） |
Rにおける時系列データの取り扱い
ts型
Rの標準パッケージに含まれる時系列データの方としてts型がある（time series）
# 例として、2000年1月から月単位のデータを36個生成してみる。
# 月単位のデータにしたい場合、freq = 12
> ts_sample <- ts(1:36, start = c(2000, 1), frequency = 12)
> ts_sample
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
2000 1 2 3 4 5 6 7 8 9 10 11 12
2001 13 14 15 16 17 18 19 20 21 22 23 24
2002 25 26 27 28 29 30 31 32 33 34 35 36
# 四半期ごとに区切る場合freq = 4とする
> ts_sample <- ts(1:24, start = c(2000, 1), frequency = 4)
> ts_sample
Qtr1 Qtr2 Qtr3 Qtr4
2000 1 2 3 4
2001 5 6 7 8
2002 9 10 11 12
2003 13 14 15 16
2004 17 18 19 20
2005 21 22 23 24
# window関数で特定の期間を抽出
> window(ts_sample, start = c(2000, 3), end = c(2000, 10))
Qtr1 Qtr2 Qtr3 Qtr4
2000 3 4
2001 5 6 7 8
2002 9 10
xts型
ts型の欠点として、日単位のデータに弱いという点がある。
Rにおいては、他にもzooなどの時系列データを扱う型があるが、それらを統合している型がxts型。
xts型のデータを作成する方法は複数あるが、行名として日付を指定したmatrix型データを引数に与えるのが簡単。
> mat <- matrix(
c(1:5),
dimnames = list(
as.character(seq(as.Date("2000-01-01"), as.Date("2000-01-05"), by = "day"))),
ncol = 1
)
> xts_sample <- as.xts(mat)
> xts_sample
[,1]
2000-01-01 1
2000-01-02 2
2000-01-03 3
2000-01-04 4
2000-01-05 5
# 日付を指定してデータ取得
> xts_sample["2000-01-03"]
[,1]
2000-01-03 3
# コロンで範囲指定も可
> xts_sample["2000-01-03::"]
[,1]
2000-01-03 3
2000-01-04 4
2000-01-05 5
> xts_sample["2000-01-02::2000-01-04"]
[,1]
2000-01-02 2
2000-01-03 3
2000-01-04 4
# 下記のようなdata.frameをxtsに変換する場合、read.zooを介すとよい
> df
time date
1 2000-01-01 1
2 2000-01-02 2
3 2000-01-03 3
4 2000-01-04 4
5 2000-01-05 5
> df_xts <- as.xts(read.zoo(df))
> df_xts
[,1]
2000-01-01 1
2000-01-02 2
2000-01-03 3
2000-01-04 4
2000-01-05 5
グラフ描画
R組み込み「EuStockMarkets」(Daily Closing Prices Of Major European Stock Indices, 1991--1998)のドイツの株価（DAX）を描画してみる。
plot(EuStockMarkets[,'DAX'],
ylab = 'DAX',
main = 'ドイツにおける株価')
ggplot2とggfortifyパッケージ読み込みでよりきれいなグラフを描ける。
autoplot(EuStockMarkets[,'DAX'],
xlim = c(1991,1999),
xlab = 'time',
ylab = 'DAX',
main = 'ドイツにおける株価') +
scale_x_continuous(breaks = seq(1991,1999,1))
decompose関数を使うことで、時系列データをトレンド・季節成分・ノイズなどに分解できる。
decomp <- decompose(EuStockMarkets[,'DAX'])
plot(decomp)
forecastパッケージのggtsdisplay関数を使うと、データ系列と相関係数（ACF)、偏相関係数（PACF）のコレログラムを描くことができる。
# 簡単のため先頭100日分のデータで表示
> ggtsdisplay(EuStockMarkets[1:100,'DAX'], main = 'ドイツの株価DAX')
左下の相関係数のグラフを見ると、長期にわたり自己相関が続いていることが分かる。
右下の偏自己相関では、およそ5日周期で大きな偏自己相関がみられる。
定常/非定常過程の判別：単位根検定
データの差分を取る必要があるかどうかを判別する方法に、単位根検定がある。
| 検定 | 帰無仮説 | 対立仮説 | 判断 |
|---|---|---|---|
| ADF検定 | 単位根あり | 単位根なし | 危険率5%で有意となったら差分は取らないと判断 |
| KPSS検定 | 単位根なし | 単位根あり | 危険率5%で有意となったら差分を取るべきと判断 |
ADF検定の例
> library(tseries)
> adf.test(EuStockMarkets[,'DAX'])
Augmented Dickey-Fuller Test
data: EuStockMarkets[, "DAX"]
Dickey-Fuller = -0.82073, Lag order = 12, p-value = 0.9598
alternative hypothesis: stationary
# p値0.9598となり、帰無仮説は棄却されない→単位根あり→差分を取るべき（非定常過程）
# 差分を取って再検証
> adf.test(diff(EuStockMarkets[,'DAX']))
Augmented Dickey-Fuller Test
data: diff(EuStockMarkets[, "DAX"])
Dickey-Fuller = -9.9997, Lag order = 12, p-value = 0.01
alternative hypothesis: stationary
Warning message:
In adf.test(diff(EuStockMarkets[, "DAX"])) :
p-value smaller than printed p-value
# p値0.01となり帰無仮説は棄却→単位根なし→差分は取らなくてよい（定常過程）
差分を取ったデータを図示してみる。
> plot(diff(EuStockMarkets[,'DAX']))
KPSS検定の例
> library(urca)
> summary(ur.kpss(EuStockMarkets[,'DAX']))
#######################
# KPSS Unit Root Test #
#######################
Test is of type: mu with 8 lags.
Value of test-statistic is: 15.4007
Critical value for a significance level of:
10pct 5pct 2.5pct 1pct
critical values 0.347 0.463 0.574 0.739
# 有意水準を5%としたときの棄却点が0.463であり、棄却点15.4007を上回っている。
# 単位根がないという帰無仮説が棄却されたので単位根があると見なす。（差分を取る必要がある、非定常過程）
# ndiffs関数で何回差分を取るべきか判定することができる。
> library(forecast)
> ndiffs(EuStockMarkets[,'DAX'])
[1] 2
# 1階差分
> summary(ur.kpss(diff(EuStockMarkets[,'DAX'])))
#######################
# KPSS Unit Root Test #
#######################
Test is of type: mu with 8 lags.
Value of test-statistic is: 0.7564
Critical value for a significance level of:
10pct 5pct 2.5pct 1pct
critical values 0.347 0.463 0.574 0.739
# 2階差分
> summary(ur.kpss(diff(diff(EuStockMarkets[,'DAX']))))
#######################
# KPSS Unit Root Test #
#######################
Test is of type: mu with 8 lags.
Value of test-statistic is: 0.0112
Critical value for a significance level of:
10pct 5pct 2.5pct 1pct
critical values 0.347 0.463 0.574 0.739
# 確かに2回差分で、0.463を下回った。
# 単位根がないという帰無仮説がされないので単位根がないと見なす。（差分を取る必要がない、定常過程）
<補足>
| 前処理の種類 | 使用シーン |
|---|---|
| 対数変換 | 時系列データの分散が一定ではない時に行うことで、分散を一定と見なせる場合がある |
| 差分 | 差分はデータのトレンドを除くために行う。トレンドが除去できない場合は、二階差分を行う場合もある。 |
| 対数差分 | データの変動を明らかにしたい時に用いる。一つ前の値と次の値の比をとったものを対数変換しているので、変動を表す。 $ \Delta(log(y_{t})) = log(y_{t}) - log(y_{t-1}) = log \dfrac{y_{t}}{y_{t-1}}$ |
時系列モデルについて
手元の時系列データから未来のデータの傾向、変化を予測するために様々なモデルが提案されている。
| モデル | 数式表現 | 特徴 |
|---|---|---|
| AR(Autoregression) | AR(p) $ y_{t} = \phi_{1}y_{t-1}+\phi_{2}y_{t-2}+\dots\phi_{p}y_{t-p} + \epsilon_{t}$ |
・定常過程が対象 ・p次前のデータで回帰 ・過去の自分のデータを説明変数として回帰モデルを作成することで自己相関を表現 |
| MA(Moving average) | MA(q) $ y_{t} = \mu + \theta_{1}\epsilon_{t-1} + \theta_{2}\epsilon_{t-2} \dots \theta_{q}\epsilon_{t-q} + \epsilon_{t}$ |
・定常過程が対象 ・q次分の自己相関で予測 ・過去と未来で共通の値を使用する事で自己相関を表現 |
| ARMA(Autoregressive moving average) | ARMA(p,q) $ y_{t} = \phi_{1}y_{t-1}+\phi_{2}y_{t-2}+\dots\phi_{p}y_{t-p} + \theta_{1}\epsilon_{t-1} + \theta_{2}\epsilon_{t-2} \dots \theta_{q}\epsilon_{t-q} +\epsilon_{t}$ |
・定常過程が対象 ・ARとMAの両方のモデルを適用することで、自己相関をより柔軟に表現する |
| ARIMA(Autoregressive integrated moving average) | ARIMA(p,d,q) $ y_{t} - y_{t-d} =c + \phi_{1}y_{t-1} + \phi_{2}y_{t-2} + \cdots + \phi_{p}y_{t-p} + \varepsilon_{t} + \theta_{1}\varepsilon_{t-1} + \cdots + \theta_{q}\varepsilon_{t-q} $ または $ y_{t} - y_{t-d}=c + \varepsilon_{t} + \displaystyle \sum_{ i = 1 }^{ p } \phi_{i}y_{t-i} +\displaystyle \sum_{ i = 1 }^{ q } \theta_{i}\varepsilon_{t-i} $ |
・時系列データのd階差分系列、$y_{t}−y_{t−d}$をARMAモデルで表現するモデル ・非定常過程も対象 ・ARMAに差分を取る操作を追加している ・d階和分過程I(d)においてARMA(p,q)を適用する |
モデルの検証
データセット
R組み込みのデータSeatbelts（イギリスの交通事故死傷者数のデータ）のfront（前席における死傷者数）を使用して予測してみる。
# データの確認
> Seatbelts %>% head()
DriversKilled drivers front rear kms PetrolPrice VanKilled law
Jan 1969 107 1687 867 269 9059 0.1029718 12 0
Feb 1969 97 1508 825 265 7685 0.1023630 6 0
Mar 1969 102 1507 806 319 9963 0.1020625 12 0
Apr 1969 87 1385 814 407 10955 0.1008733 8 0
May 1969 119 1632 991 454 11823 0.1010197 10 0
Jun 1969 106 1511 945 427 12391 0.1005812 13 0
> Seatbelts[,'front']
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
1969 867 825 806 814 991 945 1004 1091 958 850 1109 1113
1970 925 903 1006 892 990 866 1095 1204 1029 1147 1171 1299
1971 944 874 840 893 1007 973 1097 1194 988 1077 1045 1115
1972 1005 857 879 887 1075 1121 1190 1058 939 1074 1089 1208
1973 903 916 787 1114 1014 1022 1114 1132 1111 1008 916 992
1974 731 665 724 744 910 883 900 1057 1076 919 920 953
1975 664 607 777 633 791 790 803 884 769 732 859 994
1976 704 684 671 643 771 644 828 748 767 825 810 986
1977 714 567 616 678 742 840 888 852 774 831 889 1046
1978 889 626 808 746 754 865 980 959 856 798 942 1010
1979 796 643 794 750 809 716 851 931 834 762 880 1077
1980 748 593 720 646 765 820 807 885 803 860 825 911
1981 704 691 688 714 814 736 876 829 818 942 782 823
1982 595 673 660 676 755 815 867 933 798 950 825 911
1983 619 426 475 556 559 483 587 615 618 662 519 585
1984 483 434 513 548 586 522 601 644 643 641 711 721
ggtsdisplay(Seatbelts[,'front'],main='イギリスの交通事故死傷者(前席)')
12か月周期で自己相関が大きくなっており、偏自己相関も12か月単位で大きな偏自己相関がみられる。
1年単位での相関が目立つため、季節成分があることは確かなので、ggsubseriesplot関数を用いて月ごとのグラフを描いてみる。
ggsubseriesplot(Seatbelts[,'front'],main='月毎の死傷者数')
ARモデル
過去p点分のデータ予測に使用するモデルがARモデル。
front_diff <- diff(Seatbelts[,'front'], lag = 12)# 1階差分データを用いる（12か月の季節差分）
> ndiffs(front_diff)# 定常過程かどうかの検証、何回差分を取るか
[1] 0
front_diff_train <- window(front_diff, start = c(1970, 1), end = c(1983,12))# trainデータの抽出
result <- ar(front_diff_train, method = 'ols')# ols（最小二乗法）でARモデル作成
result_pr <- predict(result, front_diff_train, n.ahead = 12)# モデルをもとに先12か月分を予測
result_pred <- ts(result_pr$pred, start = c(1984, 1), frequency = 12)
plot(front_diff)
lines(result_pred, col="red")
legend("topright", legend = c('actual','predict'),lty=c(1,1),col=c(rgb(0,0,0),rgb(1,0,0)))
MAモデル
q次分のデータの自己相関からモデルを構築。
result <- Arima(front_diff_train, order = c(0, 0, 1))
result_pred <- ts(forecast(result, h = 12)$mean,c(1984, 1), frequency = 12)
result_pr <- predict(result, front_diff_train, n.ahead = 12)# モデルをもとに先12か月分を予測
result_pred <- ts(result_pr$pred, start = c(1984, 1), frequency = 12)
plot(front_diff)
lines(result_pred, col="red")
legend("topright", legend = c('actual','predict'),lty=c(1,1),col=c(rgb(0,0,0),rgb(1,0,0)))
ARIMAモデル
データの和分（差分）にARモデルとMAモデルを統合して適用するモデル。
※ARIMAモデルは内部で勝手に差分を取るため、差分系列は使用しない
#最後の1年をテストデータ、それより前を訓練データ
Seatbelts_set <- Seatbelts[,c('front', 'PetrolPrice', 'law')]
train <- window(Seatbelts, end=c(1983,12))
test <- window(Seatbelts, start=c(1984,1))
#説明変数の取り出し
oth <- train[,c('PetrolPrice', 'law')]
#モデル作成
model_sarimax <- Arima(y = train[, 'front'], # 目的変数
order = c(1, 1, 1), # SARIMA(p,d,q)(P,D,Q)における(p,d,q)
seasonal = c(1, 0, 0), # 季節成分(P,D,Q)
xreg = oth) # 説明変数
> model_sarimax
Series: train[, "front"]
Regression with ARIMA(1,1,1)(1,0,0)[12] errors
Coefficients:
ar1 ma1 sar1 PetrolPrice law
0.2837 -0.9453 0.6634 -2805.0083 -243.6760
s.e. 0.0832 0.0325 0.0554 863.4417 41.9162
sigma^2 estimated as 6725: log likelihood=-1044.18
AIC=2100.37 AICc=2100.86 BIC=2119.49
# 両説明変数の係数が負→石油価格が上がったり、法が施行されると事故死傷者が減ると解釈する。
# 作成したモデルをもとに予測
sarimax_f <- forecast(model_sarimax,
xreg = test[,c('PetrolPrice', 'law')],
h = 12, # 何時点先まで予測するか(ここでは12か月=1年)
level = c(95, 70)) # 信頼区間の表示
# 描画
autoplot(sarimax_f,series='Forecast') +
autolayer(test[, 'front'], series = 'Actual') +
autolayer(sarimax_f$mean, series = 'Predict') +
guides(colour=guide_legend(title="Data series"))
モデルの同定
次数の決定
適切なモデルの選定のためにモデルの次数を決定する。
# 自動モデル選択関数のauto.arimaを用いる
model_sarimax_auto <- auto.arima(y = train[, 'front'],
xreg = oth, # 説明変数
ic = "aic", # AICを使ってモデル選択
max.order = 7, # SARIMA(p,d,q)(P,D,Q)におけるmax(p+q+P+Q)。これを増やすほど複雑なモデルを候補に入れてモデル選択できる。
stepwise = F, # Tにすると候補となる次数の組み合わせが減る
approximation = F, # Tにすると毎回の計算で近似的な手法を使うことで計算速度を上げる。
parallel = T, # 並列計算の指定
num.cores = 4) # 4コア
> model_sarimax_auto
Series: train[, "front"]
Regression with ARIMA(2,0,0)(0,1,1)[12] errors
Coefficients:
ar1 ar2 sma1 drift PetrolPrice law
0.3561 0.1938 -0.8621 -1.4226 -3084.3851 -163.1541
s.e. 0.0763 0.0759 0.1279 0.3035 964.8893 45.7199
sigma^2 estimated as 5700: log likelihood=-970.01
AIC=1954.01 AICc=1954.71 BIC=1975.88
残差（予測値－実測値）のチェック
モデルが「良い」モデルかどうかは、残差（予測値－実測値）がきれいにばらついているかどうかを確かめることによって、ある程度判断できる。
予測値の方が実測値よりも常に大きいなどの系統的な誤差があれば失敗。また、残差に自己相関があれば、これはこれでうまくモデリングできていないことになる。誤差はあくまでも、誤差らしくあるような状態が理想。
残差の自己相関の検定には、forecast::checkresiduals関数を使う。
> checkresiduals(model_sarimax_auto)
Ljung-Box test
data: Residuals from Regression with ARIMA(2,0,0)(0,1,1)[12] errors
Q* = 37.794, df = 18, p-value = 0.004124
Model df: 6. Total lags used: 24
p<0.05のため、有意な自己相関がみられるという結果になり、あまりうまくモデリングできていないと解釈する。図示の結果を見てみると、突出した自己相関がところどころ見受けられる。
残差の正規性の検定も行う。tseries::jarque.bera.test関数を使う。残差はresid関数を使えば取得できる。
> jarque.bera.test(resid(model_sarimax_auto))
Jarque Bera Test
data: resid(model_sarimax_auto)
X-squared = 5.9597, df = 2, p-value = 0.0508
こちらはp>0.05となっており、正規分布と有意に異なっているとは言えない結果。
一方で、残差の自己相関に関しては検定をパスしなかったため、モデルの同定手順としては、再度定常過程への変換等を行いモデリングを行う必要がある。
上記検証では、一階差分を取っただけのデータを使用していますが、事前に対数変換を行ったデータで再検証すると以下のようになる。
# 事前に対数変換
Seatbelts_log <- Seatbelts[,c('front', 'PetrolPrice', 'law')]
Seatbelts_log[,'front'] <- log(Seatbelts[,'front'])
Seatbelts_log[,'PetrolPrice'] <- log(Seatbelts[,'PetrolPrice'])
train <- window(Seatbelts_log, end=c(1983,12))
test <- window(Seatbelts_log, start=c(1984,1))
petro_law <- train[,c('PetrolPrice', 'law')]
# auto.arimaによるモデル自動選択
model_sarimax_auto <- auto.arima(y = train[, 'front'],
xreg = petro_law,
ic = "aic",
max.order = 7,
stepwise = F,
approximation = F,
parallel = T,
num.cores = 4)
# 残差の自己相関検定
> checkresiduals(model_sarimax_auto)
Ljung-Box test
data: Residuals from Regression with ARIMA(2,0,1)(0,1,1)[12] errors
Q* = 20.99, df = 18, p-value = 0.2799
Model df: 6. Total lags used: 24
# 残差の正規性の検定
> jarque.bera.test(resid(model_sarimax_auto))
Jarque Bera Test
data: resid(model_sarimax_auto)
X-squared = 0.39938, df = 2, p-value = 0.819
いずれもp>0.05となり、残差に有意な自己相関は見られず、正規性についても有意に異なっているとは言えない結果となった。明らかな問題はないモデルであることが確認できた。
※ただし、今回のモデルでは、説明変数に石油価格（PetrolPrice）を用いてモデルを作成しているが、これは将来の石油価格が分かっている前提のモデルとなってしまっているため、実際は何らかの値で代用してやる必要がある。
予測
# 作成したモデルをもとに予測
sarimax_f <- forecast(model_sarimax_auto,
xreg = test[,c('PetrolPrice', 'law')],
h = 12, # 何時点先まで予測するか(ここでは12か月=1年)
level = c(95, 70)) # 信頼区間の表示
# 描画
autoplot(sarimax_f,series='Forecast') +
autolayer(test[, 'front'], series = 'Actual') +
autolayer(sarimax_f$mean, series = 'Predict') +
guides(colour=guide_legend(title="Data series"))
予測精度の評価
予測精度の評価
sarimax_rmse <- sqrt(sum((sarimax_f$mean - test[,'front'])^2) / length(sarimax_f$mean))
> sarimax_rmse
[1] 0.09674572
# accuracy関数を使うとより簡便
> accuracy(sarimax_f, x = test[, 'front'])
ME RMSE MAE MPE MAPE MASE ACF1 Theil's U
Training set -0.004353509 0.08283297 0.06430968 -0.07784816 0.9593004 0.5848660 -0.004032101 NA
Test set 0.071495376 0.09674572 0.07537019 1.10793155 1.1703728 0.6854562 0.337639488 1.03686