時系列データの分析(arima.simと季節調整法)
UPDATE: 2022-10-28 21:23:10
はじめに
時系列データへの理解、分析方法について、まとめていく。時系列データは、これまでやらないといけないと思いつつも、基礎をすっとばして、Prophet(一般化加法モデルベース)や状態空間モデルをやっていたが、やはり基礎は大事だと思うことがたくさんあったので、基礎からARIMAくらいまでをおさらいする。多分何回にわかれる。
arima.sim()でデータ生成過程
srima.sim()というARIMAモデルからデータを生成し、シミュレーションできる便利な関数がある。これを色々と触ってみる。例えば、下記の設定であれば、自己回帰(AR)モデルが1次で、移動平均(MA)モデルも1次であるARMAモデルからシミュレーションデータが生成される。階差の設定をしていないので、和分過程ではなく、定常過程。
library(tidyverse)
library(forecast)
set.seed(1989)
y <- arima.sim(
n = 500,
model = list(
order = c(1,0,1), # ARIMA(p,d,q)
ar = c(0.5), # arの強さ
ma = c(0.5) # maの強さ
)
)
ggplot(data = tibble::tibble(index = 1:length(y), y = y), aes(index, y)) +
geom_line(col = "#01489D") +
scale_y_continuous() +
ggtitle("Time Series Data ARIMA(1, 0, 1)") +
theme_bw()
下記の設定であれば、自己回帰(AR)モデルが1次で、移動平均(MA)モデルも1次で、階差が1次のARIMAモデルからシミュレーションデータが生成される。階差の設定をしているので、これは1次和分過程。
set.seed(1989)
y <- arima.sim(
n = 500,
model = list(
order = c(1,1,1), # ARIMA(p,d,q)
ar = c(0.5), # arの強さ
ma = c(0.5) # maの強さ
)
)
ggplot(data = tibble::tibble(index = 1:length(y), y = y), aes(index, y)) +
geom_line(col = "#01489D") +
scale_y_continuous() +
ggtitle("Time Series Data ARIMA(1, 1, 1)") +
theme_bw()
季節成分を入れたい場合は、Arima()でモデルを作って、それからそのモデルをsimulate()に入れることで、データを生成すれば季節成分が入ったデータをシュミレーションできる。
model <- Arima(ts(rnorm(100), freq = 12),
order = c(1,0,1),
seasonal = list(order=c(0,1,0), period=12),
fixed=c(phi = 0.5, theta = -0.4), #AR,MAのつよさ
)
y <- simulate(model, nsim=200)
forecast::ggtsdisplay(y, lag.max = 60)
季節調整法
Rには便利な関数decompose()がある。これは時期列データを、移動平均法を使って季節性、トレンド、不規則変動に分解してくれる関数。つまり、時期列データは下記の要素で構成されていると考える。
- 季節性：一定期間繰り返されるパターンのこと。
- トレンド：メトリックの基本的な傾向。
- 不規則変動：ノイズのことで、季節およびトレンドが削除された後の元の時系列の残差。
内容を把握したいので、やっていることを再現する。データを分割する方法は下記のように行われている。中心化移動平均でトレンドを推定する。周期が12ヶ月周期であれば、windowは12で、クオーターとかであれば、windowは4となる。このトレンドをもとに、トレンドを除去したデータを作り、これの月ごとの平均を得ることで、平均的な季節周期の値が得られる。そして、もとの現系列からトレンド、季節周期を除けば、ランダムノイズが得られ分解完了。
df <- tibble::as_tibble(AirPassengers) %>%
dplyr::mutate(time = seq(from = lubridate::ymd("1949-01-01"), by = "1 month", length.out = nrow(.))) %>%
dplyr::rename(y = x)
df_plt <-
df %>%
dplyr::mutate(trend_air = ma(y, order = 12, centre = TRUE),
detrend = y - trend_air) %>%
dplyr::group_by(m = lubridate::month(time)) %>%
dplyr::mutate(seasonal_air = mean(detrend, na.rm = TRUE),
random_air = y - trend_air - seasonal_air,
y_recomp = trend_air + seasonal_air + random_air) %>%
dplyr::select(time, y, trend_air, seasonal_air, random_air, y_recomp) %>%
dplyr::ungroup()
df_plt %>%
print(n = 24)## # A tibble: 144 × 7
## m time y trend_air seasonal_air random_air y_recomp
## <dbl> <date> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 1949-01-01 112 NA -25.5 NA NA
## 2 2 1949-02-01 118 NA -36.9 NA NA
## 3 3 1949-03-01 132 NA -2.99 NA NA
## 4 4 1949-04-01 129 NA -8.79 NA NA
## 5 5 1949-05-01 121 NA -5.26 NA NA
## 6 6 1949-06-01 135 NA 34.7 NA NA
## 7 7 1949-07-01 148 127. 63.1 -41.9 148
## 8 8 1949-08-01 148 127. 62.1 -41.3 148
## 9 9 1949-09-01 136 128. 15.8 -7.73 136
## 10 10 1949-10-01 119 129. -21.4 11.8 119
## 11 11 1949-11-01 104 129 -54.3 29.3 104
## 12 12 1949-12-01 118 130. -29.4 17.6 118
## 13 1 1950-01-01 115 131. -25.5 9.25 115
## 14 2 1950-02-01 126 133. -36.9 29.9 126
## 15 3 1950-03-01 141 135. -2.99 9.08 141
## 16 4 1950-04-01 135 136. -8.79 7.37 135
## 17 5 1950-05-01 125 137. -5.26 -7.16 125
## 18 6 1950-06-01 149 139. 34.7 -24.4 149
## 19 7 1950-07-01 170 141. 63.1 -34.0 170
## 20 8 1950-08-01 170 143. 62.1 -35.2 170
## 21 9 1950-09-01 158 146. 15.8 -3.48 158
## 22 10 1950-10-01 133 148. -21.4 5.98 133
## 23 11 1950-11-01 114 152. -54.3 16.8 114
## 24 12 1950-12-01 140 155. -29.4 14.7 140
## # … with 120 more rows
## # ℹ Use `print(n = ...)` to see more rowsまとめてプロットするとこんな感じ。
# pivot_longerでlongでプロットしたほうがいいけど、tsクラスの消し方がわからん…
ggplot(df_plt) +
geom_line(aes(time, y), col = "#e41749") +
geom_line(aes(time, trend_air), col = "#ff8a5c") +
geom_line(aes(time, seasonal_air), col = "#f5587b") +
geom_line(aes(time, random_air), col = "#8ac6d1") +
theme_bw()
便利な関数があるので、賢い人はこっちを使うこと。
AirPassengers %>% decompose() %>% autoplot()