はじめに

こんにちは、事業会社で働くデータサイエンティストです。

先日、業務の一環で時系列データの因果関係を調べていました。分析には、時系列データ間の因果関係を推定する手法の一つにインパルス応答関数を用いました。

インパルス応答関数やその実装方法に関するネット記事は多いものの、個人的に一つ気になった点があります。それは、インパルス応答関数を使う際の変数の並べ方に触れている記事が（管見の限り）見当たらなかったことです。

インパルス応答関数は変数の並べ方次第で推定結果が変わってしまうため、ドメイン知識に基づいた妥当な順序で変数を並べる必要があると考えています。そこで、今回はその理由を述べた上で、簡単な実装を行いました。

インパルス応答関数の概要

ベクトル自己回帰モデル

インパルス応答関数について述べる前に、ベクトル自己回帰モデル（vector autoregressive, 以下VARと表記）について簡単に説明しておく。

VARは複数の時系列データの動きを同時に分析する枠組みである。
例えば、1次の自己相関を持つ2変数のモデルを考える。以下、西山計量にならって記述する。

$Y_t = \beta_{10} + \beta_{11}Y_{t-1} + \beta_{12}X_{t-1} + \epsilon_{1t} \quad \cdots ①$

$X_t = \beta_{10} + \beta_{11}Y_{t-1} + \beta_{12}X_{t-1} + \epsilon_{2t} \quad \cdots ②$

$\epsilon_{1t}$と$\epsilon_{2t}$が独立ならば、以下のようにして今期の$X_t$から将来の$Y_t$への処置効果を評価することができる。

$E[Y_{t+1}\mid \epsilon_{2t} = 1] - E[Y_{t+1}\mid \epsilon_{2t} = 0] = \beta_{12}$

だが、実際は$\epsilon_{1t}$と$\epsilon_{2t}$が相関している可能性が高い。そこで、政策や施策による純粋なショックのみを誤差項から取り出すのがインパルス応答関数である。

インパルス応答関数

では、実際に誤差項から純粋なショックを取り出す手続きを考える。ここで扱うのは、短期制約と呼ばれる仮定である。

誤差項$\epsilon_{1t}$と$\epsilon_{2t}$を以下のように分解する。

$\epsilon_{1t} = c_{11}e_{1t} + c_{12}e_{2t}, \quad \epsilon_{1t} \sim N(0,\sigma_1^2) \quad \cdots ③$
$\epsilon_{2t} = c_{21}e_{1t} + c_{22}e_{2t}, \quad \epsilon_{2t} \sim N(0,\sigma_2^2) \quad \cdots ④$

ここで$e_{1t}$と$e_{2t}$は互いに相関しない純粋なショックである。

誤差項$\epsilon_{1t}$と$\epsilon_{2t}$の共分散が$\sigma_{12}$で既知のとき、$c_{11}$, $c_{12}$, $c_{21}$, $c_{22}$が分かれば③と④を$e_{1t}$と$e_{2t}$について解くことで、純粋なショックを抽出できる。

いま、以下の3本の式を考えることができる。
$\sigma_1^2 = V(c_{11}e_{1t} + c_{12}e_{1t}) = c_{11}^2 + c_{12}^2$
$\sigma_2^2 = V(c_{21}e_{2t} + c_{22}e_{2t}) = c_{21}^2 + c_{22}^2$
$\sigma_{12} = Cov(c_{11}e_{1t} + c_{12}e_{1t}, ~ c_{21}e_{2t} + c_{22}e_{2t}) = c_{11}c_{21} + c_{12}c_{22}$

上記の3本の式に対して、求めたい係数は$c_{11}$, $c_{12}$, $c_{21}$, $c_{22}$の4つなので、識別するには係数が1つ多い。

そこで、$c_{21}=0$という制約を課すのが短期制約の考え方である。

①と②に③と④を代入する。さらに、①と②に過去の$Y_t$,$X_t$を逐次代入すれば、今期の$Y_t$と$X_t$は純粋なショックの加重無限和で表現できる。

