概要

これは、BigQuery Advent Calendar 2023 の 3 日目の記事の仮原稿として作成したものである。

まとめ直した正式版の記事はこちら：

この Scrap は？

BigQuery ML の ARIMA_PLUS は、ARIMA モデルをベースにしたアルゴリズムで時系列のモデリングを自動で行う。

しかし、そのアルゴリズムの詳細は公開されていないため、
公式ドキュメントとそのパイプライン図を手がかりに推測するしかない。

ARIMA_PLUS のアルゴリズムのうち、トレンド成分と季節成分の分解については、前述のドキュメントとパイプライン図に基づくと、以下のように行なっているものと推測できる（したがって、 SARIMA モデルとは異なるという認識）：

STL によって季節成分を推定する
「1.」で得られた季節成分を差し引いた残差成分から ARIMA モデルを作成し、その予測値をトレンド成分と見做す

もしこの推測が正しければ、季節成分とトレンドの他に祝日効果やスパイク・外れ値のないデータであれば、STL と ARIMA モデルを組み合わせることで R など他の実行環境でも BQML と同様の実行結果を生成できるはずである。

従って、この一連の Scrap では、上記の推測通りにトレンド・季節成分分解を行っているか検証すべく、 R でモデリングを行い、実行結果が BQML の ARIMA_PLUS によるものと一致するかを検証した。

Summary

~~R による STL と ARIMA を組み合わせたモデリングの結果は、BQML ARIMA_PLUS による結果と一致せず、季節成分・トレンド成分共に異なる値が出力された。~~

~~不一致の原因としては、前処理の違いや、STL 時のパラメータの与え方の違いが考えられる。~~

※2023/12/10追記
上記の正式版の記事を参照ください。

実行コード (R)

実験に用いた R コードはこちら

畳屋民也 5ヶ月前に更新

サンプルデータの生成

y = （トレンド成分）+（周期成分）+（ノイズ）

y：データ点が100個の時系列
トレンド成分：次数(1,1,2)のARIMAモデル
周期成分：周期7の sine 曲線

set.seed(13)

tr <- 1.2 * arima.sim(model=list(order=c(1,1,2), ar=c(0.8), ma=c(0.2, -0.8)), n=99)
s <- 10 * sin(seq(0, 99, length.out=100) * 2 * pi / 7)
y <- tr  + s + rnorm(100)

time.points <- seq.Date(as.Date("2014-01-01"), by=1, length.out=100)

data <- data.frame(day=time.points, y=y, trend=tr, seasonal=s)

p1 <- ggplot(data=data) +
  geom_line(aes(x=day, y=y))+
  labs(x="day", y="y")

p2 <- ggplot(data=data) +
  geom_line(aes(x=day, y=trend))+
  labs(x="day", y="trend")

p3 <- ggplot(data=data) +
  geom_line(aes(x=day, y=seasonal))+
  labs(x="day", y="seasonal")


grid.arrange(p1, p2, p3)

畳屋民也 5ヶ月前

BQML による ARIMA_PLUS モデル

モデル作成

CREATE OR REPLACE MODEL `project_id.dataset.arima_plus_exploration`
 OPTIONS(MODEL_TYPE='ARIMA_PLUS',
         time_series_timestamp_col='day',
         time_series_data_col='y',
         DECOMPOSE_TIME_SERIES=TRUE
         ) AS
SELECT
  day,
  y
FROM
  `project_id.dataset.sample_for_arima_plus`

作成結果

係数

4次の MA 過程として推定された。

SELECT
  *
FROM
  ML.ARIMA_COEFFICIENTS(MODEL `project_id.dataset.arima_plus_exploration` )

ar_coefficients	ma_coefficients	intercept_or_drift
	0.183275788118965	0.0
	0.44292870561309977
	-0.48301393033169393
	-0.42382755350887391

モデル説明（成分分解）

weekly の周期成分以外に holiday_effect, spikes_and_dips, step_changes などは持たない：

SELECT
  *
FROM
  ML.EXPLAIN_FORECAST(MODEL `project_id.dataset.arima_plus_exploration`,
                      STRUCT(0.95 AS confidence_level))
WHERE time_series_type="history"

可視化

Looker Studio 使用

観測値と予測値（全成分）

季節成分

入力データ生成時に与えた真の季節成分は（同じ曜日なら）一定であったが、BQML では時間変化して推定されている。

トレンド成分

入力データ生成時に与えた真のトレンドと比べて、なめらかでない。

Hidden comment

畳屋民也 5ヶ月前に更新

2023/12/09追記： R によるモデル作成（STL で季節成分を固定にしなかった場合）

