高级计量经济学 11：最大似然估计(上)

此文内容为《高级计量经济学及STATA应用》的笔记，陈强老师著，高等教育出版社出版。

我只将个人会用到的知识作了笔记，并对教材较难理解的部分做了进一步阐述。为了更易于理解，我还对教材上的一些部分（包括代码和正文）做了修改。

仅供学习参考，请勿转载，侵删！

6 最大似然估计法
- 6.5 最大似然法的大样本性质
  - 6.5.1 估计量的一致性
  - 6.5.2 渐近有效和渐近正态
- 6.6 最大似然估计量的渐近协方差矩阵
  - 6.6.1 期望值法
  - 6.6.2 观测信息矩阵法
  - 6.6.3 梯度向量外积或BHHH法

$\S \text{ 第 6 章 } \S$

$\text{最大似然估计}$

6 最大似然估计法

6.1 最大似然估计法的定义

假设随机变量 $y$ 的概率密度函数为 $f(y;\pmb \theta)$ ，其中 $\pmb \theta$ 为 $K$ 维位置参数向量， $\pmb \theta \in \Theta$ 表示参数向量处于参数空间中，即所有 $\pmb \theta$ 可能取值所构成的集合。我们现在通过抽取随机样本 $\{y_1,y_2, \cdots,y_n \}$ 来估计 $\pmb \theta$ 。假设 $\{y_1,y_2, \cdots,y_n \}$ 为 $\text{i.i.d.}$ ，那么样本数据的联合密度函数为： $f(y_1;\pmb\theta)\cdot f(y_2;\pmb\theta) \cdots f(y_n;\pmb\theta)$ 。

注意，使用最大似然法估计，我们相当于就知道了（或者说，假设了） $y$ 的概率密度函数。

这通常被称为“先验分布”

在抽样之前， $\{y_1,y_2, \cdots,y_n \}$ 被视为随机过程；抽样后， $\{y_1,y_2, \cdots,y_n \}$ 就有了特定的样本值。因此，可以将样本的联合密度函数看作在 $\{y_1,y_2, \cdots,y_n \}$ 给定下，关于 $\pmb \theta$ 的函数。于是我们定义似然函数（likelihood function）为：
$L(\pmb \theta;y_1,\cdots,y_n) = \prod_{i=1}^n f(y_i;\pmb \theta)$
由此可见，似然函数与联合密度函数完全相等，只是把原本是参数的 $\pmb \theta$ 换成了自变量；把原来是自变量的 $\{y_1,y_2, \cdots,y_n \}$ 换成了参数。为了运算方便，通常把似然函数取对数，将乘积转换为求和：
$\ln L(\pmb \theta;y_1,\cdots,y_n) = \ln \prod_{i=1}^n f(y_i;\pmb \theta) = \sum_{i=1}^n \ln f(y_i;\pmb \theta)$
最大似然估计法（Maximum Likelihood Estimation，MLE or ML）来源于一个简单而深刻的想法：给定样本取值后，样本最有可能来自于参数 $\pmb\theta$ 为何值的总体。换言之，寻找 $\hat{\pmb \theta}_{ML}$ 使得观察到样本数据的可能性最大，于是我们的目标是最大化对数似然函数（Loglikelihood Function，LLF）：
$\max_{\pmb\theta \in \Theta} \ln L(\pmb \theta;\pmb y)$

说人话，也就是MLE是“存在即真理”的方法论。既然我现在观察到 $\{y_1,y_2, \cdots,y_n \}$ 出现了，所以参数 $\pmb \theta$ 也应该尽可能让 $\{y_1,y_2, \cdots,y_n \}$ 出现的概率最大。

想到学概率论的一个例子：A和B一起射箭，A和B是否射中应该服从一个两点分布，分别设A和B射中的概率为 $p_A$ 和 $p_B$ 。现在观察到A的射箭结果为 $(0,0,0,0,0)$ 而B的为 $(1,0,0,0,1)$ 。于是我们可以从A和B射箭的结果中反推A和B射中的概率为0和0.2。在这个的例子中：

