12.1 複雑な季節性

ここまでは、四半期データや月次データでの比較的単純な季節パターンだけを考えてきました。しかし、より高い頻度の時系列はより複雑な季節パターンをしばしば示します。例えば、日次データには週次パターンも年次パターンもあり得ます。時間ごとのデータには通常、日次パターン、週次パターン、年次パターンの3つのタイプの季節性があります。週次データでさえ、1年間の週の数が整数でなく、年次パターンの周期が平均で\(365.25/7\approx 52.179\)となるため、予測するのが難しくなり得ます。ここまでに検討した手法のほとんどは、こうした複雑な季節性を扱えません。

モデルにあり得る全ての季節周期を含めたいわけではなく、データに存在していそうな季節周期だけを含めたいのです。例えば、180日分のデータしかなければ、年次季節性は無視できます。データが自然現象の計測値(例えば、気温)なら、週次季節性を無視しても問題はないはずです。

12.1は、毎平日午前7時から午後9時5分までに北米のある商業銀行へかかってきた電話の5分間隔ごとの数、33週間分です。下のパネルは、同じ時系列の最初の4週間分です。強い周期169の日次季節パターン(一日当たり169回の5分間隔があります)と、弱い周期\(169 \times 5=845\)の週次季節パターンがあります。(月曜日の電話数が週の残りの曜日より多い傾向があります。)もっと長い系列のデータが入手可能なら、年次季節パターンも観測したかもしれません。

bank_calls %>%
  fill_gaps() %>%
  autoplot(Calls) +
  labs(y = "電話数",
       title = "銀行への5分間隔の電話数")
北米のある大商業銀行の平日午前7時から午後9時5分までに取り扱った5分間隔の電話数。上のパネル: 2003年3月3日から10月24日までのデータ。下のパネル: データの最初の4週間分

図 12.1: 北米のある大商業銀行の平日午前7時から午後9時5分までに取り扱った5分間隔の電話数。上のパネル: 2003年3月3日から10月24日までのデータ。下のパネル: データの最初の4週間分

季節周期が複数あること以外にも、この系列には勤務日と勤務日の間に欠損値があるという追加の複雑さがあります。

複数の季節周期を持つSTL

STL()関数は複数の季節周期を扱えるよう設計されています。トレンド成分、残余成分と一緒に、複数の季節成分を返します。この場合、欠損値を避けるには、tsibbleを再インデックス化してから、明示的に季節周期を与える必要があります。

calls <- bank_calls %>%
  mutate(t = row_number()) %>%
  update_tsibble(index = t, regular = TRUE)
calls %>%
  model(
    STL(sqrt(Calls) ~ season(period = 169) +
                      season(period = 5*169),
        robust = TRUE)
  ) %>%
  components() %>%
  autoplot() + labs(x = "観測")
電話数データへの複数の季節性でのSTL分解

図 12.2: 電話数データへの複数の季節性でのSTL分解

日次パターン(3番目のパネル)と週次パターン(4番目のパネル)の2つの季節パターンが表示されています。このグラフを適切に解釈するには、縦軸の尺度に注意することが大切です。このケースでは、トレンド成分と週次季節性成分に、他の成分と比べて長いバー(と比較的狭い範囲)がありますが、データにトレンドがほとんど見られない上に、週次季節性が弱いためです。

季節成分のそれぞれは季節ナイーブ法で予測し、季節調整値はETSで予測する形で、分解を予測に使うことができます。

以下のコードはいつもより若干複雑になっています。観測値が欠損している期間を扱うためにtsibbleを再インデックス化した際に失った時間スタンプを付け戻しているせいです。STL分解の際に平方根変換しているのは、予測を確実に正値にするためです。

# STL分解 + ETSで予測
my_dcmp_spec <- decomposition_model(
  STL(sqrt(Calls) ~ season(period = 169) +
                    season(period = 5*169),
      robust = TRUE),
  ETS(season_adjust ~ season("N"))
)
fc <- calls %>%
  model(my_dcmp_spec) %>%
  forecast(h = 5 * 169)