$Y_t = \sum_{j =0}^\infty c_{11}(j)e_{1t-j} ~ + \sum_{j =0}^\infty c_{12}(j)e_{2t-j} \quad \cdots ①^{'}$

$X_t = \sum_{j =0}^\infty c_{21}(j)e_{1t-j} ~ + \sum_{j =0}^\infty c_{22}(j)e_{2t-j} \quad \cdots ②^{'}$

いま、$c_{21}=0$という制約を課しているので、

$②^{'} \Leftrightarrow X_t = \sum_{j =0}^\infty c_{22}(j)e_{2t-j} \quad \cdots ②^{''}$

つまり、純粋なショック$e_{1t}$は今期の$Y_t$に影響を与えるが、今期の$X_t$には影響を与えない。逆に言えば、①と②は外生性の高い順に並んでおかなければならない、ということである。

実際にインパルス応答関数を推定する際は、変数の並べ方に注意する必要がある。インパルス応答関数は時系列の因果関係を調べるのに便利だが、知りたい因果関係の順序を仮定しなければならないというジレンマがある。

Rで実装

TVCMの費用が何期先の売上まで影響しているかが気になるような場面を考える。データを生成するプロセスは「因果推論基礎から機械学習・時系列分析・因果探索を用いた意思決定のアプローチ」を参考にした。本記事においてデータを生成する過程は本質的でないので、読み飛ばしていただいて構わない。データを生成した上で、インパルス応答関数を推定する。

データの生成

売上は以下の回帰式に従って生成されるとする。

$Y_t = \beta_0 + \sum\beta_1\times Hill(T_t^*;K,S) + Trend_t + Seasonality_t + \epsilon_t$

ただし、$Hill(T_t^*;K,S)$はヒル関数で、

$Hill(T_t;K,S)=\frac{1} {1 + (\frac {T_t}{K})^{-S}},$

$where \quad T_t\geq0.$

また、$T_t^*$はアドストック関数で、

$T_t^* = adstock(T_{t-L+1}, \cdots, T_t;w,L) = \frac {\sum_{l = 0}^{L-1}w(l)\times T_{t-l}} {\sum_{l=0}^{L-1}w(l)}.$

$w(l)$は重み関数で、

$w(l; \alpha, \theta) = \alpha^{(l- \theta_m)^2},$

$where \quad l = 0,\cdots L-1, ~ 0 < \alpha < 1, ~ 0 \leq \theta \leq L-1.$

重み関数を以下のように定義しておく。

#重み関数を定義

weight_function <- function(l, alpha = 0.2, theta = 0) {
  return(alpha^((l - theta)^2))
}

アドストック関数の分母は以下のようになる。

#アドストック関数の分母を定義しておく
denominator <- function(alpha, theta, L) {
  l <- 0:L-1
  total_weight <- sum(alpha^((l - theta)^2))
  return(total_weight)
}

では、実際に売上データを生成していこう。

df_sales <- data.frame(date = seq(from = as.Date("2020-04-01"),
                                  to = as.Date("2023-09-01"),
                                  by = "week")) |>
  dplyr::mutate(id = seq(from = 1,
                         to = dplyr::n())) |>
  dplyr::mutate(tv = runif(n = dplyr::n(), min = 0, max = 1),
                tv = dplyr::case_when(
                  tv > 0.9 ~ tv,
                  TRUE ~ tv/2
                )) |>
  dplyr::mutate(tv1 = dplyr::lag(tv,1),
                tv2 = dplyr::lag(tv,2),
                tv3 = dplyr::lag(tv,3),
                tv4 = dplyr::lag(tv,4),
                tv5 = dplyr::lag(tv,5),
                tv6 = dplyr::lag(tv,6),
                tv7 = dplyr::lag(tv,7)) |>
  dplyr::mutate(dplyr::across(c(tv, tv1, tv2, tv3, tv4, tv5, tv6, tv7), 
                              ~tidyr::replace_na(., 0))) |>
  dplyr::mutate(adstock_tv_numerator = tv * weight_function(l = 0) +
                  tv1 * weight_function(l = 1) +
                  tv2 * weight_function(l = 2) +
                  tv3 * weight_function(l = 3) +
                  tv4 * weight_function(l = 4) +
                  tv5 * weight_function(l = 5) +
                  tv6 * weight_function(l = 6) +
                  tv7 * weight_function(l = 7)) |>
  dplyr::mutate(adstock_tv_denominator = denominator(alpha = 0.2, theta = 0, L = 8)) |>
  dplyr::mutate(adstock_tv = adstock_tv_numerator/adstock_tv_denominator) |>
  dplyr::mutate(hill_tv = 1/(1 + ((adstock_tv/0.3)^(-0.5)))) |>
  dplyr::mutate(day_of_year = lubridate::yday(date),
                seasonality = 0.5*(-sin(2*2*pi*day_of_year/365.5) + 
                                     cos(1*2*pi*day_of_year/365.5))) |>
  dplyr::mutate(trend0 = seq(from = 0,
                             to = 50,
                             length.out = dplyr::n()),
                trend = (trend0 + 10)^(1/4) - 1) |>
  dplyr::mutate(intercept = 2.0) |>
  dplyr::mutate(beta_tv = 3.0) |>
  dplyr::mutate(epsilon = rnorm(n = dplyr::n(), mean = 0, sd = 0.25)) |>
  dplyr::mutate(sales = intercept + trend + seasonality + beta_tv * hill_tv + epsilon)

データの可視化

売上データを可視化する。

ggplot2::ggplot(df_sales, ggplot2::aes(x = date, y = sales)) +
  ggplot2::geom_line() +
  ggplot2::scale_x_date(date_labels = "%y-%m", date_breaks = "1 month") +
  ggplot2::theme_bw() +
  ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust = 0.8))

ここで売上額の単位は重要でないので、とくに考慮していない。

インパルス応答関数の推定

本来は定常性の仮定を確認することが望ましい（変数が階差変換によってはじめて定常性を満たす可能性があるから）が、ここでは割愛する。

VARモデルを準備するために、データを整形する。

df_sales_var <- df_sales |>
  dplyr::select(c(tv, sales))

インパルス応答関数は、{vars}パッケージを用いて簡単に推定することができる。最適なラグ数を選択し、VARを推定。それをirf関数に投げ込む。

#最適なラグ数の選択
var_select <- vars::VARselect(df_sales_var, lag.max = 10, type = "both", season = 4)

#VARの推定
sales_var <- vars::VAR(df_sales_var, p = var_select$selection[1], lag.max = 10, type = "both", season = 4)

#インパルス応答関数の推定
sales_irf <- vars::irf(sales_var, n.ahead = 8)

インパルス応答関数の可視化

{ggplot2}パッケージを用いて可視化した。

data.frame(impulse = as.vector(sales_irf$irf$tv),
           upper = as.vector(sales_irf$Upper$tv),
           lower = as.vector(sales_irf$Lower$tv),
           to = rep(c("tv", "sales"), c(9,9))) |>
  ggplot2::ggplot(ggplot2::aes(x = rep(1:9, 2), y = impulse)) +
  ggplot2::geom_line(ggplot2::aes(color = to)) +
  ggplot2::geom_ribbon(ggplot2::aes(ymax = upper, ymin = lower), alpha = 0.2, linetype = 2) +
  ggplot2::facet_wrap(.~ to, nrow = 2) +
  ggplot2::labs(title = "TV→売上の影響", x = "lag") +
  ggplot2::scale_x_continuous(breaks = seq(2, 8, by = 2)) +
  ggplot2::theme_bw()

おわりに

残念ながら、今回生成したデータではTVCMに与えられたショックが1期以上先の売上に影響を与えているかどうかは判断できませんでした（信頼区間が0をまたいでいるので）。

ですが、本記事で重要なのは変数の並べ方です。今回はTVCMの費用の方が売上より外生性が高いという判断のもと、インパルス応答関数を推定しました。ビジネスの現場ではより多くのKPIを分析に投入しなければならないでしょう。その際、できるだけドメイン知識に基づいた妥当な変数の並べ方を探る必要があります。

参考文献

沖本竜義. 2010. 「経済・ファイナンスデータの計量時系列分析」. 朝倉書店.
金本拓. 2024. 「因果推論基礎から機械学習・時系列分析・因果探索を用いた意思決定のアプローチ」. オーム社.
西山慶彦ら. 2019. 「計量経済学」. 有斐閣.