A和B射中的概率就是这个两点分布的参数 $\theta$

射箭结果就是 $\{y_1,\cdots,y_n\}$

之所以是0和0.2，是因为他们是让已经观察到的结果出现的概率最大的参数。

数学上，我们通常把最大似然估计量 $\hat{\pmb \theta}_{ML}$ 写成：
$\hat{\pmb \theta}_{ML} \equiv \mathop{\rm argmax}_{\pmb\theta} \ln L(\pmb \theta;\pmb y)$
假设存在唯一内点解，那么无约束极值问题的一阶条件可以写成：
$s(\pmb \theta;\pmb y) \equiv \frac{\partial \ln L(\pmb \theta; \pmb y)}{\partial \pmb \theta} \equiv \left(\begin{array}{c} \frac{\partial \ln L(\pmb \theta; \pmb y)}{\partial \theta_1} \\ \vdots \\ \frac{\partial \ln L(\pmb \theta; \pmb y)}{\partial \theta_K} \end{array}\right) = \pmb 0$

此一阶条件要求LLF的梯度向量（gradient） $s(\pmb \theta;\pmb y)$ 为 $\pmb 0$ 。这实际上是 $K$ 个位置参数 $(\theta_1,\cdots,\theta_K)$ 和 $K$ 各方程构成的方程组。该梯度向量也被称为得分函数（score function）或得分向量（score vector）。 $s(\pmb \theta;\pmb y)$ 本身是 $\pmb y$ 的函数，从而也是一个随机变量。记真实的参数为 $\pmb\theta_0$ ，那么我们有：

命题1 得分函数的期望为 $\pmb 0$ 。如果似然函数正确（correctly specified），则 ${\rm E}[s(\pmb\theta_0;\pmb y)]=\pmb 0$

证明1：如果似然函数正确， ${\rm E}[s(\pmb\theta_0;\pmb y)]=\pmb 0$

因为似然函数 $L(\pmb \theta;\pmb y)$ 是概率密度函数，所以 $-\infty \rightarrow +\infty$ 的积分为1，即：
$\int L(\pmb \theta;\pmb y) {\rm d}\pmb y = 1$
为了引入LLF，把上面的积分写成：
$\int \exp[\ln L(\pmb \theta;\pmb y)] {\rm d}\pmb y = 1$
两边对 $\pmb \theta$ 求导，有：
$\int \exp[\ln L(\pmb \theta;\pmb y)] \cdot \frac{\partial \ln L(\pmb \theta;\pmb y)}{\partial \pmb \theta} {\rm d}\pmb y = \pmb 0$
也就是：
$\int L(\pmb \theta;\pmb y) \cdot \frac{\partial \ln L(\pmb \theta;\pmb y)}{\partial \pmb \theta} {\rm d}\pmb y = \pmb 0$
由于似然函数正确，于是在 $\pmb \theta = \pmb \theta_0$ 处， $L(\pmb \theta;\pmb y)$ 就是真实的概率密度函数。于是令上式 $\pmb \theta = \pmb \theta_0$ ，就有：
${\rm E}\left[ \frac{\partial \ln L(\pmb \theta;\pmb y)}{\partial \pmb \theta} \right] = \int L(\pmb \theta;\pmb y) \cdot \frac{\partial \ln L(\pmb \theta;\pmb y)}{\partial \pmb \theta} {\rm d}\pmb y = \pmb 0$
中间的等号用了期望的定义： ${\rm E}[g(x)] = \int g(x)f(x){\rm d}x$ ，其中 $f(x)$ 为概率密度函数。

证毕。

进一步，可以把得分函数 $s(\pmb \theta;\pmb y)$ 写成：
$s(\pmb \theta;\pmb y) \equiv \frac{\partial \ln L(\pmb \theta; \pmb y)}{\partial \pmb \theta} = \frac{\partial \sum\limits_{i=1}^n \ln f(y_i;\pmb \theta)}{\partial \pmb\theta} = \sum_{i=1}^n \frac{\partial f(y_i;\pmb \theta)}{\partial \pmb\theta} \equiv \sum_{i=1}^n s_i(\pmb \theta;y_i)$

