【時系列】広告効果をモデリングする ～その4 ローカルレベルモデル～

久々に時系列の話。
広告効果をモデリングする手法は、簡単に線形回帰モデルなどを使って対処する方法もあるのですが、何かと問題がありました。

そこで王道となる状態空間モデルを使って広告効果をモデリングしてみようと思います。

まずは、ツールの話から。

■ dlm、KFASパッケージを利用
線形モデルで、誤差に正規ホワイトノイズを仮定する場合に使うことができます。
動的線形モデル(Dynamic Linear Model)と呼ばれています。
⇒ 最尤推定＋カルマンフィルタによる方法
1）最尤法によるモデルパラメータの推定
2）カルマンフィルタによる予測

■ Stanを利用
非線形なモデルや、正規分布以外の分布を扱う場合に使うことができます。
⇒ ベイズ推定（MCMC法）による方法
MCMCにより、モデルパラメータ分布と予測分布を推定

個人的には、ちょっとstanのコードを書くのが面倒なのですが、覚えてしまえば色々なモデルを解くことができるstanがおススメです。

状態空間モデルは様々なモデルがあるのですが、まずは、一番シンプルなローカルレベルモデルを適用するとどうなるか見てみたいと思います。

ローカルレベルモデルは、ランダムウォーク・プラス・ノイズとも呼ばれ、トレンドや季節性が見られない系列に用いることができます。



■ local-level_20201029.r

library(rstan)

# データの準備
Sales <- c(19.76, 17.68 ,17.56, 21.41, 21.04,
           21.79, 24.08, 24.84, 18.57, 12.62,
           31.72, 27.01, 26.48, 31.62, 27.79,
           23.97, 30.20, 29.56, 28.87, 31.84,
           19.65, 22.02, 20.28, 22.34, 24.79,
           22.74, 20.29, 24.02, 22.96, 21.26,
           19.83, 17.46, 22.49, 19.57, 28.46,
           34.42, 28.80, 26.21, 28.29, 30.90,
           29.58, 29.09, 27.17, 21.24, 18.11,
           16.74, 18.18, 23.00, 18.52, 22.16,
           15.70, 19.61)
CM <- c(0, 0, 0, 0, 0,
        0, 0, 0, 0, 0,
        1, 1, 1, 1, 1,
        1, 1, 1, 1, 1,
        0, 0, 0, 0, 0,
        0, 0, 0, 0, 0,
        0, 0, 0, 0, 0,
        1, 1, 1, 1, 1,
        1, 1, 1, 0, 0,
        0, 0, 0, 0, 0,
        0, 0)

dat <- data.frame(Sales, CM)
head(dat, 10)

# CMあり/なし期間別の統計量
by(dat$Sales, INDICES=dat$CM, FUN=summary)
ts.plot(dat$Sales, type="o", ylab="Sales")

# データの準備
data_list <- list(
  y = dat$Sales, 
  T = nrow(dat)
)

# モデルの推定
local_level_stan <- stan(
  file = "D:\\Projects\\LM\\arima\\local-level_20201029.stan",
  data = data_list,
  seed = 123
)

## 各chainの平均値
get_posterior_mean(local_level_stan)

## Trace Plot（収束の確認）
stan_trace(local_level_stan, inc_warmup=TRUE)

## サンプリングしたパラメータのデータのヒストグラム
stan_dens(local_level_stan, pars="mu")
stan_plot(local_level_stan, pars="mu")

## サンプリング結果
summary(rstan::extract(local_level_stan)$mu)

# 結果の表示
print(local_level_stan,
      pars = c("s_w", "s_v","lp__"),
      probs = c(0.025, 0.5, 0.975))

# サンプリング結果を保存
mcmc_sample <- rstan::extract(local_level_stan)

# 95%ベイズ信用区間と中央値を得る
state_name <- "mu"
quantile(mcmc_sample[[state_name]][, 1], 
         probs=c(0.025, 0.5, 0.975))

# 95%ベイズ信用区間と中央値
result_df <- data.frame(t(apply(
  X = mcmc_sample[[state_name]],# 実行対象となるデータ
  MARGIN = 2,                   # 列を対象としてループ
  FUN = quantile,               # 実行対象となる関数
  probs=c(0.025, 0.5, 0.975)    # 上記関数に入れる引数
)))
colnames(result_df) <- c("lwr", "fit", "upr")

# 結果のプロット
ts.plot(cbind(dat$Sales, result_df$fit), 
        type="o",
        col=c("black","red"),
        xlab="Week", ylab="Sales")

■ local-level_20201029.stan

data {
  int T;         // データ取得期間の長さ
  vector[T] y;   // 観測値
}

parameters {
  vector[T] mu;       // 状態の推定値(水準成分)
  real s_w;  // 過程誤差の標準偏差
  real s_v;  // 観測誤差の標準偏差
}

model {
  // 状態方程式に従い、状態が遷移する
  for(i in 2:T) {
    mu[i] ~ normal(mu[i-1], s_w);
  }
  
  // 観測方程式に従い、観測値が得られる
  for(i in 1:T) {
    y[i] ~ normal(mu[i], s_v);
  }
}

■ 結果のプロット

結果を見ると、 CMあり/なしの変数が入っていないので、CMありの部分で水準が増えており、CMなしの部分で水準が低下しているような結果となっています。
知りたいのはCMの効果なのですが、ローカルレベルモデルだとイマイチしっくりこない結果となってしまいます。

そこで、次回はCMあり/なしを外生変数（説明変数）に入れてモデルを作ってみたらどうなるかやってみたいと思います。

～過去の記事～
【時系列】広告効果をモデリングする

その1 ABテスト（t検定）
https://skellington.blog.ss-blog.jp/2020-07-20

その2 簡単な線形回帰モデル
https://skellington.blog.ss-blog.jp/2020-07-21

その3 線形回帰モデルにラグを入れる
https://skellington.blog.ss-blog.jp/2020-07-22