前回、R によるモデル作成時に STL を stl_results <- stl(data_ts, s.window="per") のように実行していた。

これにより、季節成分は（同じ曜日なら）時間変化せず一定になるように成分推定が行われていた。

しかし、s.windows=7 とすることで時間変化を許したところ、季節成分に関しては BQML による出力結果とほぼ一致する結果が得られた。

STL による季節成分検出

data_ts <- ts(data$y, frequency = 7)
stl_results <- stl(data_ts, s.window=7)

data_stl_results <- data %>%
  mutate(
  seasonal_stl=stl_results$time.series[,"seasonal"],
  trend_stl=stl_results$time.series[,"trend"],
  residual_stl=stl_results$time.series[,"remainder"]
)


p_s <- ggplot(data=data_stl_results) +
  geom_line(aes(x=day, y=seasonal, color="actual", linetype="actual"))+
  geom_line(aes(x=day, y=seasonal_stl, color="STL", linetype="STL"))+
  scale_color_manual(name="", values = c("STL"="blue", "actual"="black")) +
  scale_linetype_manual(name="", values = c("STL"="solid", "actual"="dashed")) +
  labs(x="day", y="seasonal")

p_t <- ggplot(data=data_stl_results) +
  geom_line(aes(x=day, y=trend, color="actual", linetype="actual"))+
  geom_line(aes(x=day, y=trend_stl, color="STL", linetype="STL"))+
  scale_color_manual(name="", values = c("STL"="blue", "actual"="black")) +
  scale_linetype_manual(name="", values = c("STL"="solid", "actual"="dashed")) +
  labs(x="day", y="trend")

grid.arrange(p_s, p_t)

ARIMA によるトレンド成分の検出＆モデル化

MA(4) 係数の値は、1.834006, 2.117289, 1.834006, 0.999999 であった（BQML AUTO_ARIMA と一致せず）。

arima_fit <- arima(data_after_stl$y2, order=c(0,0,4), include.mean = FALSE)
coef(arima_fit)
#       ma1       ma2       ma3       ma4 
# 1.834006 2.117289 1.834006 0.999999 

arima_results <- data.frame(day=time_points,
                            residuals=arima_fit$residuals,
                            fitted=data_after_stl$y2 - arima_fit$residuals,
                            actual_trend=data_after_stl$actual_trend
                            )

p <- ggplot(data=arima_results) +
  geom_line(aes(x=day, y=actual_trend, color="actual trend", linetype="actual trend"))+
  geom_line(aes(x=day, y=fitted, color="fitted", linetype="fitted"))+
  scale_color_manual(name="", values = c("fitted"="blue", "actual trend"="black")) +
  scale_linetype_manual(name="", values = c("fitted"="solid", "actual trend"="dashed")) +
  labs(x="day", y="y2")
print(p)

BQML AUTO_ARIMA との比較

季節成分

見た目、ほぼ一致。

p_s <- ggplot(data=data_all) +
  geom_line(aes(x=day, y=seasonal, color="actual", linetype="actual"))+
  geom_line(aes(x=day, y=seasonal_stl, color="STL", linetype="STL"))+
  geom_line(aes(x=day, y=seasonal_bqml, color="BQML", linetype="BQML"))+
  scale_color_manual(name="", values = c("BQML"="blue", "STL"="black", "actual"="black")) +
  scale_linetype_manual(name="", values = c("BQML"="solid", "STL"="solid","actual"="dotted")) +
  labs(x="day", y="seasonal")

print(p_s)

平均絶対誤差で見ても、7.295137 $\times 10^{-15}$ と極めて差が小さい。

mean((data_all %>% mutate(seasonal_diff=abs(seasonal_bqml-seasonal_stl)))$seasonal_diff)
# 7.295137e-15

トレンド成分

p_t <- ggplot(data=data_all) +
  geom_line(aes(x=day, y=trend, color="actual", linetype="actual"))+
  geom_line(aes(x=day, y=trend_arima, color="arima", linetype="arima"))+
  geom_line(aes(x=day, y=trend_bqml, color="BQML", linetype="BQML"))+
  scale_color_manual(name="", values = c("BQML"="blue", "arima"="black", "actual"="black")) +
  scale_linetype_manual(name="", values = c("BQML"="solid", "arima"="solid","actual"="dotted")) +
  labs(x="day", y="trend")

print(p_t)

考察：なぜトレンド成分が一致しないか

考えられる要因：

ARIMA の推定アルゴリズムが異なる
- 最適化方法が異なる、など。例えば誤差項分散はニュートン法や勾配降下法などを用いて推定する必要がある。このアルゴリズムに差があるために、R と BQML で異なる推定結果となった可能性が考えられる。