其中， $s_i(\pmb \theta;y_i) \equiv \frac{\partial f(y_i;\pmb \theta)}{\partial \pmb\theta}$ 为第 $i$ 个观测值对得分函数的贡献。

6.2 线性回归模型的MLE

假设线性回归模型为：

$\pmb y = {\bf X}\pmb\beta+\pmb\varepsilon$

为了使用MLE，首先要对扰动线做一些假设（球形扰动项）： $\pmb\varepsilon|{\rm X}\sim N(\pmb 0,\sigma^2 {\bf I}_n)$ ，于是被解释变量的条件分布为 $\pmb y|{\rm X} \sim N({\bf X}\pmb\beta,\sigma^2 {\bf I}_n)$ ，于是其条件概率密度函数为：

$f(\boldsymbol{y} | \bf{X})=\left(2 \pi \sigma^{2}\right)^{-n / 2} \exp \left\{-\frac{1}{2 \sigma^{2}}(\boldsymbol{y}-\bf{X} \boldsymbol{\beta})^{\prime}(\boldsymbol{y}-\bf{X} \boldsymbol{\beta})\right\}$

用 $\hat{\boldsymbol{\beta}}$ 和 $\hat{\sigma}^2$ 代入，取对数，于是LLF为：

$\ln L(\hat{\boldsymbol{\beta}},\hat{\sigma}^2) = -\frac{n}{2} \ln 2\pi - \frac{n}{2} \ln \hat{\sigma}^2 - \frac{1}{2\hat{\sigma}^2}(\boldsymbol y - {\bf X}\hat{\boldsymbol \beta})^\prime(\boldsymbol y - {\bf X}\hat{\boldsymbol \beta})$

