RStan:時変係数モデル

# データ生成の設定
seed <- 20250702
set.seed(seed)
T <- 52 # 観測期間：52週間

# --- 真のパラメータを生成 ---

# 1. 基礎売上レベル（時変切片 alpha_t）
#    ランダムウォークで生成
alpha_true <- numeric(T)
alpha_true[1] <- 100 # 初期値
sigma_w_alpha <- 2 # システムノイズの標準偏差
for (t in 2:T) {
  alpha_true[t] <- alpha_true[t - 1] + rnorm(1, mean = 0, sd = sigma_w_alpha)
}

# 2. 広告効果（時変係数 beta_t）
#    シナリオに沿って生成
beta_true <- numeric(T)
beta_true[1] <- 2.5 # 初期値
sigma_w_beta <- 0.1 # システムノイズの標準偏差
for (t in 2:T) {
  beta_true[t] <- beta_true[t - 1] + rnorm(1, mean = 0, sd = sigma_w_beta)
}
# 摩耗効果を表現するために、全体的に少しずつ減少させる
beta_true <- beta_true - seq(0, 1, length.out = T)
# 30週目に競合が登場し、効果が急落するイベント
beta_true[30:T] <- beta_true[30:T] - 1.0

# 3. 説明変数：広告投下量（GRP）
#    毎週のGRPをランダムに生成
X_grp <- round(runif(T, min = 50, max = 150))

# 4. 観測ノイズ
sigma_v <- 20 # 観測ノイズの標準偏差

# 5. 観測データ（売上 Y_sales）の生成
#    観測方程式: Y_t = alpha_t + beta_t * X_t + v_t
Y_sales <- numeric(T)
for (t in 1:T) {
  Y_sales[t] <- alpha_true[t] + beta_true[t] * X_grp[t] + rnorm(1, mean = 0, sd = sigma_v)
}
Y_sales <- round(Y_sales) # 売上台数なので整数に丸める

# --- 生成したデータをデータフレームにまとめる ---
sim_data <- tidyr::tibble(
  time = 1:T,
  sales = Y_sales,
  grp = X_grp,
  alpha_true = alpha_true,
  beta_true = beta_true
)

# --- 生成したデータの可視化 ---
library(ggplot2)
p1 <- ggplot(sim_data, aes(x = time, y = sales)) +
  geom_line() +
  labs(title = "週次売上台数の推移", x = "週", y = "売上台数") +
  theme_bw()

p2 <- ggplot(sim_data, aes(x = time, y = beta_true)) +
  geom_line(color = "red", linewidth = 1) +
  geom_vline(xintercept = 30, linetype = "dashed") +
  annotate("text", x = 30.5, y = 1.2, label = "競合製品の登場", hjust = 0, size = 4) +
  labs(title = "真の広告効果（beta）の推移", x = "週", y = "広告効果（係数）") +
  theme_bw()

# patchworkでグラフを並べる
library(patchwork)
p1 + p2

# stanfit オブジェクトの読み込み
setwd(stan_output)
fit <- readRDS("stan_fit.rds")

# MCMCサンプルを抽出
mcmc_samples <- rstan::extract(fit)

# 推定結果をデータフレームにまとめる
library(dplyr)
result_df <- tidyr::tibble(
  time = 1:T,
  # alpha（時変切片）の事後統計量
  alpha_mean = apply(mcmc_samples$alpha, 2, mean),
  alpha_median = apply(mcmc_samples$alpha, 2, median),
  alpha_l95 = apply(mcmc_samples$alpha, 2, quantile, probs = 0.025),
  alpha_u95 = apply(mcmc_samples$alpha, 2, quantile, probs = 0.975),
  # beta（時変係数）の事後統計量
  beta_mean = apply(mcmc_samples$beta, 2, mean),
  beta_median = apply(mcmc_samples$beta, 2, median),
  beta_l95 = apply(mcmc_samples$beta, 2, quantile, probs = 0.025),
  beta_u95 = apply(mcmc_samples$beta, 2, quantile, probs = 0.975)
) %>%
  # シミュレーションの真の値を結合
  left_join(sim_data, by = "time")


# --- 結果のプロット ---

# 1. 時変切片 alpha のプロット
p_alpha <- ggplot(result_df, aes(x = time)) +
  geom_ribbon(aes(ymin = alpha_l95, ymax = alpha_u95), fill = "skyblue", alpha = 0.5) +
  geom_line(aes(y = alpha_median), color = "dodgerblue", linewidth = 1) +
  geom_line(aes(y = alpha_true), color = "black", linetype = "dashed") +
  labs(
    title = "時変切片 (alpha) の推定結果",
    subtitle = "実線: 推定中央値, 網掛け: 95%信用区間, 破線: 真値",
    x = "週", y = "基礎売上レベル"
  ) +
  theme_bw() +
  theme(plot.subtitle = element_text(size = 9))

# 2. 時変係数 beta のプロット
p_beta <- ggplot(result_df, aes(x = time)) +
  geom_ribbon(aes(ymin = beta_l95, ymax = beta_u95), fill = "salmon", alpha = 0.5) +
  geom_line(aes(y = beta_median), color = "firebrick", linewidth = 1) +
  geom_line(aes(y = beta_true), color = "black", linetype = "dashed") +
  geom_vline(xintercept = 30, linetype = "dashed") +
  annotate("text", x = 30.5, y = 1.2, label = "競合製品の登場", hjust = 0, size = 4) +
  labs(
    title = "時変係数 (beta) の推定結果",
    subtitle = "実線: 推定中央値, 網掛け: 95%信用区間, 破線: 真値",
    x = "週", y = "広告効果（係数）"
  ) +
  theme_bw() +
  theme(plot.subtitle = element_text(size = 9))

p_alpha + p_beta