# fableに正しい時間スタンプを付け戻す
fc_with_times <- bank_calls %>%
  new_data(n = 7 * 24 * 60 / 5) %>%
  mutate(time = format(DateTime, format = "%H:%M:%S")) %>%
  filter(
    time %in% format(bank_calls$DateTime, format = "%H:%M:%S"),
    wday(DateTime, week_start = 1) <= 5
  ) %>%
  mutate(t = row_number() + max(calls$t)) %>%
  left_join(fc, by = "t") %>%
  as_fable(response = "Calls", distribution = Calls)

# 結果をデータの最後の3週間と合わせてプロット
fc_with_times %>%
  fill_gaps() %>%
  autoplot(bank_calls %>% tail(14 * 169) %>% fill_gaps()) +
  labs(y = "電話数", level = "区間予測",
       title = "銀行への5分間隔の電話数")
#> Warning: Removed 100 rows containing missing values or values outside the
#> scale range (`geom_line()`).
STL分解後、季節成分予測には季節ナイーブ法を使い、季節調整値予測にはETSを使った、電話数データの予測

図 12.3: STL分解後、季節成分予測には季節ナイーブ法を使い、季節調整値予測にはETSを使った、電話数データの予測

複数の季節周期を持つ動学的ハーモニック回帰

複数の季節性があっても、以前の章で行った(7.4節と10.5節、参照)ようにフーリエ項が使えます。複数の季節性があるので、季節周期ごとにフーリエ項を加える必要があります。このケースでは、季節周期は169と845なので、フーリエ項は以下表記のようになります。 \[ \sin\left(\frac{2\pi kt}{169}\right), \quad \cos\left(\frac{2\pi kt}{169}\right), \quad \sin\left(\frac{2\pi kt}{845}\right), \quad \cos\left(\frac{2\pi kt}{845}\right) \] ただし、\(k=1,2,\dots\)です。いつもように、fourier()関数があなたに代わってこれらを生成してくれます。

ARIMA誤差構造を持つ動学的ハーモニック回帰モデルを適合させましょう。季節周期ごとのフーリエ項の総数はAICcを最小化するよう選択することもできます。しかし、季節周期が大きい場合、そうすると必要な項の数を過大推計する傾向があるので、日次季節性には10項、週次季節性には5項と、より主観的に選択することにします。ここでも、平方根変換を使って予測と区間予測が正値になるよう確保しています。非定常性を回帰項に対処させるため\(D=d=0\)、季節性を回帰項に対処させるため\(P=Q=0\)、と設定します。

fit <- calls %>%
  model(
    dhr = ARIMA(sqrt(Calls) ~ PDQ(0, 0, 0) + pdq(d = 0) +
                  fourier(period = 169, K = 10) +
                  fourier(period = 5*169, K = 5)))

fc <- fit %>% forecast(h = 5 * 169)

# fableに正しい時間スタンプを付け戻す
fc_with_times <- bank_calls %>%
  new_data(n = 7 * 24 * 60 / 5) %>%
  mutate(time = format(DateTime, format = "%H:%M:%S")) %>%
  filter(
    time %in% format(bank_calls$DateTime, format = "%H:%M:%S"),
    wday(DateTime, week_start = 1) <= 5
  ) %>%
  mutate(t = row_number() + max(calls$t)) %>%
  left_join(fc, by = "t") %>%
  as_fable(response = "Calls", distribution = Calls)
# 結果を最近3週と合わせてプロット
fc_with_times %>%
  fill_gaps() %>%
  autoplot(bank_calls %>% tail(14 * 169) %>% fill_gaps()) +
  labs(y = "電話数", level = "区間予測",
       title = "銀行への5分間隔の電話数")
電話数データに適用した動学的ハーモニック回帰からの予測

図 12.4: 電話数データに適用した動学的ハーモニック回帰からの予測

これは大きなモデルで、ARMA係数4個、周期169のためのフーリエ係数20個、周期845のためのフーリエ係数8個、切片も加えると合計33個のパラメータがあります。周期845のフーリエ項は、周期169の項と重複する(\(845=5\times169\)なので)ものがあるので、全ては使われていません。

事例: 電力需要

こうしたモデルを適用することが多いのが、電力需要のモデル化です。図12.5は、オーストラリア・ビクトリア州での2012–2014年の30分ごとの電力需要(MWh)とメルボルン(ビクトリア州最大の都市)の同時期の気温(摂氏)です。

