1. はじめに
効果検証において、ユーザーごとや店舗ごとの詳細な履歴を追跡したパネルデータが利用可能だとします。このとき、差の差法（DID）や合成コントロール法、固定効果モデルといった因果推論の枠組みに落とし込むことでマーケティング施策の因果効果を推定できる可能性が高まります。
しかし、実際のマーケティング実務においては、テレビCMや交通広告（OOH）といったマス媒体で誰がいつ広告に接触したかをユーザー単位でトラッキングできません。したがって「週ごとの総出稿量（GRP）」や「日ごとの広告費」というマクロな時系列データとしてしか観測できないという限界があります。
このように、マーケティングのデータ集計が日次や週次の単一の時系列データに丸められてしまうと、パネルデータで使えた個体・ユニットごとの違いを利用した因果推論的アプローチが使えなくなります。
時系列データの世界では、全ての要因（競合の動き、季節性、自社の複数媒体の同時出稿）が時間軸上で複雑に絡み合って同時に変化します。だからこそ、パネルデータ的な因果推論が使えない環境において、変数同士が互いに与え合う動的なインパクトを捉えるVARモデル（ベクトル自己回帰）のような多変量時系列解析のアプローチが実務における主要な選択肢となると考えます。
2. なぜ通常の回帰分析（OLS）では不十分なのか？
マーケティング効果の測定において、一般的な最小二乗法（OLS）による回帰分析（「サイト流入数 ～ 広告費」など）をそのまま時系列データに適用することには、統計学的なリスクが伴います。通常の回帰分析が前提とする「誤差項の独立性（互いに無関係であること）」という仮定から、時系列データは構造的に大きく逸脱するためです。
典型的な統計学の枠組みでは、ランダムサンプル$X_1, X_2, \cdots, X_n$は$\text{i.i.d}(\mu, \sigma^2)$であると想定されます（ただし、$0 <\sigma^2<\infty$）。この仮定のもとでは、いくつかのよい性質が得られます。
-
一致性
$$
\bar X_n=\frac{1}{n}\Sigma_{t=1}^nX_i \xrightarrow{p}\mu=E(X_t) \quad (n\rightarrow \infty)
$$ -
漸近正規性
$$
\sqrt{n}(\bar X_n-\mu)\xrightarrow{d}N(0,\sigma^2) \quad (n\rightarrow \infty)
$$
しかしながら、時系列データは$X_1, X_2, \cdots, X_n$が相互に依存しあうことを想定するため、i.i.d.のフレームワークを適用することは不適切となります。
古典的な線形回帰モデルを考えます。
$$
y = X \beta + u
$$
$y$は$n\times1$,$X$は$n\times k$, $\beta$は$k\times1$, $u$は$n\times1$のベクトルです。
ここで、以下を仮定します。
- $\text{rank}(X) = k \quad (k \leq n)$
- $E(u | X) = 0$
- $E(uu' | X) = \sigma^2I$
これらの仮定のもとで、OLS推定量$\hat \beta$は不偏性を満たします。
$$
E(\hat\beta) = \beta
$$
続いて、線形回帰モデルを動的時系列で置き換えます。
$$
y = \beta y_{-1} + u
$$
ただし、
$y = \begin{pmatrix}
y_2 \\
y_3 \\
\vdots \\y_n
\end{pmatrix}, \quad
y_{-1} = \begin{pmatrix}
y_1 \\
y_2 \\
\vdots \\y_{n-1}
\end{pmatrix}, \quad
u = \begin{pmatrix}
u_2 \\
u_3 \\
\vdots \\u_n
\end{pmatrix}.
$
以下の仮定のもとで、
- $|\beta|<1$
- $\lbrace y_t \rbrace$は定常
$$
y_t = \Sigma_{j=0}^{\infty}\beta^j u_{t-j}.
$$
$y_t$は当期と過去の誤差項に依存するため、古典的な線形回帰モデルが満たしていた不偏性を満たさなくなります。
3. VARモデル
VAR（ベクトル自己回帰）モデルは、自己回帰モデルを多変量に拡張したもので、変数間の動学的関係の分析を行うことができます。
VARモデルはすべての方程式が同一の説明変数をもちます。すべての方程式を同時に推定する必要はなく、各方程式を個別にOLS推定することで、漸近正規性をもった推定量を得ることができます。
3.1 VAR(p)モデル
過去 $p$ 期前までのラグを考慮した$\text{VAR}(p)$ モデルを考えます。
$Y_t = \mu + A_1 Y_{t-1} + \dots + A_p {Y}_{t-p} + \epsilon_t$
- $Y_t$ : $m \times 1$ のベクトル
- $\mu$ : $m \times 1$ のベクトル
- $A_l$ ($l = 1, 2, \dots, p$) : $m \times m$ の行列
- $\epsilon_t \equiv \text{i.i.d.}(0, \Sigma_{\epsilon})$
高次のVARを表記するために、より便利な表現を考えます。
$\tilde Y_t = \begin{pmatrix}
Y_t \\
Y_{t-1} \\
\vdots \\
Y_{t-p+1}
\end{pmatrix}, \quad
\tilde \mu = \begin{pmatrix}
\mu \\
0 \\
\vdots \\
0
\end{pmatrix}, \quad
A = \begin{pmatrix}
A_1 & A_2 & \cdots & A_{p-1} & A_p \\
I_m & 0 & \cdots & 0 & 0 \\
0 & I_m & \ddots & \vdots & 0 \\
\vdots & \ddots & \ddots & 0 & \vdots \\
0 & \cdots & \cdots & I_m & 0 \\
\end{pmatrix}, \quad
\tilde \epsilon_t = \begin{pmatrix}
\epsilon_t \\
0 \\
\vdots \\
0 \\
\end{pmatrix}.
$
ただし、$\tilde Y_t$,$\tilde \mu$,$\tilde \epsilon_t$は$mp\times1$のベクトル、$A$は$mp \times mp$の行列です。
これにより、$\text{VAR}(p)$ モデルは以下のように書き表され、$\text{VAR}(1)$ モデルで表現できることがわかります（companion form）。
$\tilde Y_t = \begin{pmatrix}
Y_t \\
Y_{t-1} \\
\vdots \\
Y_{t-p+1}
\end{pmatrix} =
\begin{pmatrix}
\mu \\
0 \\
\vdots \\
0
\end{pmatrix}
+
\begin{pmatrix}
A_1 & A_2 & \cdots & A_{p-1} & A_p \\
I_m & 0 & \cdots & 0 & 0 \\
0 & I_m & \ddots & \vdots & 0 \\
\vdots & \ddots & \ddots & 0 & \vdots \\
0 & \cdots & \cdots & I_m & 0 \\
\end{pmatrix}
\begin{pmatrix}
Y_{t-1} \\
Y_{t-2} \\
\vdots \\
Y_{t-p}
\end{pmatrix}
+
\begin{pmatrix}
\epsilon_t \\
0 \\
\vdots \\
0
\end{pmatrix}.
$
$\tilde Y_t = \tilde \mu + A \tilde Y_{t-1} + \tilde \epsilon_t$
3.2 VARパラメータのOLS推定量の漸近正規性
以下の2つの仮定を考えます。
1.（定常性）$\text{det}(I_{mp} - A_z) = \text{det}(I_{mp} -A_1z - \cdots - A_pz^p)=0 \Rightarrow |z|>1.$
2. $\lbrace \epsilon_t\rbrace \equiv \text{i.i.d.}(0, \Sigma_{\epsilon}).\quad$ただし、正の定数$C$が存在して、$0 < 1/C <\lambda_{min}(\Sigma_{\epsilon})<\lambda_{max}(\Sigma_{\epsilon})\leq C < \infty.$
これらの仮定のもとで、VARパラメータのOLS推定量は漸近正規性を満たします（証明略）。
$$\sqrt T(\beta_{OLS}-\beta) \xrightarrow{d} N(0, V) \quad(T \rightarrow \infty)$$
ただし、$\beta = \text{vec}B, V = (\Sigma_{\epsilon}\otimes M^{-1})$で、
$M = p \lim_{T\rightarrow \infty}\frac{X'X}{T},$
$
B = \begin{pmatrix}
A_1^{'} \\
A_2^{'} \\
\vdots \\
A_{p}^{'}
\end{pmatrix}
,$
$
X = \begin{pmatrix}
X_1^{'} \\
X_2^{'} \\
\vdots \\
X_{T}^{'}
\end{pmatrix}
,$
$
X_t = \begin{pmatrix}
Y_{t-1} \\
Y_{t-2} \\
\vdots \\
Y_{t-p}
\end{pmatrix}
.$
$B$は$mp\times m$,$X$は$T\times mp$, $X_t$は$mp\times1$のベクトルです。
4. 分析の全体フロー
時系列データを用いたVARモデルの構築は、データの「定常性」と変数間の「長期的な関係（共和分関係）」を正しく評価しながら、以下のフローに沿って進めます。
① ステップ1：単位根（Unit Root）過程のチェック
投入する各時系列データが定常か非定常（単位根過程）かをADF（拡張ディッキー・フラー）検定などで確認します。
- 単位根なし（定常系列）: 原系列をVARモデルに投入します。
- 単位根あり（非定常系列）: 共和分チェックへ進みます。
② ステップ2：共和分（Cointegration）関係のチェック
単位根がある変数同士の中に、長期的な均衡関係が存在するかをヨハンセン（Johansen）検定などで確認します。検定結果のランク（行列の階数）によって分岐します。
- ランク $> 0$（共和分関係あり）: 変数間に長期的なつながりがあるため、単に差分をとってVARに投入すると、その長期的な情報が失われてしまいます。そのため、短期的なズレの修正メカニズムを組み込んだベクトル誤差修正モデル（VECM）を経由してVARモデルへと落とし込みます。
- ランク $= 0$（共和分関係なし）: 長期的な関係はないため、データを1期前の差分（前月比など）に変換して定常系列にした上で、通常のVARモデルに投入します。
※実務上の重要な注意点
2〜3年の月別データ（サンプルサイズ $T = 24 \sim 36$ 件程度）だと、時系列解析の文脈ではデータ数が少なすぎます。
この規模のデータ量では、共和分検定（ヨハンセン検定）の検出力が著しく低下し、実際には共和分関係があってもランク=0（関係なし）と判定されてしまう恐れがあります。
可能であれば週次データに粒度を細かくしてサンプルサイズを確保するか、検定結果を過信せず、ビジネスドメインの知識を加味して慎重に判断する必要があります。
5. インパルス応答関数
5.1 インパルス応答関数
インパルス応答関数は、「ある変数にショック（インパルス）が加わったとき、他の変数が複数期間にわたってどのように反応（応答）するか」を定量的に分析する手法です。施策の即効性だけでなく、効果が何ヶ月間持続するのか（残存効果）、どの時点でピークを迎えるのかを直感的に把握できます。
再度、$\text{VAR}(p)$モデルを$\text{VAR}(1)$モデルで表したcompanion formを考えます。
$$\tilde Y_t = \tilde \mu + A \tilde Y_{t-1} + \dots + \tilde \epsilon_t$$
とくに、以下のような状況を考えます。
- $t<0$なるすべての$t\in \mathbb{Z}$について、$Y_t=0.$
- $t>0$なるすべての$t\in \mathbb{Z}$について、$\epsilon_t=0.$
- $Y_0 = \epsilon_t = e_l$
- $t=0$における$l$番目の要素の1単位当たりのショック
ただし、$Y_t, \epsilon_t, Y_0, \epsilon_t, e_l$はそれぞれ$m \times 1$のベクトルで、
$e_l = \begin{pmatrix}
0 \\
\vdots \\
0 \\
1 \\
0 \\
\vdots \\
0
\end{pmatrix}.$
よって、companion formで表すと、
$$\tilde Y_0 = \begin{pmatrix}
Y_0 \\
0 \\
\vdots \\
0
\end{pmatrix} =
\tilde \epsilon_0 =
\begin{pmatrix}
\epsilon_0 \\
0 \\
\vdots \\
0
\end{pmatrix} =
\begin{pmatrix}
e_l \\
0 \\
\vdots \\
0
\end{pmatrix}.$$
同様にして、
$\tilde Y_1 = A \tilde Y_0 \quad(t=1)$
$\tilde Y_2 = A \tilde Y_1 = A^2 \tilde Y_0 \quad(t=2)$
$\vdots$
よって、
$\tilde Y_j = A^j Y_0 \quad(t=j)$
ただし$j$は正の整数。
$\tilde Y_0 =
\begin{pmatrix}
e_l \\
0 \\
\vdots \\
0
\end{pmatrix}$であることを考慮すると、$A^j \tilde Y_0$は$A^j$の$l$番目の列であることがわかります。
いま、$m\times mp$ のベクトル
$$J = (I_m \ 0\cdots \ 0)$$
を考えます。
さらに、$j = 0,1,2\cdots$について$\Psi_j = JA^JJ'$が$\text{VMA}(\infty)$表現：
$$Y_t = \Sigma_{j=0}^{\infty}\Psi_j\epsilon_{t-j} $$
の係数行列だとします。
5.2 直交化
ここですべての$k \neq l$とすべての$t$について$\epsilon_{kt}$が$\epsilon_{lt}$と無相関であるとするのは現実的でなく、実際は相関しているものと考えられます。
$\Sigma_{\epsilon}$は正定値行列だと考えられるので、コレスキー分解を行います：
$$\Sigma_{\epsilon} = PP'.$$
ただし、$P$は正の対角要素をもつ下三角行列です：
$$
P =
\begin{pmatrix}
p_{11} & 0 & \cdots & 0 \\
p_{21} & p_{22} & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
p_{m1} & p_{m2} & \cdots & p_{mm} \\
\end{pmatrix} \quad p_{kk} >0 \ \text{for} \ k\in\lbrace1,\cdots, m \rbrace.
$$
この分解により、$\text{VMA}(\infty)$表現は以下のように書き換えられます：
$Y_t = \Sigma_{j=0}^{\infty}\Psi_j PP^{-1}\epsilon_{t-j} = \Sigma_{j=0}^\infty \Theta_ju_{t-j}.$
ただし、$\Theta_j = \Psi_j P$で、$u_{t-j} = P^{-1}\epsilon_{t-j}.$
この変形のもとで、
$$E(u_t u_t') = E(P^{-1}\epsilon_t\epsilon_t'P'^{-1}) = P^{-1}E(\epsilon_t\epsilon_t')P'^{-1} = P^{-1}\Sigma_{\epsilon}P'^{-1} = P^{-1}PP'P'^{-1} = I_m.$$
よって、すべての$k \neq l$とすべての$t$について$u_{kt}$は$u_{lt}$と無相関となります。
$u_t$と$\epsilon_t$は以下のように関係づけられます：
$$\epsilon_t = Pu_t.$$
あるいは、
$$
\begin{pmatrix}
\epsilon_{1t} \\
\epsilon_{2t} \\
\vdots \\
\epsilon_{mt}
\end{pmatrix} =
\begin{pmatrix}
p_{11} & 0 & \cdots & 0 \\
p_{21} & p_{22} & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
p_{m1} & p_{m2} & \cdots & p_{mm} \\
\end{pmatrix}
\begin{pmatrix}
u_{1t} \\
u_{2t} \\
\vdots \\
u_{mt}
\end{pmatrix}.
$$
したがって、以下が成立します：
$\epsilon_{1t} = p_{11}u_{1t}$
$\epsilon_{2t} = p_{21}u_{1t} + p_{22}u_{2t} = \frac{p_{21}}{p_{11}}\epsilon_{1t} + p_{22}u_{2t}$
$\quad \quad \vdots$
$\epsilon_{mt} = p_{m1}u_{1t} + \cdots + p_{mm}u_{mt} = \frac{p_{m1}}{p_{11}}\epsilon_{1t} + \cdots + p_{mm}u_{mt}.$
注意すべきは、例えば$\epsilon_{1t}$の変動は$\epsilon_{2t}$の変動に影響しうる一方、逆は成立しないことです。コレスキー分解は上に位置する変数ほどドミナントで、下にある変数へ即座に影響を与え、逆は成り立たないという構造を強制します。
インパルス応答の結果は変数の並び順によって変化するため、「指名検索 $\rightarrow$ サイト訪問 $\rightarrow$ 応募」といったマーケティングのファネル構造（ドメイン知識）に基づいて、妥当な変数の順序を明示的に定義する必要があります。
6. 実装イメージ
# ─── 1. ファネル構造に沿った擬似データの生成 ───
set.seed(42)
n <- 48 # 4年分（48ヶ月）の月次データを想定
data_raw_df <- tibble::tibble(month = seq(from = as.Date("2022-01-01"),
to = as.Date("2025-12-01"),
by = "month"),
search = cumsum(rnorm(n, mean = 2.0, sd = 5.0)) + 1000) |> # 指名検索数
dplyr::mutate(visit = 0.7 * search + cumsum(rnorm(n, mean = 1.0, sd = 3.0)) + 500, # サイト訪問数
cv = 0.3 * visit + 0.1 * search + cumsum(rnorm(n, mean = 0.5, sd = 1.5)) + 100) # 応募数
# ─── 2. 単位根検定（ADF）と差分化（定常化） ───
tseries::adf.test(data_raw_df$search)
tseries::adf.test(data_raw_df$visit)
tseries::adf.test(data_raw_df$cv)
# 指名検索数とサイト訪問数が非定常系列であるため、前月比（差分系列）を作成して定常化
data_diff_df <- data_raw_df |>
dplyr::mutate(dplyr::across(-c(month), ~ . - dplyr::lag(.))) |>
dplyr::slice(-c(1))
# ─── 3. VARモデルのラグ選定と推定 ───
# コレスキー分解の再帰的構造（ファネルの順序）を明示的に指定
# 順序: [1] 指名検索(最も外生) -> [2] サイト訪問 -> [3] 応募 (最も内生)
data_var_input <- data_diff_df |>
dplyr::select(search, visit, cv)
# AIC基準で最適なラグ長を選択
lag_select <- vars::VARselect(data_var_input, lag.max = 3, type = "const")
optimal_lag <- lag_select$selection["AIC(n)"]
# VARモデルの推定
var_estimated <- vars::VAR(data_var_input, p = optimal_lag, type = "const")
summary <- summary(var_estimated)
# ─── 4. インパルス応答関数（IRF）の計算 ───
# 「指名検索」の1単位ショックが最終的な「応募数」に与える影響を計算
imp <- vars::irf(var_estimated, n.ahead = 10, boot = TRUE, ci = 0.95, runs = 100)
data.frame(month = 0:10,
impulse = as.vector(imp$irf$search),
upper = as.vector(imp$Upper$search),
lower = as.vector(imp$Lower$search),
to = rep(c("指名検索数", "サイト訪問数", "応募数"), c(11,11,11))) |>
dplyr::filter(to == "応募数") |>
ggplot2::ggplot(ggplot2::aes(x = month, y = impulse)) +
ggplot2::geom_ribbon(ggplot2::aes(ymax = upper, ymin = lower), alpha = 0.2, linetype = 2) +
#ggplot2::facet_wrap(.~ to, nrow = 5) +
ggplot2::geom_line(ggplot2::aes(color = to)) +
ggplot2::scale_x_continuous(breaks = 0:10) +
ggplot2::labs(x = "何カ月先の影響か", y = "応募数への影響") +
#ggplot2::ylim(-1000, 1000) +
ggplot2::geom_hline(yintercept = 0, linetype = 2) +
ggplot2::theme_minimal()
なお、VARモデルを推定してインパルス応答関数を出す際は、誤差項1標準偏差分のショックが加えられた時の、各期の目的変数の変動を見ています。例えば指名検索数に与えられた誤差項1標準偏差のサイズは以下で確認できます。
summary$varresult$search$sigma
[1] 5.001317
つまり、当月に約5件の指名検索が増えれば、当月の応募数が約1.4件増えるということになります（現実的ではないですが）。