这是用python做计量经济学实证分析系列第二篇。
原文地址：https://python.quantecon.org/mle.html
[scode type="green"]本系列文章对原文进行了一定的修改，主要是作图工具从matplotlib变为plotly.express，因为后者更为简洁明了，如果你不习惯plotly，可以看原文的代码[/scode]
[scode type="yellow"]本系列文章需要你对计量经济学有一定的了解，并且熟悉python、pandas、numpy[/scode]
引入
在上一节课中，我们使用线性回归来估计因变量和解释变量之间的关系。但如果我们的模型并不是线性关系呢？一种广泛使用的替代方法是最大似然估计，它包括指定一类由未知参数索引的分布，然后使用数据来确定这些参数值。相对于线性回归的好处是，它使变量之间的概率关系具有更大的灵活性。
这里我们通过复制丹尼尔·特雷斯曼（Daniel Treisman，2016）的论文《俄罗斯的亿万富翁》（Russia's Billionals）来说明最大似然估计，该论文将一个国家的亿万富翁数量与其经济特征联系起来。本文的结论是，俄罗斯的亿万富翁数量高于市场规模和税率等经济因素的预测。
文献地址：http://pubs.aeaweb.org/doi/pdfplus/10.1257/aer.p20161068
导入库：
import numpy as np
from numpy import exp
from scipy.special import factorial
import pandas as pd
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import statsmodels.api as sm
from statsmodels.api import Poisson
from scipy import stats
from scipy.stats import norm
from statsmodels.iolib.summary2 import summary_col
import plotly.express as px
import plotly
plotly.offline.init_notebook_mode(connected=True)开始
简介
最大似然估计的第一步是选择数据的概率分布。更准确地说，我们需要假设参数服从什么分布类。如正态分布或Gamma分布。
每一类都是由有限个参数索引的分布族。如正态分布是一类均值为$\mu \in(-\infty, \infty)$，标准差$\sigma \in(0, \infty)$的正态分布。
我们将通过固定参数的分布让数据挑选出特定元素。由此产生的参数估计称为最大似然估计。
亿万富翁
Treisman对估算不同国家的亿万富翁数量很感兴趣。亿万富翁的数量是整数值。
因此我们考虑只取非负整数的分布。（这就是为什么最小二乘回归不是解决当前问题的最佳工具，因为线性回归中的因变量不限于整数值）
一个整数分布是Poisson分布，其概率质量函数（probability mass function，pmf）为：
$${f(y)}=\frac{\mu^{y}}{y !} e^{-\mu}, \quad y=0,1,2, \ldots, \infty$$
对于u的不同值，我们可以绘制y上的Poisson分布，如下所示
poisson_pmf = lambda y, μ: μ**y / factorial(y) * exp(-μ)
df_fin = pd.DataFrame()
for μ in [1, 5, 10]:
df = pd.DataFrame()
distribution = [poisson_pmf(y_i, μ) for y_i in range(0, 25)]
df['y'] = list(range(0, 25))
df['f'] = distribution
df['μ'] = μ
df_fin = df_fin.append(df)
fig = px.line(df_fin, x="y", y="f",color="μ",labels = {'f':'f(y/μ)'})
fig.update_traces(mode='markers+lines')
plotly.offline.plot(fig)
注意，随着y的平均值的增加，Poisson分布开始类似于正态分布。
让我们来看看我们在这堂课中要处理的数据的分布情况。Treisman的主要数据来源是《福布斯》对亿万富翁及其估计净值的排名。原数据集在这：https://www.aeaweb.org/articles?id=10.1257/aer.p20161068
df = pd.read_stata('https://cdn.jsdelivr.net/gh/QuantEcon/QuantEcon.lectures.code/mle/fp.dta')
print(df.head())
使用柱状图，我们可以看到2008年每个国家亿万富翁人数的分布，(numbil0)（出于作图目的舍弃了美国的数据）
numbil0_2008 = df[(df['year'] == 2008) & (df['country'] != 'United States')]
fig = px.histogram(numbil0_2008, x="numbil0",labels={'numbil0':'Number of billionaires in 2008'})
plotly.offline.plot(fig)
从直方图来看，泊松假设并非不合理（尽管u值很低，且存在一些异常值）
条件分布
在Treisman的论文中，因变量——i国亿万富翁的数量$y_i$——它被模拟为人均国内生产总值、人口规模以及加入关贸总协定和世贸组织的年数的函数。因此，yi的分布需要以解释变量$x_i$的向量为条件。标准的poisson回归模型如下：
$$\begin{array}{c}f\left(y_{i} | \mathbf{x}_{i}\right)=\frac{\mu_{i}^{y_{i}}}{y_{i} !} e^{-\mu_{i}} ; \quad y_{i}=0,1,2, \ldots, \infty \qquad(1)\\ \text { where } \mu_{i}=\exp \left(\mathbf{x}_{i}^{\prime} \beta\right)=\exp \left(\beta_{0}+\beta_{1} x_{i 1}+\ldots+\beta_{k} x_{i k}\right)\end{array}$$
为了说明yi的分布依赖于xi，让我们做一个简单的拟合。使用上面的poisson_pmf函数，并代入$\beta$和$x_i$的任意值。
poisson_pmf = lambda y, μ: μ**y / factorial(y) * exp(-μ)
df_fin = pd.DataFrame()
# Define a parameter vector with estimates
β = np.array([0.26, 0.18, 0.25, -0.1, -0.22])
# Create some observations X
datasets = [np.array([0, 1, 1, 1, 2]),
np.array([2, 3, 2, 4, 0]),
np.array([3, 4, 5, 3, 2]),
np.array([6, 5, 4, 4, 7])]
for X in datasets:
μ = round(exp(X @ β),3)
distribution = [poisson_pmf(y_i, μ) for y_i in range(0, 25)]
df = pd.DataFrame()
df['y'] = list(range(0, 25))
df['f'] = distribution
df['μ'] = μ
df_fin = df_fin.append(df)
fig = px.line(df_fin, x="y", y="f",color="μ",
labels = {'f':r'$f(y \mid x_i; \beta)$',,'y':'$y \mid x_i$'})
fig.update_traces(mode='markers+lines')
plotly.offline.plot(fig,include_mathjax='cdn')
我们可以看到yi的分布是以xi为条件的；（ui不再是常数）
最大似然估计
在我们的亿万富翁模型中，条件分布包含4个我们需要估计的参数（k=4）。我们将整个参数向量标记为$\beta$，其中
$$\beta=\left[\begin{array}{l}\beta_{0} \\ \beta_{1} \\ \beta_{2} \\ \beta_{3}\end{array}\right]$$
为了用极大似然估计法估计模型，我们希望使我们的估计$\hat\beta$是真参数$\beta$的可能性最大。首先，我们需要构造似然函数（$\mathcal{L}(\beta)$），它类似于联合概率密度函数(probability density function，pdf)。
假设我们有一些数据$y_{i}=\left\{y_{1}, y_{2}\right\} \text { and } y_{i} \sim f\left(y_{i}\right)$。如果y1和y2是独立的，则这些数据的联合pmf是$f\left(y_{1}, y_{2}\right)=f\left(y_{1}\right) \cdot f\left(y_{2}\right)$。如果$y_i$服从Poisson分布，$\lambda=7$，我们可以可视化联合pmf如下：
import matplotlib.pyplot as plt
def plot_joint_poisson(μ=7, y_n=20):
poisson_pmf = lambda y, μ: μ**y / factorial(y) * exp(-μ)
yi_values = np.arange(0, y_n, 1)
# Create coordinate points of X and Y
X, Y = np.meshgrid(yi_values, yi_values)
# Multiply distributions together
Z = poisson_pmf(X, μ) * poisson_pmf(Y, μ)
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z.T, cmap='terrain', alpha=0.6)
ax.scatter(X, Y, Z.T, color='black', alpha=0.5, linewidths=1)
ax.set(xlabel='$y_1$', ylabel='$y_2$')
ax.set_zlabel('$f(y_1, y_2)$', labelpad=10)
plt.show()
plot_joint_poisson(μ=7, y_n=20)不会用plotly做类似的图，只能照抄源码用matplotlib.pyplot了。
类似地，我们数据的联合pmf（作为条件Poisson分布）可以写成
$$f\left(y_{1}, y_{2}, \ldots, y_{n} | \mathbf{x}_{1}, \mathbf{x}_{2}, \ldots, \mathbf{x}_{n} ; \beta\right)=\prod_{i=1}^{n} \frac{\mu_{i}^{y_{i}}}{y_{i} !} e^{-\mu_{i}}$$
yi取决于$x_i$和参数$\beta$的值。似然函数与联合PMF相同，但将参数$\beta$作为随机变量对待，并给出观测值（$y_i$，$x_i$）。
$$\begin{aligned} \mathcal{L}\left(\beta | y_{1}, y_{2}, \ldots, y_{n} ; \mathbf{x}_{1}, \mathbf{x}_{2}, \ldots, \mathbf{x}_{n}\right) &=\prod_{i=1}^{n} \frac{\mu_{i}^{y_{i}}}{y_{i} !} e^{-\mu_{i}} \\ &=f\left(y_{1}, y_{2}, \ldots, y_{n} | \mathbf{x}_{1}, \mathbf{x}_{2}, \ldots, \mathbf{x}_{n} ; \beta\right) \end{aligned}$$
既然我们有了似然函数，我们想找到产生最大似然值的$\beta$
$$\max _{\beta} \mathcal{L}(\beta)$$
log-likelihood通常更容易最大化（考虑区分${f(x)}=x \exp (x) \text { vs. } {f(x)}=\log (x)+x$）
假设取对数是单调递增变换，则似然函数的极大值也是对数似然函数的极大值。在我们的例子中，log-likelihood是：
$$\begin{aligned} \log \mathcal{L}(\beta)=& \log \left(f\left(y_{1} ; \beta\right) \cdot f\left(y_{2} ; \beta\right) \cdot \ldots \cdot f\left(y_{n} ; \beta\right)\right) \\ &=\sum_{i=1}^{n} \log f\left(y_{i} ; \beta\right) \\ &=\sum_{i=1}^{n} \log \left(\frac{\mu_{i}^{y_{i}}}{y_{i} !} e^{-\mu_{i}}\right) \\ &=\sum_{i=1}^{n} y_{i} \log \mu_{i}-\sum_{i=1}^{n} \mu_{i}-\sum_{i=1}^{n} \log y ! \end{aligned}$$
然而，上述问题不存在解析解-要找到最大似然估计，我们需要使用数值方法。
数值方法
许多分布没有很好的解析解，因此需要数值方法来求解参数估计。其中一种数值方法是牛顿-拉斐逊算法。我们的目标是找到最大似然估计$\hat\beta$。在$\hat\beta$处，log-likelihoo函数的一阶导数等于0。我们假设：
$$\log \mathcal{L}(\beta)=-(\beta-10)^{2}-10$$
β = np.linspace(1, 20)
logL = -(β - 10) ** 2 - 10
dlogL = -2 * β + 20
df = pd.DataFrame({'logL':logL,'dlogL':dlogL,'β':β})
df_long=pd.melt(df, id_vars=['β'], value_vars=['logL', 'dlogL'])
fig = px.line(df_long, x="β", y="value", facet_row="variable",line_shape="spline")
fig['layout']['yaxis']['title']['text']=r'$\frac{dlog \mathcal{L(\beta)}}{d \beta}$'
fig['layout']['yaxis2']['title']['text']=r'$log \mathcal{L(\beta)}$'
fig['layout']['annotations'][0]['text'] = ''
fig['layout']['annotations'][1]['text'] = ''
fig.update_yaxes(matches=None)
plotly.offline.plot(fig,include_mathjax='cdn')
上图显示了最大似然值（顶部图表）出现在$\frac{d \log \mathcal{L}(\beta)}{d \beta}=0$时。因此，当$\beta=10$时，似然值最大。
我们还可以通过检查二阶导数（底部曲线图的斜率）是否为负值来确保该值是最大值（而不是最小值）。
Newton-Raphson算法找到一个一阶导数为0的点。为了使用该算法，我们在最大值处进行初始猜测（OLS参数估计可能是合理的猜测），然后
1.使用更新规则迭代算法
$$\beta_{(k+1)}=\beta_{(k)}-H^{-1}\left(\beta_{(k)}\right) G\left(\beta_{(k)}\right)$$
其中：
$$\begin{aligned} G\left(\beta_{(k)}\right) &=\frac{d \log \mathcal{L}\left(\beta_{(k)}\right)}{d \beta_{(k)}} \\ H\left(\beta_{(k)}\right) &=\frac{d^{2} \log \mathcal{L}\left(\beta_{(k)}\right)}{d \beta_{(k)} d \beta_{(k)}^{\prime}} \end{aligned}$$
2.检查是否$\beta_{(k+1)}-\beta_{(k)}<t o l$
- 如果为true，则停止迭代并设置$\hat{\beta}=\beta_{(k+1)}$
- 如果为false，则更新$\beta_{(k+1)}$
从更新方程可以看出，只有当$G\left(\beta_{(k)}\right)=0$即一阶导数等于0时，$\beta_{(k+1)}=\beta_{(k)}$，（实际上，当差异低于一个可容忍的阈值时，我们就停止迭代）
让我们开始实现Newton-Raphson算法。首先，我们将创建一个名为poissonrege的类，这样我们就可以很容易地重新计算每次迭代的对数似然、梯度和Hessian值
class PoissonRegression:
def __init__(self, y, X, β):
self.X = X
self.n, self.k = X.shape
# Reshape y as a n_by_1 column vector
self.y = y.reshape(self.n,1)
# Reshape β as a k_by_1 column vector
self.β = β.reshape(self.k,1)
def μ(self):
return np.exp(self.X @ self.β)
def logL(self):
y = self.y
μ = self.μ()
return np.sum(y * np.log(μ) - μ - np.log(factorial(y)))
def G(self):
X = self.X
y = self.y
μ = self.μ()
return X.T @ (y - μ)
def H(self):
X = self.X
μ = self.μ()
return -(X.T @ (μ * X))我们的函数newton_raphson将取一个PoissonRegression 对象，该对象具有对参数向量$\beta_0$的初始猜测。该算法根据更新规则更新参数向量，并在新的参数估计下重新计算梯度矩阵和Hessian矩阵。
迭代将在下列任一情况下结束：
- 参数和更新参数之间的差异低于可容忍水平。
- 已达到最大迭代次数（意味着未达到收敛）
因此，我们可以了解算法运行时的情况，在每次迭代时添加一个选项display=True来打印值。
def newton_raphson(model, tol=1e-3, max_iter=1000, display=True):
i = 0
error = 100 # Initial error value
# Print header of output
if display:
header = f'{"Iteration_k":<13}{"Log-likelihood":<16}{"θ":<60}'
print(header)
print("-" * len(header))
# While loop runs while any value in error is greater
# than the tolerance until max iterations are reached
while np.any(error > tol) and i < max_iter:
H, G = model.H(), model.G()
β_new = model.β - (np.linalg.inv(H) @ G)
error = β_new - model.β
model.β = β_new
# Print iterations
if display:
β_list = [f'{t:.3}' for t in list(model.β.flatten())]
update = f'{i:<13}{model.logL():<16.8}{β_list}'
print(update)
i += 1
print(f'Number of iterations: {i}')
print(f'β_hat = {model.β.flatten()}')
# Return a flat array for β (instead of a k_by_1 column vector)
return model.β.flatten()让我们用一个X包含5个观测值和3个变量的小数据集来验证我们的算法。
X = np.array([[1, 2, 5],
[1, 1, 3],
[1, 4, 2],
[1, 5, 2],
[1, 3, 1]])
y = np.array([1, 0, 1, 1, 0])
# Take a guess at initial βs
init_β = np.array([0.1, 0.1, 0.1])
# Create an object with Poisson model values
poi = PoissonRegression(y, X, β=init_β)
# Use newton_raphson to find the MLE
β_hat = newton_raphson(poi, display=True)
由于这是一个观测较少的简单模型，6次迭代后就收敛了。可以看到，随着每次迭代，log似然值都会增加。记住，我们的目标是最大化对数似然函数，这是算法一直致力于实现的。
另外，注意随着迭代次数的增加，$\log \mathcal{L}\left(\beta_{(k)}\right)$的增加值变小了。这是因为当我们达到最大值时，梯度接近0，因此我们的更新方程中的分子变得更小。梯度向量在$\hat\beta$处应该接近于0
poi.G()
>>>
[[-3.952e-07]
[-1.001e-06]
[-7.731e-07]]迭代过程可以在下图中可视化，其中最大值出现在$\beta=10$处
logL = lambda x: -(x - 10) ** 2 - 10
def find_tangent(β, a=0.01):
y1 = logL(β)
y2 = logL(β+a)
x = np.array([[β, 1], [β+a, 1]])
m, c = np.linalg.lstsq(x, np.array([y1, y2]), rcond=None)[0]
return m, c
β = np.linspace(2, 18)
fig, ax = plt.subplots(figsize=(12, 8))
ax.plot(β, logL(β), lw=2, c='black')
for β in [7, 8.5, 9.5, 10]:
β_line = np.linspace(β-2, β+2)
m, c = find_tangent(β)
y = m * β_line + c
ax.plot(β_line, y, '-', c='purple', alpha=0.8)
ax.text(β+2.05, y[-1], f'$G({β}) = {abs(m):.0f}$', fontsize=12)
ax.vlines(β, -24, logL(β), linestyles='--', alpha=0.5)
ax.hlines(logL(β), 6, β, linestyles='--', alpha=0.5)
ax.set(ylim=(-24, -4), xlim=(6, 13))
ax.set_xlabel(r'$\beta$', fontsize=15)
ax.set_ylabel(r'$log \mathcal{L(\beta)}$',rotation=0,labelpad=25,fontsize=15)
ax.grid(alpha=0.3)
plt.show()
请注意，我们对Newton-Raphson算法的实现是相当基本的—有关更健壮的实现，请参见scipy.optimize。
使用statsmodels进行最大似然估计
知道了背后的原理，我们可以将MLE应用到一个有趣的应用程序中。我们将在statsmodels中使用Poisson回归模型来获得更丰富的输出，包括标准误差、测试值等等。
statsmodels使用与上面相同的算法来寻找最大似然估计。在我们开始之前，让我们用stats模型重新估计我们的简单模型，以确认我们得到了相同的系数和对数似然值。
X = np.array([[1, 2, 5],
[1, 1, 3],
[1, 4, 2],
[1, 5, 2],
[1, 3, 1]])
y = np.array([1, 0, 1, 1, 0])
stats_poisson = Poisson(y, X).fit()
print(stats_poisson.summary())
现在，让我们复制丹尼尔·特雷斯曼（Daniel Treisman）在报告中提到的俄罗斯亿万富翁论文的结果。
Treisman从估计等式（1）开始，式中：
- yi为亿万富翁数
- $x_{i1}$是人均GDP的对数
- $x_{i2}$是人口对数
- $x_{i3}$是加入GATT和世贸组织时长（代表进入国际市场）
本文只考虑2008年的估计。我们将像这样设置我们的估计变量
df = pd.read_stata('https://cdn.jsdelivr.net/gh/QuantEcon/QuantEcon.lectures.code/mle/fp.dta')
df = df[df['year'] == 2008]
df['const'] = 1
# Variable sets
reg1 = ['const', 'lngdppc', 'lnpop', 'gattwto08']
reg2 = ['const', 'lngdppc', 'lnpop','gattwto08', 'lnmcap08', 'rintr', 'topint08']
reg3 = ['const', 'lngdppc', 'lnpop', 'gattwto08', 'lnmcap08','rintr', 'topint08', 'nrrents', 'roflaw']然后利用statsmodels中的Poisson函数对模型进行拟合。在作者的论文中，我们将使用健壮的标准错误
# Specify model
poisson_reg = sm.Poisson(df[['numbil0']], df[reg1],missing='drop').fit(cov_type='HC0')
print(poisson_reg.summary())
成功！算法在9次迭代中收敛。
我们的结果表明，人均国内生产总值（GDP）和加入关税及贸易总协定（GATT）的年限与一个国家的亿万富翁数量正相关，正如预期的那样。
让我们来估计一下作者的更全功能的模型，并将它们显示在一个表中
regs = [reg1, reg2, reg3]
reg_names = ['Model 1', 'Model 2', 'Model 3']
info_dict = {'Pseudo R-squared': lambda x: f"{x.prsquared:.2f}",
'No. observations': lambda x: f"{int(x.nobs):d}"}
regressor_order = ['const','lngdppc','lnpop','gattwto08','lnmcap08','rintr',
'topint08','nrrents','roflaw']
results = []
for reg in regs:
result = sm.Poisson(df[['numbil0']], df[reg],
missing='drop').fit(cov_type='HC0',maxiter=100, disp=0)
results.append(result)
results_table = summary_col(results=results,float_format='%0.3f',stars=True,
model_names=reg_names,info_dict=info_dict,
regressor_order=regressor_order)
results_table.add_title('Table 1 - Explaining the Number of Billionaires in 2008')
print(results_table)
研究结果表明，亿万富翁的数量与人均GDP、人口规模、股市市值呈正相关，与最高边际所得税率呈负相关。为了按国家分析我们的结果，我们可以绘制出预测值和实际值之间的差异，然后从最高值到最低值排序，并绘制出前15个国家
data = ['const', 'lngdppc', 'lnpop', 'gattwto08', 'lnmcap08','rintr',
'topint08', 'nrrents', 'roflaw','numbil0','country']
results_df = df[data].dropna()
# Use last model (model 3)
results_df['prediction'] = results[-1].predict()
# Calculate difference
results_df['difference'] = results_df['numbil0'] - results_df['prediction']
# Sort in descending order
results_df = results_df.sort_values('difference', ascending=False)[:15]
# Plot the first 15 data points
fig = px.bar(results_df, x="country", y="difference")
plotly.offline.plot(fig,include_mathjax='cdn')
正如我们所看到的，俄罗斯的亿万富翁人数远远超过了该模型的预测（比预期多出约50人）Treisman利用这一实证结果讨论了俄罗斯亿万富翁过剩的可能原因，包括俄罗斯财富的来源、政治气候以及苏联解体后的私有化历史。
总结
在本课程中，我们使用最大似然估计来估计泊松模型的参数。statsmodels包含其他内置的似然模型，如Probit和Logit。可以在这里找到示例：https://www.statsmodels.org/dev/examples/notebooks/generated/generic_mle.html
版权属于：作者名称
本文链接：https://www.sitstars.com/archives/89/
转载时须注明出处及本声明
太高森了
Ryan 2020-06-08