操作変数法との付き合い方を考えてみた
#はじめに
因果推論の方法で操作変数法という手法がある。
この手法は、未観測の欠落変数バイアスにによる処置の効果のバイアス取り除くための手法だ。
ただし、前提とする仮定が厳しく、実務において単純に実施するのは難しい。
#難しいポイント
操作変数法は前提とする仮定が厳しく特に下記の点が実務への応用を難しくさせている
・結果変数とは処置を通してのみ影響を与える変数が良い
・他の変数とは独立な変数が必要
・全ての候補を洗い出すのは事実上不可能
この中でも、特に1つ目の結果変数とは処理を通してのみ影響を与える変数が良いと言うのが難しい。
#価格変更による弾力性で考えると
例えば、EC商品における、月次の需要に対する価格弾力性の因果効果を推定すると仮定し、操作変数として利用可能な変数を考えると下記のようなものが考えられる。
操作変数候補
・社内による価格改定ルール
・機械的な価格改定ルール自体は価格に直接影響を与えるが需要へは価格の変化を通じてのみ影響をあたえるため
・輸送コストの増加
・輸送費は材料の調達などのコストであり、価格に反映されるが、購入点数に直接影響を与えることが少ないためコントロール変数候補
結果変数である需要へ直接影響を与えるが、価格弾力性の推定に利用可能な変数
・原価、在庫、競合他社価格、同ブランドの価格帯、同カテゴリの価格帯、キャンペーン、季節性、地域別供給コスト、天気気温
しかし、機械的な価格変更ルールも顧客が知っていた場合、価格変更前の駆け込み需要の発生など結果変数である需要に直接影響を与えるケースや、輸送費が増加することにより、在庫が減少し需要が下がる可能性が考えられるなど、実際に結果変数に直接影響を与えない完璧な操作変数を選択するのは難しい。
#コラム:価格変更による弾力性の2SLSモデル
$$
\log y_i = \alpha + \pi \log P_i + X_i^\top \beta + \varepsilon_i
$$
$$
\log P_i = a + \gamma Z_i + X_i^\top \delta + \nu_i
$$
$${y_i}$$:購入点数
$${P_i}$$:価格
$${X_i}$$:コントロール変数
$${\pi}$$:価格弾力性
$${\varepsilon_i}$$:未観測需要ショック
$${Z_i}$$:操作変数
$${\gamma}$$:操作変数が価格に与える影響
$${\delta}$$:first stageにおけるコントロール変数の係数
$${\nu_i}$$:価格のうち、$${Z_i}$$ と $${X_i}$$ で説明できない部分
#対処法
完璧な操作変数を探すのではなく、仮定を緩め現実的に用意が可能な変数を用いて評価すること
・仮定を満たしていない場合の影響を因果効果に含めた評価
・ドメイン知識を元に仮定の妥当性を議論
#手法.Conley
操作変数の結果変数への直接影響を許容した上で、影響の程度を組み込むことで、妥当な因果効果の範囲を推定する手法
$$
Y = \beta D + \gamma Z + u
$$
$${\gamma}$$:操作変数の直接影響
$$
\beta = \frac{Cov(Y,Z) - \gamma Var(Z)}{Cov(D,Z)}
$$$${\gamma = 0}$$:操作変数が結果変数への直接影響がない状態
$${\gamma \neq 0}$$:操作変数が結果変数への直接影響がある状態
#サンプルデータで実装
CigarettesSWデータの概要
・件数:96件
・カラム;州、年、消費者物価指数、1人あたりのタバコ消費量、州税、所得、タバコ価格等
・概要:アメリカの州別・年次のパネルデータであり、タバコの価格と需要の関係に加え、税という価格に影響を与える外生的な変動要因が含まれている点が特徴であり、操作変数法の典型的な例として利用できる
○実施内容
・タバコの価格弾力性を州税を操作変数として推定
・Conleyの手法を用いて、操作変数の結果変数への直接影響を許容した妥当な因果効果の範囲を推定
○結果

排除制約をわずかに緩和した場合(γ ≈ -0.02)で推定値がゼロを横切り、 γ ≈ -0.05で符号が反転することが確認される。つまり、小さな条件緩和でも結果が変わる脆弱性があることが伺える。この場合、今回使用した税率が操作変数に適していたかを考え直す必要がある。
#まとめ
IV法は「使えない手法」ではない。ただし、変数を入れれば因果効果が出る手法でもない。実務では完全な識別は困難なため、仮定を満たしていない場合の影響を因果効果の評価やドメイン知識を元に仮定の妥当性を担保した上で、「意思決定」する必要がある
#参考
・Plausibly Exogenous
・因果効果入門-ミックステープ-
#サンプルコード
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
from linearmodels.iv import IV2SLS
#. Cigarette Demand datasetを利用
df = sm.datasets.get_rdataset("CigarettesSW", "AER").data
#. 物価指数を価格/所得に反映
df["real_price"] = df["price"] / df["cpi"]
df["real_income"] = df["income"] / df["cpi"]
#. log変換
df["log_packs"] = np.log(df["packs"])
df["log_price"] = np.log(df["real_price"])
df["log_income"] = np.log(df["real_income"])
# ダミー変数作成
df_tf = pd.get_dummies(df, columns=["state", "year"], drop_first=True)
#. OLS推定
#. 説明変数
X = df_tf[["log_price", "log_income"] +
[c for c in df_tf.columns if c.startswith("state_") or c.startswith("year_")]]\
.astype(float)
X_sm = sm.add_constant(X)
ols = sm.OLS(df["log_packs"], X_sm).fit()
print(ols.summary())
#. IV
#. 説明変数
X = df_tf[["log_income"] +
[c for c in df_tf.columns if c.startswith("state_") or c.startswith("year_")]]\
.astype(float)
iv = IV2SLS(
dependent=df["log_packs"],
exog=X,
endog=df["log_price"],
instruments=df["tax"]
).fit()
print(iv.first_stage)
print(iv.summary)
#. 仮定を違反した場合の影響
def conley_bounds(df, gamma_grid):
X = df[["log_income"] +
[c for c in df_tf.columns if c.startswith("state_") or c.startswith("year_")]]\
.astype(float)
X_sm = sm.add_constant(X)
y_res = sm.OLS(df["log_packs"], X_sm).fit().resid
d_res = sm.OLS(df["log_price"], X_sm).fit().resid
z_res = sm.OLS(df["tax"], X_sm).fit().resid
cov_yz = np.cov(y_res, z_res, bias=True)[0,1]
cov_dz = np.cov(d_res, z_res, bias=True)[0,1]
var_z = np.var(z_res)
betas = []
for gamma in gamma_grid:
beta = (cov_yz - gamma * var_z) / cov_dz
betas.append(beta)
return np.array(betas)
gamma_grid = np.linspace(-0.5, 0.1, 200)
betas = conley_bounds(df_tf, gamma_grid)
result_df = pd.DataFrame({
"gamma": gamma_grid,
"beta": betas
})
#. βが0になるγ
zero_cross = result_df.iloc[(result_df["beta"].abs()).argsort()[:1]]
print(zero_cross)
#. 図の描画
plt.plot(gamma_grid, betas)
plt.axhline(0, linestyle="--")
plt.xlabel("gamma")
plt.ylabel("price elasticity")
plt.title("Conley sensitivity")
plt.show()