这是用python做计量经济学实证分析系列第二篇。

原文地址:https://python.quantecon.org/mle.html

本系列文章对原文进行了一定的修改,主要是作图工具从matplotlib变为plotly.express,因为后者更为简洁明了,如果你不习惯plotly,可以看原文的代码

本系列文章需要你对计量经济学有一定的了解,并且熟悉python、pandas、numpy

引入

在上一节课中,我们使用线性回归来估计因变量和解释变量之间的关系。但如果我们的模型并不是线性关系呢?一种广泛使用的替代方法是最大似然估计,它包括指定一类由未知参数索引的分布,然后使用数据来确定这些参数值。相对于线性回归的好处是,它使变量之间的概率关系具有更大的灵活性。

这里我们通过复制丹尼尔·特雷斯曼(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

Last modification:May 12th, 2020 at 11:46 am