Table of Contents
対数変換と回帰モデル
\[\mathbb E[y|\mathbf x] = \exp(\beta_0 + \beta_1 \ln(x_1) + \beta_2x_2)\]という回帰モデルについて弾力性を考えた場合, $x_1$についての弾力性は以下で表現されます
\[\text{elasticity} = \frac{\partial\ln(\mathbb E[y|\mathbf x])}{\partial\ln (x_1)} = \beta_1\]実務で上記を推定する場合は, 線形関数に直したほうが扱いやすいという理由で以下の推定式を用いるケースがあります
\[\begin{align} &\ln(\mathbb E[y|\mathbf x]) = \beta_0 + \beta_1 \ln(x_1) + \beta_2x_2\\ &\Rightarrow \ln(y) = b_0 + b_1 \ln(x_1) + b_2x_2 + \epsilon \end{align}\]| (1), (2)は似た推定式ですが, 厳密に言うと(2)は$\mathbb E[\ln(y) | x]$をtargetとして推定しているので, 一般的には |
となってしまいます. このstatementは（厳密ではないですが）Jensen’s inequalityで確かめられます
Jensen’s inequality
今回の対数関数はconcaveなのでconcave functionに限ってJensen’s inequalityを紹介します.
Theorem: Jensen’s Inequality
If $g(x)$ is a concave function on $\mathbf R_X$, and $\mathbb E[g(X)]$ and $\mathbb g(\mathbb E[X])$ are finite, then
\[\mathbb E[g(X)] \leq g(\mathbb E[X])\]Sketch of Proof
関数$g$がconcaveならば, 定義域の任意の点$x_0$で以下の関係式が成立します
\[g(x) \leq g(x_0) + g'(x_0)(x - x_0), \forall x \in R_x\]上記の式について$x_0 = \mathbb E[x]$とした上で, 期待値をとると
\[\begin{align*} \mathbf E[g(x)] &\leq \mathbb E[g(\mathbb E[x]) + g'(\mathbb E[x])(x - \mathbb E[x])]\\ &= g(\mathbb E[x]) + g'(\mathbb E[x])(\mathbb E[x] - \mathbb E[x])\\ &= g(\mathbb E[x]) \end{align*}\]なお, strictly concave & $x$が「almost surelyの意味でconstant」ではないならば
\[\mathbb E[g(X)] < g(\mathbb E[X])\]が成立します.
対数変換でパラメーターが識別できるときはあるのか？
= \mathbb E[u] Jensen’s Inequalityによって, 対数変換のようなstrictly concaveな非線形関数のとき
\[\mathbb E[\ln(y)|\mathbf x] \neq \ln(\mathbb E[y|\mathbf x])\]となりますが, それでも弾力性を識別できるときはあります. Intuitionだけ説明したいので, 上記のTrue form回帰モデルを 以下のようにして説明します.
\[\ln(y) = \beta_0 + \beta_1 \ln(x_1) + \beta_2x_2 + u\]このとき, 両辺を指数変換 & $X$で条件づけた期待値をとると
\[\mathbb E[y|X] = \exp(\beta_0 + \beta_1 \ln(x_1) + \beta_2x_2)\times\mathbb E[\exp(u)|X]\]となります. このとき,
\[\mathbb E[\exp(u)|X] = \mathbb E[\exp(u)] = 1\]が成り立つならば, 簡単な計算で
\[\frac{\partial\ln(\mathbb E[y|\mathbf x])}{\partial\ln (x_1)} = \frac{\partial\mathbb E[\ln(y)|\mathbf x]}{\partial\ln (x_1)}\]が確かめられます.
識別条件
\[\mathbb E[\exp(u)|X] = \mathbb E[\exp(u)] = 1\]この関係式は以下のとき成立する
\[\begin{align*} &X\text{と}u\text{が独立} \\ &\mathbb E[u|X] = \mathbb E[u] = 0 \end{align*}\]conditional meanが0, i.e., $\mathbb E[u|X] = 0$の条件だけでは不十分です. 例として, $u$の分散が$x$の大きさに比例するなどのとき(Heteroskedasticity)は, 「$X$と$u$が独立」とはいえません.
実務上差し障りがあるのか？
実務上は対数変換したとしても, 推定係数が大きく変わるということはありません. Appendixのコードで確かめたところ 小さなズレはありましたが, あまり差し障りがないと判断しても良いと思います.
ただし, True formで推定したほうが推定値の分散が小さくなるのは間違いないことです.
\[y = \exp(\beta_0 + \beta_1 x_1) + \epsilon\]上記をtrue formとしたとき,
- error term$\epsilon$がregressorの逆数に比例する条件で不均一性をもたせ
- その上で分散の大きさを変化させてみた
という形でsimulation dataを作成して回帰したところ, $\beta_!$の推定値の分散はnoiseの大きさに対して以下のような結果になりました. ただし, orderがnanoの単位なので現実的にはあまり気にしなくてもよいと思われます.
Appendix: Python Code
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
import plotly.express as px
import numpy as np
import statsmodels.api as sm
from scipy import stats
from scipy.optimize import minimize
import pandas as pd
class LogExpReg:
def __init__(self, n, v, max_counter:int=100, seed:int=42):
self.ob_size = n
self.noise_level = v
self.max_counter = max_counter
self.exp_form_beta = []
self.log_form_beta = []
np.random.seed(seed)
def dgp(self):
CONST = 20
beta = 1
self.x = np.random.uniform(0, 1, self.ob_size)
error = np.array(list(map(lambda x: np.random.normal(0, x, 1),
self.noise_level/np.abs(self.x)))).flatten()
self.y = np.exp(CONST + beta*self.x) + error
def log_reg(self):
outcome = np.log(self.y)
regressor = np.array([np.array([1]*len(self.x)), self.x]).T
self.log_beta = np.linalg.inv(regressor.T @ regressor) @ (regressor.T @ outcome)
def gaussian_mle(self, params):
beta_0 = params[0]
beta_1 = params[1]
sd = params[2]
ce = np.exp(beta_0 + beta_1 * self.x)
obj = -stats.norm.logpdf(self.y, loc=ce, scale=sd).sum()
return obj
def exp_reg(self):
init_sigma = np.sqrt(np.mean((self.y - np.exp(self.log_beta[0] + self.log_beta[1] * self.x))**2))
initParams = [self.log_beta[0], self.log_beta[1], init_sigma]
results = minimize(self.gaussian_mle, initParams, method='Nelder-Mead')
self.exp_beta = results.x[:2]
def fit(self):
counter = 0
while counter < self.max_counter:
counter += 1
self.dgp()
self.log_reg()
self.exp_reg()
self.log_form_beta.append(self.log_beta[1])
self.exp_form_beta.append(self.exp_beta[1])
def plot(self):
df = pd.DataFrame({'log_beta':self.log_form_beta,
'exp_beta':self.exp_form_beta})
fig = px.histogram(df, x=['log_beta', 'exp_beta'],
labels={'value': 'beta_1'},
histnorm='probability',
barmode='overlay')
return fig
def summary(self):
exp_beta_mean, log_beta_mean = np.mean(self.exp_form_beta), np.mean(self.log_form_beta)
exp_beta_var, log_beta_var = np.var(self.exp_form_beta), np.var(self.log_form_beta)
return exp_beta_mean, log_beta_mean, exp_beta_var, log_beta_var
#-------------------
# main
#-------------------
from os import cpu_count
from multiprocessing import Pool
def main(var_size):
Solver = LogExpReg(50, var_size)
Solver.fit()
return Solver.summary()
num_workers = int(0.8 * cpu_count())
with Pool(num_workers) as p:
result = np.array(p.map(main, np.arange(10, 1000, 10)))
df = pd.DataFrame(result, columns=['exp_mean', 'log_mean', 'exp_var', 'log_var'])
df['noise_sigma'] = np.arange(10, 1000, 10)
fig = px.line(df, x = 'noise_sigma', y = ['exp_var', 'log_var'])
fig.show()
References
(注意：GitHub Accountが必要となります）