vic_elec %>%
  pivot_longer(Demand:Temperature, names_to = "Series") %>%
  ggplot(aes(x = Time, y = value)) +
  geom_line() +
  facet_grid(rows = vars(Series), scales = "free_y") +
  labs(y = "")
オーストラリア・ビクトリア州での2012--2014年の30分ごとの電力需要と対応する気温

図 12.5: オーストラリア・ビクトリア州での2012–2014年の30分ごとの電力需要と対応する気温

電力需要を気温に対してプロットする(図12.6)と、需要は気温が低いと(暖房のため)増加し、気温が高いと(冷房のため)増加する、という非線形の関係が両者の間にあることが見えてきます。

elec <- vic_elec %>%
  mutate(
    DOW = wday(Date, label = TRUE),
    Working_Day = !Holiday & !(DOW %in% c("Sat", "Sun")),
    Cooling = pmax(Temperature, 18)
  )
elec %>%
  ggplot(aes(x=Temperature, y=Demand, col=Working_Day)) +
  geom_point(alpha = 0.6) +
  labs(x="気温 (摂氏度)", y="需要 (MWh)")
ビクトリア州の30分ごとの電力需要とビクトリア州最大都市メルボルンの同時期の気温の散布図

図 12.6: ビクトリア州の30分ごとの電力需要とビクトリア州最大都市メルボルンの同時期の気温の散布図

気温の区分線形関数(18度にノットを1つ含む)と日次季節パターンに対応するハーモニック回帰項を有する回帰モデルを適合させます。ここでも、フーリエ項の次数は主観的に選択する一方、ARIMA誤差の次数はAICcを使って選択させます。

fit <- elec %>%
  model(
    ARIMA(Demand ~ PDQ(0, 0, 0) + pdq(d = 0) +
          Temperature + Cooling + Working_Day +
          fourier(period = "day", K = 10) +
          fourier(period = "week", K = 5) +
          fourier(period = "year", K = 3))
  )

予測変数の将来値が必要なので、このモデルでの予測は困難です。フーリエ項の将来値は計算が容易ですが、将来の気温はもちろん未知です。一週間先までの予測だけに興味があるなら、気候モデルから得た気温予測を使うこともできるでしょう。替わりに、シナリオに基づく予測(6.5節)を使って、あり得る気温パターンを入れ込むこともできるでしょう。以下の例では、最後の2日の気温が繰り返すとして、将来あり得る需要値を生成します。

elec_newdata <- new_data(elec, 2*48) %>%
  mutate(
    Temperature = tail(elec$Temperature, 2 * 48),
    Date = lubridate::as_date(Time),
    DOW = wday(Date, label = TRUE),
    Working_Day = (Date != "2015-01-01") &
                   !(DOW %in% c("Sat", "Sun")),
    Cooling = pmax(Temperature, 18)
  )
fc <- fit %>%
  forecast(new_data = elec_newdata)

fc %>%
  autoplot(elec %>% tail(10 * 48)) +
  labs(title="30分ごとの電力需要: ビクトリア州",
       level = "区間予測",
       y = "需要 (MWh)", x = "時間 [30分]")
動学的ハーモニック回帰モデルを30分ごとの電力需要データに適用して得た予測

図 12.7: 動学的ハーモニック回帰モデルを30分ごとの電力需要データに適用して得た予測

短期予測は理にかなったものに見えますが、これは複雑な過程に対しては粗野なモデルです。図12.8のプロットから分かるように、このモデルに取り込まれなかった多くの情報が残差に出ています。

fit %>% gg_tsresiduals()
動学的ハーモニック回帰モデルの残差診断

図 12.8: 動学的ハーモニック回帰モデルの残差診断

このモデルよりずっと良い予測を生成するさらに洗練されたバージョンについては、 Hyndman & Fan (2010)Fan & Hyndman (2012) に記述があります。

参考文献

Fan, S., & Hyndman, R. J. (2012). Short-term load forecasting based on a semi-parametric additive model. IEEE Transactions on Power Systems, 27(1), 134–141. [DOI]
Hyndman, R. J., & Fan, S. (2010). Density forecasting for long-term peak electricity demand. IEEE Transactions on Power Systems, 25(2), 1142–1153. [DOI]