极大似然估计就是要求 $\hat{\boldsymbol{\beta}}$ 和 $\hat{\sigma}^2$ 使得 $\ln L(\hat{\boldsymbol{\beta}},\hat{\sigma}^2)$ 最大。由于 $\hat{\boldsymbol{\beta}}$ 只出现在第三项，于是 $\hat{\boldsymbol{\beta}}$ 要使得 $(\boldsymbol y - {\bf X}\hat{\boldsymbol \beta})^\prime(\boldsymbol y - {\bf X}\hat{\boldsymbol \beta})$ 最小，这正正好是OLS的目标函数 $\boldsymbol{e'e}$ ，于是：
$\hat{\boldsymbol{\beta}}_{ML} = \hat{\boldsymbol{\beta}}_{OLS} = ({\bf X'X})^{-1}{\bf X'}\boldsymbol y$
下一步就是求解 $\hat{\sigma}^2$ ，此时LLF变为集中对数似然函数（concentrated log likelihood function）：
$\ln L(\hat{\boldsymbol{\beta}},\hat{\sigma}^2) = -\frac{n}{2} \ln 2\pi - \frac{n}{2} \ln \hat{\sigma}^2 - \frac{1}{2\hat{\sigma}^2} \boldsymbol{e'e}$
对 $\hat{\sigma}^2$ 求导，有：
$-\frac{n}{2} \frac{1}{\tilde{\sigma}^{2}}+\frac{1}{2 \tilde{\sigma}^{4}} \boldsymbol{e'e}=0$
于是可以求解：
$\hat{\sigma}_{ML}^{2}=\frac{\boldsymbol{e^{\prime} e}}{n} \neq \hat{\sigma}_{OLS}^{2}=\frac{\boldsymbol{e^{\prime} e}}{n-K}=s^{2}$
于是我们发现，MLE对回归系数 $\boldsymbol \beta$ 的估计与OLS是一致的，不过对扰动项方差 $\sigma^2$ 的估计则不同（受制于自由度的问题），不过这个问题在大样本下消失。另外，由于我们知道 $s^2$ 是 $\sigma^2$ 的无偏估计，从而 $\hat{\sigma}^2_{ML}$ 对 $\sigma^2$ 是有偏估计。

6.3 MLE的数值解

如果模型存在非线性，MLE通常没有解析解，只能寻找数值解（numerical solution），不过这部分不是我们的重点，可以自己翻看教材。

实际上Matlab、Python、R之类的都提供了非线性求解器，所以这不需要我们操心

6.4 信息矩阵与无偏估计的最小方差

为了研究MLE的大样本性质，定义信息矩阵（information matrix）为对数似然函数的黑塞矩阵之期望值的负数，即：

在一维情形下，就是LLF的二阶导数的负数。由于LLF为凹函数，所以其二阶导数为负数，于是一维情形下的信息“矩阵”为正。更一般地，就是LLF在空间的高维曲率（curvature），取期望以后的其实就是表示平均曲率（对平均）。

在求解时，如果LLF曲率大、很陡峭，就比较容易分辨的位置；如果LLF曲率小、很平坦，就很难分辨的位置，如图6.3。更极端地，如果LLF是完全平坦的，那么MLE没有唯一解。

image

信息矩阵涉及二阶偏导数，不易计算，所以通常把它表达为一阶偏导数的乘积的形式：

命题2 在 $\boldsymbol \theta = \boldsymbol \theta_0$ ，信息矩阵等式（information matrix equality）成立：
$\boldsymbol{I}\left(\boldsymbol{\theta}_{0}\right)=-\mathrm{E}\left[\frac{\partial^{2} \ln L\left(\boldsymbol{\theta}_{0} ; \boldsymbol{y}\right)}{\partial \boldsymbol{\theta} \partial \boldsymbol{\theta}^{\prime}}\right]=\mathrm{E}\left[\frac{\partial \ln L\left(\boldsymbol{\theta}_{0} ; \boldsymbol{y}\right) }{\partial \boldsymbol{\theta}} \cdot \frac{\partial \ln L\left(\boldsymbol{\theta}_{0} ; \boldsymbol{y}\right)}{\partial \boldsymbol{\theta}^\prime}\right]=\mathrm{E}\left[\boldsymbol{s}\left(\boldsymbol{\theta}_{0} ; \boldsymbol{y}\right) \boldsymbol{s}\left(\boldsymbol{\theta}_{0} ; \boldsymbol{y}\right)^{\prime}\right]$

证明2：信息矩阵等式

从证明1我们已经知道：
$\int \exp[\ln L(\pmb \theta;\pmb y)] \cdot \frac{\partial \ln L(\pmb \theta;\pmb y)}{\partial \pmb \theta} {\rm d}\pmb y = \pmb 0$
该方程两边继续对 $\boldsymbol \theta^\prime$ 求导，有：
$\int\left\{\frac{\partial^{2} \ln L(\boldsymbol{\theta} ; \boldsymbol{y})}{\partial \boldsymbol{\theta} \partial \boldsymbol{\theta}^{\prime}} \exp [\ln L(\boldsymbol{\theta} ; \boldsymbol{y})]+\frac{\partial \ln L(\boldsymbol{\theta} ; \boldsymbol{y})}{\partial \boldsymbol{\theta}} \exp [\ln L(\boldsymbol{\theta} ; \boldsymbol{y})] \frac{\partial \ln L(\boldsymbol{\theta} ; \boldsymbol{y})}{\partial \boldsymbol{\theta}^{\prime}}\right\} \mathrm{d} \boldsymbol{y}=\mathbf{0}$
移项，然后令 $L(\boldsymbol{\theta} ; \boldsymbol{y}) := \exp [\ln L(\boldsymbol{\theta} ; \boldsymbol{y})]$ ，那么就有：
$-\int\left\{\frac{\partial^{2} \ln L(\boldsymbol{\theta} ; \boldsymbol{y})}{\partial \boldsymbol{\theta} \partial \boldsymbol{\theta}^{\prime}} L(\boldsymbol{\theta} ; \boldsymbol{y})\right\} \mathrm{d} \boldsymbol{y}=\int\left\{\frac{\partial \ln L(\boldsymbol{\theta} ; \boldsymbol{y})}{\partial \boldsymbol{\theta}} \cdot \frac{\partial \ln L(\boldsymbol{\theta} ; \boldsymbol{y})}{\partial \boldsymbol{\theta}^{\prime}} L(\boldsymbol{\theta} ; \boldsymbol{y})\right\} \mathrm{d} \boldsymbol{y}$
于是，根据期望的定义，就有：
$-\mathrm{E}\left[\frac{\partial^{2} \ln L\left(\boldsymbol{\theta}_{0} ; \boldsymbol{y}\right)}{\partial \boldsymbol{\theta} \partial \boldsymbol{\theta}^{\prime}}\right]=\mathrm{E}\left[\frac{\partial \ln L\left(\boldsymbol{\theta}_{0} ; \boldsymbol{y}\right)}{\partial \boldsymbol{\theta}} \cdot \frac{\partial \ln L\left(\boldsymbol{\theta}_{0} ; \boldsymbol{y}\right)}{\partial \boldsymbol{\theta}^{\prime}}\right]$
证毕。

其实证明1和证明2的思路都挺简单的，最关键的技巧是从LF写成 $\exp(\ln(\cdot))$ 的形式以凑出LLF。然后都是求导（或者二阶导），交换积分号和求导符号，移项，把LLF重新写成LF的形式，用期望的定义，即可证出来。

命题3** 在 $\boldsymbol \theta = \boldsymbol \theta_0$ ，信息矩阵 $\boldsymbol I(\boldsymbol\theta_0)$ 就是得分函数的协方差矩阵 ${\rm Var}[s(\boldsymbol\theta_0;\boldsymbol y)]$ 。

证明3：得分函数的方差为信息矩阵
$\begin{split} {\rm Var}[s(\boldsymbol\theta_0;\boldsymbol y)] &= {\rm E}\left[s(\boldsymbol\theta_0;\boldsymbol y)s(\boldsymbol\theta_0;\boldsymbol y)^\prime\right] - {\rm E}[s(\boldsymbol\theta_0;\boldsymbol y)] \cdot {\rm E}[s(\boldsymbol\theta_0;\boldsymbol y)]^\prime\\ &= {\rm E}\left[s(\boldsymbol\theta_0;\boldsymbol y)s(\boldsymbol\theta_0;\boldsymbol y)^\prime\right] -0\cdot0\\ &=\boldsymbol I(\boldsymbol\theta_0) \end{split}$
证毕。

在统计学中有一个著名的结论：假设 $\hat{\boldsymbol\theta}$ 是对真实参数 $\boldsymbol\theta_0$ 的任意无偏估计，那么在一定的正则条件（regularity conditions）下， $\hat{ \boldsymbol \theta}$ 的方差不会小于 $[\boldsymbol I(\boldsymbol \theta_0)]^{-1}$ ，即 ${\rm Var}\left(\hat{\boldsymbol \theta}\right) \geqslant [\mathrm I(\boldsymbol \theta_0)]^{-1}$ 。其中 $[\boldsymbol I(\boldsymbol \theta_0)]^{-1}$ 被称为克莱默-劳下限（Cramer-Rao Lower Bound）。因此，无偏估计所能达到的最小方差与信息矩阵（即LLF的平均曲率）有关——曲率 $\boldsymbol I(\boldsymbol \theta_0)$ 越大，则 $[\boldsymbol I(\boldsymbol \theta_0)]^{-1}$ 越小，那么无偏估计可能达到的方差就越小。

证明*：Cramer-Rao Lower Bound

为简单期间，只证明一维情形。由于 $\hat{\boldsymbol \theta}$ 是对真实参数 $\boldsymbol \theta_0$ 的任意无偏估计，于是：
$\boldsymbol \theta_0 = {\rm E}\left[ \hat{\boldsymbol \theta}(\boldsymbol y) \right] = \int \hat{\boldsymbol \theta}(\boldsymbol y) L(\boldsymbol \theta_0;\boldsymbol y) \mathrm d \boldsymbol y$
将上式两边同时对 $\boldsymbol \theta_0$ 求导可得：
$\begin{split} 1 = \int \hat{\boldsymbol \theta}(\boldsymbol y) \frac{\partial L(\boldsymbol \theta_0;\boldsymbol y) }{\partial \boldsymbol\theta}\mathrm d \boldsymbol y &= \int \hat{\boldsymbol \theta}(\boldsymbol y) \frac{\partial \ln L(\boldsymbol \theta_0;\boldsymbol y) }{\partial \boldsymbol\theta} \cdot L(\boldsymbol \theta_0;\boldsymbol y) \mathrm d \boldsymbol y \\ (期望的定义)& = \mathrm{E}\left[\boldsymbol{\hat{\theta}}(\boldsymbol{y}) \frac{\partial \ln L\left(\boldsymbol{\theta}_{0} ; \boldsymbol{y}\right)}{\partial \boldsymbol{\theta}}\right]\\ ({\rm E}(AB)的计算+命题1)&= \operatorname{Cov}\left[\boldsymbol{\hat{\theta}}(\boldsymbol{y}) , \frac{\partial \operatorname{ln} L\left(\boldsymbol\theta_{0} ;\boldsymbol y\right)}{\partial\boldsymbol \theta}\right]\\ (相关系数\leqslant1)&\leqslant \operatorname{Var}\left[\boldsymbol{\hat{\theta}}(\boldsymbol{y}) \right] \cdot \operatorname{Var}\left[\frac{\partial\operatorname{ln} L\left(\boldsymbol\theta_{0} ; \boldsymbol{y}\right)}{\partial \boldsymbol\theta}\right]\\ (信息矩阵等式)&=\operatorname{Var}\left[\boldsymbol{\hat{\theta}}(\boldsymbol{y}) \right] \cdot \mathbf{E}\left[\frac{\partial \ln L\left(\boldsymbol{\theta}_{0} ; \boldsymbol{y}\right)}{\partial\boldsymbol \theta}\right]^{2} \end{split}$
于是：
$\operatorname{Var}\left[\boldsymbol{\hat{\theta}}(\boldsymbol{y}) \right] \geqslant\left[E\left(\frac{\partial \ln L\left(\boldsymbol\theta_{0} ;\boldsymbol y\right)}{\partial \boldsymbol\theta}\right)^{2}\right]^{-1}=[\boldsymbol I\left(\boldsymbol\theta_{0}\right)]^{-1}$
证毕。

可以证明，在古典线性回归模型中，根据信息矩阵的定义：
$[\boldsymbol I\left(\boldsymbol\theta_{0}\right)]^{-1} = \left(\begin{array}{c} \sigma^2({\bf X'X})^{-1} & \boldsymbol 0\\ \boldsymbol 0 & \frac{2\sigma^4}{n} \end{array}\right)$

其中， $\boldsymbol\theta_0 = \left(\boldsymbol\beta \sigma^2\right)^\prime$ 。由于 ${\rm Var}(\hat{\boldsymbol\beta}_{OLS}) = \sigma^2 ({\bf X'X})^{-1}$ ，于是 $\hat{\boldsymbol\beta}_{ML} = \hat{\boldsymbol\beta}_{OLS}$ 均达到了无偏估计的最小方差。

命题4 在高斯-马尔可夫定理中，如果加上扰动项为正态分布的假定，那么OLS是达到了最佳无偏估计而非仅仅最佳线性无偏估计

高级计量经济学 11：最大似然估计(上)

高级计量经济学 11：最大似然估计(上)

6 最大似然估计法

6.1 最大似然估计法的定义

6.2 线性回归模型的MLE

6.3 MLE的数值解

6.4 信息矩阵与无偏估计的最小方差

友情链接更多精彩内容