ほくそ笑む

R言語と統計解析について

ARIMAX で祝日効果を盛り込んだ時系列予測モデルの作成

書籍紹介 統計

1. はじめに

うちのブログは平日のアクセス数が休日の 2 倍くらいあります。
みなさんお仕事で必要になったときに検索されて、このブログたどり着くのでしょうか。お疲れ様です。

さて、『データサイエンティスト養成読本 R活用編』という書籍で、ARIMAX モデルを用いた時系列分析のやり方が載っています。

データサイエンティスト養成読本 R活用編 【ビジネスデータ分析の現場で役立つ知識が満載! 】 (Software Design plus)

データサイエンティスト養成読本 R活用編 【ビジネスデータ分析の現場で役立つ知識が満載! 】 (Software Design plus)

今日は、この書籍を参考に、うちのブログのアクセス数を ARIMAX モデルを用いて予測してみようと思います。

2. データの準備

まず、データですが、2014/04/01 〜 2014/11/30 までのブログのアクセス数(セッション数)を Dropbox に置いておきます。

ここからデータを取得し、読み込んでください。

library(xts)
data <- read.csv("https://dl.dropboxusercontent.com/u/432512/20150127_arimax/sessions.csv", row.names=1)
data.xts <- as.xts(data)
plot(data.xts["2014-06-01::2014-06-30"], ylim=c(0,600), type="b", main="セッション数(6月)")

6 月のデータをプロットすると、休日のアクセス数が平日に比べて落ちていることが見てとれると思います。

3. ARIMA モデルによる予測

まずは祝日を考慮せずに ARIMA モデルによる予測をやってみようと思います。
予測は、4月〜10月までのデータを用いて11月のセッション数を予測します。
最初に、訓練データとテストデータに分けます。

data.train <- data.xts["::2014-10-31"]
data.train.ts <- ts(data.train, frequency=7)
data.test <- data.xts["2014-11-01::"]
date.test <- as.Date(index(data.test))

訓練データを用いて ARIMA モデルを作成し、予測値を算出します。

library(forecast)
fit.arima <- auto.arima(data.train.ts)
fc.arima <- forecast(fit.arima, h=30)
pred.arima <- as.vector(fc.arima$mean)

res <- data.frame(date=date.test, pred.arima=pred.arima, act=as.numeric(data.test))

データを整形し、ggplot2 によるグラフの描画を行います。

library(tidyr)
library(ggplot2)

res.m <- gather(res, type, value, -date) 
ggplot(data=res.m, aes(x=date, y=value, color=type)) + geom_line(size=1.5) + ggtitle("ARIMA")

青が実測値、赤が予測値です。
あんまり予測精度はよくありませんが、そんなことより注意して見てほしいのは、11月には祝日が 2日あり(11/3, 11/24)、その日は実測値と予測値が大きく外れているということです。
これは、祝日による効果を ARIMA モデルが考慮できていないことから生じます。
次に見る ARIMAX モデルでは、この祝日効果をモデルに取り入れます。
とりあえずモデルの評価指標として RMSE を計算しておきます。

rmse <- sqrt(mean((res$pred.arima - res$act)^2))
print(rmse)

[1] 73.16669

4. ARIMAX モデルによる祝日効果の取り込み

それでは、ARIMAX モデルによる予測をやってみます。
日本の祝日判定は、Nippon パケージの is.jholiday() 関数を使って行います。

library(Nippon)
is.holiday <- function(date) !(weekdays(date) %in% c("土曜日", "日曜日")) & is.jholiday(date)

dates <- as.Date(index(data.xts), tz="Japan")
jholidays <- data.frame(row.names=dates, is.holiday=is.holiday(dates))
jholidays.xts <- as.xts(jholidays)

祝日効果を取り込んだ ARIMAX モデルを作成し、予測値を出します。

fit.arimax <- auto.arima(data.train.ts, xreg=jholidays.xts["::2014-10-31"])
fc.arimax <- forecast(fit.arimax, xreg=jholidays.xts["2014-11-01::"], h=30)
pred.arimax <- as.vector(fc.arimax$mean)

res <- data.frame(date=date.test, pred.arimax=pred.arimax, act=as.numeric(data.test))

データを整形し、ggplot2 でグラフを描画します。

res.m <- gather(res, type, session, -date)
ggplot(data=res.m, aes(x=date, y=session, color=type)) + geom_line(size=1.5) + ggtitle("ARIMAX")

おー、11/3 と 11/24 の予測値がちゃんと変化しましたね！

モデルの精度を比較するために RMSE を計算します。

rmse <- sqrt(mean((res$pred.arimax - res$act)^2))
print(rmse)

[1] 63.5615

RMSE が下がり、予測精度が良くなったことがわかります。

4. まとめ

ARIMAX モデルを用いて祝日効果を考慮した予測モデルを作成しました。
書籍では、さらに年月日や曜日などの影響も考慮したモデルを作成しています。
気になる方はぜひ書籍の方を読んでみてください。

データサイエンティスト養成読本 R活用編 【ビジネスデータ分析の現場で役立つ知識が満載! 】 (Software Design plus)

データサイエンティスト養成読本 R活用編 【ビジネスデータ分析の現場で役立つ知識が満載! 】 (Software Design plus)

以上です。

追記

auto.arima() を使わず地道にやるとこんな感じに。

auto.arima() はあまり評判が良くないようです。