【時系列】広告効果をモデリングする ~その7 周期性をモデル化~ [時系列解析 / 需要予測]
ARIMAモデルとかだと、元系列をトレンド+季節成分+不規則変動に分けたりします。
状態空間モデルも同じように周期性(季節成分)をモデル化することができます。
ポイントは、季節成分の部分。
x_seasonal[i] ~ normal(-sum(x_seasonal[(i-11):(i-1)]), sd_sea);
例えば、月別データ(周期が12か月)の場合は、
a1 + a2 + ... + a12 = ε, ε~N(0, σ2)
と周期成分の合計が0になるように書くことができ、
a1 = -Σ(a2 + ... + a12)
となっています。
■ rのコード
■ stanファイル
状態空間モデルも同じように周期性(季節成分)をモデル化することができます。
ポイントは、季節成分の部分。
x_seasonal[i] ~ normal(-sum(x_seasonal[(i-11):(i-1)]), sd_sea);
例えば、月別データ(周期が12か月)の場合は、
a1 + a2 + ... + a12 = ε, ε~N(0, σ2)
と周期成分の合計が0になるように書くことができ、
a1 = -Σ(a2 + ... + a12)
となっています。
■ rのコード
AirPassengers <- c(112, 118, 132, 129, 121, 135, 148, 148, 136, 119, 104, 118,
115, 126, 141, 135, 125, 149, 170, 170, 158, 133, 114, 140,
145, 150, 178, 163, 172, 178, 199, 199, 184, 162, 146, 166,
171, 180, 193, 181, 183, 218, 230, 242, 209, 191, 172, 194,
196, 196, 236, 235, 229, 243, 264, 272, 237, 211, 180, 201,
204, 188, 235, 227, 234, 264, 302, 293, 259, 229, 203, 229,
242, 233, 267, 269, 270, 315, 364, 347, 312, 274, 237, 278,
284, 277, 317, 313, 318, 374, 413, 405, 355, 306, 271, 306,
315, 301, 356, 348, 355, 422, 465, 467, 404, 347, 305, 336,
340, 318, 362, 348, 363, 435, 491, 505, 404, 359, 310, 337,
360, 342, 406, 396, 420, 472, 548, 559, 463, 407, 362, 405,
417, 391, 419, 461, 472, 535, 622, 606, 508, 461, 390, 432)
head(AirPassengers, 10)
ts.plot(log(AirPassengers), type="o", ylab="AirPassengers")
# データの準備
data_list <- list(
y = log(AirPassengers),
T = length(AirPassengers)
)
# モデルの推定
seasonal_model <- stan(
file = "D:\\Projects\\LM\\arima\\seasonal_model_1.stan",
data = data_list,
chains = 3,
seed = 1234
)
■ stanファイル
data {
int T; // データ取得期間の長さ
vector[T] y; // 観測値
}
parameters {
vector[T] mu; // 水準成分の推定値(過程誤差)
vector[T] x_seasonal; // 季節成分の推定値(過程誤差)
real sd_mu; // 水準成分の標準偏差(過程誤差)
real sd_sea; // 季節成分の標準偏差(過程誤差)
real sd_v; // 観測誤差の標準偏差
}
transformed parameters {
vector[T] alpha;
for (i in 1:T) {
alpha[i] = mu[i] + x_seasonal[i];
}
}
model {
// 水準成分
for(i in 3:T) {
mu[i] ~ normal(2*mu[i-1] - mu[i-2], sd_mu);
}
// 季節成分
for(i in 12:T) {
x_seasonal[i] ~ normal(-sum(x_seasonal[(i-11):(i-1)]), sd_sea);
}
// 観測モデル
for(i in 1:T) {
y[i] ~ normal(alpha[i], sd_v);
}
}