RでARモデルの勉強
時系列分析の勉強をはじめることにしました。まずはARモデルです。
以下の記事を見ながらARモデルの構築と、モデル検証の勉強をしました。
データ
Google Trendsのヒット数を使用します。検索ワードは”温泉”にしました。
期間は2016年1月1日～2019年1月1日のdailyデータで、データの取得には、以下で作成した自作ツールを使用しました。
hadadada00.hatenablog.com
基礎分析
取得したcsvファイルをRに読み込んで簡単にトレンドと季節性がないか確認しました。
最初にトレンドを確認します。
# import pacages ---------------------------------------------------------- library(lubridate) library(tidyverse) # read -------------------------------------------------------------------- onsen <- read_csv("./onsen_20160101-20190101.csv") %>% mutate(date = as.Date(date)) # check data -------------------------------------------------------------- # 1. trend onsen %>% ggplot(aes(x = date, y = hits_adj)) + geom_line() + scale_x_date(date_labels = "%Y%m%d")
全体的に少し右上がりか横ばいです。
次に曜日による季節性を確認しました。
# 2. weekly seasonality onsen %>% mutate(week_day = wday(date, label = TRUE)) %>% ggplot(aes(x = week_day, y = hits_adj)) + stat_summary(fun.y = "mean", color = "red", geom = "point") + stat_summary(fun.y = "median", color = "blue", geom = "point")
明らかに曜日による違いがあります。日曜日と土曜日が多く、次に月曜と金曜が少し多めで、火曜、水曜、木曜が少ないです。
温泉というレジャーに関する検索ヒット数のため、休日に数が増えるのは妥当だと思います。
最後に月による季節性を確認します。
# 3. monthly seasonality onsen %>% mutate(month = month(date, label = TRUE)) %>% ggplot(aes(x = month, y = hits_adj)) + stat_summary(fun.y = "mean", color = "red", geom = "point") + stat_summary(fun.y = "median", color = "blue", geom = "point")
こちらも月による季節性がありそうです。１月、８月が多いのは冬休みと夏休みがあり温泉旅行に関心が高まるからだと思います。（夏でも多いのが少し意外です。）それ以外は概ね暑い時期は少なく、寒い時期は多いという傾向が見られます。
モデル構築
まずはAR(p)のp(次数)の目安を得るために、自己相関係数と偏自己相関係数を見ます。
stats::acfとstats::pacf関数をそれぞれ使用します。
自己相関係数を見ると、自己相関が存在し、特に７日前、１４日前、と１週間区切りで相関が高いことが確認されます。
次に偏自己相関係数を確認します。AR(p)モデルの次数はこちらの結果が参考になります。
（例えば、自己相関係数だとt=1とt=7の相関を見るときにt=2, t=3など他の時点の影響が含まれてしまいます。本当はt=1, t=7間に直接相関がなかったとしてもt=1, 2, 3, 4, 5, 6, 7と影響が連鎖していき、見た目上t=1とt=7に相関があるように見えてしまいます。そのため次数の参考には偏自己相関係数を使用します。）
結果を見るとt=-1に大きなプラスの相関があります。
以降はt=-6, -12, -13, -20と概ね6～8日おきに正の相関がある一方、t=-8, -15, -22は負の相関がでております。
前日と、概ね1週間単位のラグに影響を受けていることがわかります。
モデリング
stats::ar関数でモデルを構築します。デフォルトでは次数がAICを基準に自動的に選ばれます。引数としては最大次数のみ個別に設定します。
デフォルトでは以下の記述にあるように、min(N-1, 10*log10(N) N=number of observations の値がセットされるため、2016年1月1日～2019年1月1日のdailyデータの場合、N=1097で、デフォルトの最大次数は30になります。
order.max
maximum order (or order) of model to fit. Defaults to the smaller of N-1 and 10*log10(N) where N is the number of non-missing observations except for method = "mle" where it is the minimum of this quantity and 12.
基礎分析で月の季節性があることが確認できたため、最大1年分の影響を検討対象にするべく、最大次数は400としました。
# model building ---------------------------------------------------------- # AR model ar_model <- ar(onsen$hits_adj, order.max = 400) summary(ar_model) > summary(ar_model) Length Class Mode order 1 -none- numeric ar 30 -none- numeric var.pred 1 -none- numeric x.mean 1 -none- numeric aic 401 -none- numeric n.used 1 -none- numeric n.obs 1 -none- numeric order.max 1 -none- numeric partialacf 400 -none- numeric resid 1097 -none- numeric method 1 -none- character series 1 -none- character frequency 1 -none- numeric call 3 -none- call asy.var.coef 900 -none- numeric
summary()の実行結果のうち、arが次数を表わしています。arの値は30のため、AR(30)のモデルが選択されました。
モデル検証
以下の2点を確認します。
- 誤差項の独立性
- 誤差項の正規性
これらはいずれもAR(p)モデルの前提になっている条件です。
まずは誤差項の独立性です。Ljung-Box検定を行います。
# check model ------------------------------------------------------------- # 1. Ljung-Box test(1978) (independency of error term) # H0 : data doesn't show auto correlation # small p-value means there is a autocorrelation Box.test(ar_model$res, type = "Ljung") >Box.test(ar_model$res, type = "Ljung") Box-Ljung test data: ar_model$res X-squared = 3.7343, df = 1, p-value = 0.05331
Ljung-Box検定の帰無仮説は「誤差項に自己相関が存在しない」です。
p値はは0.053のため、5%有意水準で帰無仮説は棄却されませんが、10%有意水準では棄却されます。そのため、検定結果としては誤差項に自己相関がないとはいえないです。
次に誤差項が正規分布に従うか検証します。Jauque-Bera検定を行います。
このとき、入力データに欠損がでないよう気をつけます。
# 2. Jauque-Bera test(1987) (normality of error term) # H0 : data has a normal distribution # small p-value means the data doen't follow normal distribution # before testing model, remove NA data in dataset jarque.bera.test(ar_model$res %>% tail(n = -30)) >jarque.bera.test(ar_model$res %>% tail(n = -30)) Jarque Bera Test data: ar_model$res %>% tail(n = -30) X-squared = 1533.9, df = 2, p-value < 2.2e-16
p値は小さいため、有意水準1%以下で棄却され、誤差項は正規分布に従わないという結果がでました。
検証結果を見ると、誤差項にまだ自己相関がのこっており、正規分布にも従っていないため、モデリングにはまだ足りない要素があるといえます。
今回はモデリングの手順確認が主目的のため、追加のモデル構築は行いません。
予測
構築された（イマイチな）モデルを使用して、予測を行ってみます。
ただし今回はモデルの精度を検証するというよりはARモデルで予測するとどのような予測値が得られ、どのような限界があるかを確かめることに重点をおきます。
そのため以下のようなやり方でできるだけ長期間の予測値を得られるようにします。
- 2016-01-01～2016-01-30の実際の値(hits_adj)を使用して以降の期間の予測（hits_hat）を作成する。
- 2016-01-31の値は2016-01-01～2016-01-30の実測値を使って予測する。
- 2016-02-01の値は2016-01-02～2016-01-30の実測値と2016-01-31の予測値を使用して予測する。
- 以降同様に予測に使用するデータを実測値から予測値にずらしていく。
# predict ----------------------------------------------------------------- input <- head(onsen$hits_adj, n = 30) pred <- predict(ar_model, newdata = input, n.ahead = (length(onsen$hits_adj) - 30)) onsen_hits_hat <- union_all(rep(NA, 30), pred$pred) # data frame for plot df <- data.frame(date = onsen$date, hits_adj = onsen$hits_adj, hits_hat = onsen_hits_hat) # plot true value and predicted value df %>% ggplot() + geom_line(aes(x = date, y = hits_adj)) + geom_line(aes(x = date, y = hits_hat), color = "red") + xlim(as.Date("2016-01-01"), as.Date("2019-01-01"))
実際の値と予測値をプロットすると以下のようになります。
と予測値の比較（３年分）
まず3年後は明らかに水準も違いますし、上下の変動幅も全く異なります。
描画期間を変更して、短期間での予測値の動きを見てみます。
# plot true value and predicted value df %>% ggplot() + geom_line(aes(x = date, y = hits_adj)) + geom_line(aes(x = date, y = hits_hat), color = "red") + xlim(as.Date("2016-01-01"), as.Date("2019-07-01"))
描画期間を2016-01-01～2016-07-01に変更しています。予測値は2016-01-31からのため、概ね半年間の予測をしています。
結果をみると、最初の3ヶ月あたりは、上昇の動きと下落の動きを捉えられており、かつピークとボトムの水準の誤差はそれほど大きくないため、あれば予測値としてある程度は使えると思います。
しかし、3ヶ月を過ぎると、上昇、下落のタイミングは捉えられているものの、全体的な水準が異なっているため、予測値としては使えません。
今回構築されたモデルがAR(30)のため、曜日の季節性はとらえられているものの、月の季節性が捉えられていないことが原因と思われます。
とはいえ、必要な予測期間が3ヶ月程度であればARモデルというシンプルなモデルでも、それなりに使える予測値が得られるということがわかりました。
