background-image: url("../pic/slide-front-page.jpg") class: center,middle # 高级计量暑期班</br></br>(Seminar of Advanced Econometrics) <!--- chakra: libs/remark-latest.min.js ---> ### 胡华平 ### 西北农林科技大学 ### 经济管理学院数量经济教研室 ### huhuaping01@hotmail.com ### 2022-06-28
--- class: center, middle, duke-green,hide_logo name:chapter-RDD01 # RDD PART 01:非参数估计 ### [1.均值估计(Means Estimator)](#mean-est) ### [2.局部回归(Local Regression)](#local-reg) ### [3.估计效果(Performance Analysis)](#performance) ### [4.群组分析(Cluster observations)](#cluster) --- layout: false class: center, middle, duke-orange,hide_logo name: mean-est # 1.均值估计(Means Estimator) ### [1.1 箱组均值估计(Binned Means Estimator)](#me-binned) ### [1.2 滚动箱组均值估计(Rolling Binned Means Estimator)](#me-rolling) ### [1.3 核估计(Kernel Estimator)](#me-kernel) --- layout: true <div class="my-header-h2"></div> <div class="watermark1"></div> <div class="watermark2"></div> <div class="watermark3"></div> <div class="my-footer"><span><a href="https://home.huhuaping.com/">https://home.huhuaping.com </a>    <a href="#chapter-RDD01"> RDD Part01 非参数估计 </a>                     <a href="#mean-est"> 1.均值估计(Means Estimator) </a> </span></div> --- ### (引子)为什么要进行均值估计? 对于计量模型: `$$\begin{align} Y &= m(X) +e \end{align}$$` 研究者首先需要关注的是 - **条件期望函数**(Conditional Expectation Function ,CEF): `$$\begin{align} \mathbb{E}[Y|X=x] \equiv m(x) \end{align}$$` - 此时: `$$\begin{align} Y & = m(X) +e \\ & = \mathbb{E}[Y|X=x] + e \end{align}$$` --- ### (引子)什么是非参数估计? - 理论上,条件期望函数 `\(m(x)\)`可以表现为明确的**参数化**形式(parametric function),也可以表现为任意的**非参数化**形式(non-parametric function)。 .case[ - 常见的参数化条件期望函数,例如线性形式: `$$Y= m(x) +e = \beta_0 +\beta_1 X +e$$` ] --- ### (引子)什么是非参数估计? - **非参数回归模型**(nonparametric regression model):假定**条件期望函数**表现为任意的**非参数化**形式的回归模型。 .case[ - **非参数回归模型**可以表达为: `$$\begin{align} Y = \mathbb{E}(Y|X=x) + e &=m(x) +e \\ \mathbb{E}(e|X=x) &= 0 \\ \mathbb{E}(e^2|X=x) &= \sigma^2(X) \end{align}$$` ] - 此时,我们的目标就是估计得到条件期望函数 `\(\widehat{m}(x)\)`。 --- exclude: true ## Code script: Hensen bruce fig 19.1-a - 此代码仅用于生成本地图片,为提高渲染效率代码文件仅需运行一次。 - 如更新代码,请更改`eval = TRUE` --- ### (示例)模拟数据集 为了更好地进行数据验证,我们将根据如下规则生成蒙特卡洛模拟数据集: `$$\begin{align} Y_i &= m(X) +e_i =\frac{sin(\frac{\pi}{4}\cdot(X_i-2))}{\frac{\pi}{4}\cdot(X_i-2)} +e_i\\ X_i &\sim U(0,10)\\ e_i &\sim N(0, 2)\\ n &=100 \end{align}$$` - 此时,我们具有**上帝视角**,实际上已经知道**数据生成机制**(DGP) - 此时,我们心里面已知真实模型为**非线性的** --- ###(示例)模拟的样本数据集 .pull-left[
] .pull-right[ - 样本数据的描述性统计如下: ``` index X Y Min. : 1 Min. :0.1 Min. :-0.65 1st Qu.: 26 1st Qu.:2.3 1st Qu.:-0.06 Median : 50 Median :4.6 Median : 0.27 Mean : 50 Mean :4.9 Mean : 0.38 3rd Qu.: 75 3rd Qu.:7.5 3rd Qu.: 0.81 Max. :100 Max. :9.9 Max. : 1.46 ``` ] --- ### (示例)样本数据散点图 <img src="seminar-RDD-part01_files/figure-html/unnamed-chunk-5-1.svg" height="450" style="display: block; margin: auto;" /> --- exclude: true ## Code chunk: 生成箱组区隔及相关计算 --- name: me-binned ### 1.1 箱组均值估计:表达式 对于**非参数模型**: `$$\begin{align} Y = \mathbb{E}(Y|X=x) + e &=m(x) +e \end{align}$$` 我们可以直接把数据集划分为不同**箱组**(bins),然后简单地计算各个箱组中 `\(Y_i\)`的均值。 `$$\begin{align} \hat{m}(x)=\frac{\sum_{i=1}^{n} \mathbb{1}\left\{\left|X_{i}-x\right| \leq h\right\} \cdot Y_{i}}{\sum_{i=1}^{n} \mathbb{1}\left\{\left|X_{i}-x\right| \leq h\right\}} \end{align}$$` 其中: - `\(\mathbb{1}\left\{\left|X_{i}-x\right| \leq h\right\}\)`为**指示函数**,取值为 `\(\{0,1\}\)`,以表明 `\(X_i\)`是否落在特定**箱组**内 - 以上公式可以简单视作为箱组内的**简单算数平均数**公式 --- ### 1.1 箱组均值估计:操作步骤 **箱组均值估计**(Binned Estimator)的操作步骤如下: - 根据计算点 `\(X=x_j\)`,按照特定谱宽 `\(h\)`,划分出若干**箱组**(bins): `\(\{b_1,b_2,\cdots,b_q\}\)` $$b_{j}=[x_j-h,x_j +h] $$ - 根据样本数据集,以及 `\(X_i\)`的实际情况,确定数据对 `\(\{X_i, Y_i\}\)`的箱组归属: `$$\mathbb{1}\left\{\left|X_{i}-x_j\right| \leq h\right\}$$` - 最后计算不同箱组的 `\(Y_i\)`的均值 `\(\widehat{m}(x_{b_j}),j \in (1,2,\cdots,q)\)`: `$$\begin{align} \hat{m}(x)=\frac{\sum_{i=1}^{n} \mathbb{1}\left\{\left|X_{i}-x\right| \leq h\right\} \cdot Y_{i}}{\sum_{i=1}^{n} \mathbb{1}\left\{\left|X_{i}-x\right| \leq h\right\}} \end{align}$$` --- ### (示例)箱组均值估计:设定箱组划分规则 下面我们分别设定如下**箱组**划分规则: - 设定箱组取值中心点 `\((X=x_i),x_i \in (1, 3, 5, 7, 9)\)`,以及谱宽 `\(h=1\)` - 然后得到箱组区块bins=[0,2)、[2,4)、[4,6)、[6,8)、[8,10],箱组数为5。 --- ### (示例)箱组均值估计:数据计算表 利用箱组均值估计公式,我们可以计算得到不同箱组的均值估计:
- 箱组内因变量观测值的求和`sum_ky` `\(=\sum_{i=1}^{n} \mathbb{1}\left\{\left|X_{i}-x\right| \leq h\right\} \cdot Y_{i}\)` - 箱组的样本数`sum_k` `\(=\sum_{i=1}^{n} \mathbb{1}\left\{\left|X_{i}-x\right| \leq h\right\}\)` - 箱组的均值估计`m0` `\(=\widehat{m}(x_j),j \in (1,2,\cdots,5)\)` --- ### (示例)箱组均值估计:图形表达1/3 - 首先我们展示的是5个箱组的划分: <img src="seminar-RDD-part01_files/figure-html/unnamed-chunk-8-1.svg" height="450" style="display: block; margin: auto;" /> .footnote[ **说明**:a)垂直虚线表示箱组中心取值点 `\(x_j\)`;b)不同矩形颜色区块表示不同箱组。 ] --- ### (示例)箱组均值估计:图形表达2/3 - 根据箱组均值计算值,我们展示在散点图中: <img src="seminar-RDD-part01_files/figure-html/unnamed-chunk-9-1.svg" height="450" style="display: block; margin: auto;" /> --- ### (示例)箱组均值估计:图形表达3/3 - 简单地,可将**箱组均值**作为对这一箱组**条件期望函数**CEF的近似: <img src="seminar-RDD-part01_files/figure-html/unnamed-chunk-10-1.svg" height="450" style="display: block; margin: auto;" /> --- name: me-rolling ### 1.2 滚动箱组均值估计:定义及表达式 **滚动箱组均值估**(The rolling binned means estimator):以系列数值 `\(x\)`为中心,以 `\(h\)`为谱宽,**滚动**构建一系列箱组(箱组会有重叠),并分别计算出系列箱组的均值。 `$$\begin{align} \hat{m}(x)=\frac{\sum_{i=1}^{n} \mathbb{1}\left\{\left|X_{i}-x\right| \leq h\right\} \cdot Y_{i}}{\sum_{i=1}^{n} \mathbb{1}\left\{\left|X_{i}-x\right| \leq h\right\}} \end{align}$$` .fyi[ - 我们后面马上会介绍,**滚动箱组均值估**实际上是一类特殊的**核估计**(kernel)情形,具体为Nadaraya-Watson 矩形和估计(NW rectangular kernel estimator)。 ] --- ### 1.2 滚动箱组均值估计:操作过程 **滚动箱组均值估计**(Rolling Binned Estimator)的操作步骤如下: - 根据计算点 `\(X=x_j\)`,按照特定谱宽 `\(h\)`,划分出若干**箱组**(bins)(箱组会有重叠): `\(\{b_1,b_2,\cdots,b_q\}\)` $$b_{j}=[x_j-h,x_j +h] $$ - 根据样本数据集,以及 `\(X_i\)`的实际情况,确定数据对 `\(\{X_i, Y_i\}\)`的箱组归属: `$$\mathbb{1}\left\{\left|X_{i}-x_j\right| \leq h\right\}$$` - 最后计算不同箱组的 `\(Y_i\)`的均值 `\(\widehat{m}(x_{b_j}),j \in (1,2,\cdots,q)\)`: `$$\begin{align} \hat{m}(x)=\frac{\sum_{i=1}^{n} \mathbb{1}\left\{\left|X_{i}-x\right| \leq h\right\} \cdot Y_{i}}{\sum_{i=1}^{n} \mathbb{1}\left\{\left|X_{i}-x\right| \leq h\right\}} \end{align}$$` --- ### (示例)滚动箱组均值估计:设定箱组划分规则 下面我们分别设定如下**箱组**划分规则: - 设定箱组取值中心点 `\((X=x_i),x_i \in (0.00, 0.01, 0.02, 0.03, 0.04, 0.05, \cdots, 9.95, 9.96, 9.97, 9.98, 9.99, 10.00)\)`,以及谱宽<sup>a</sup> `\(h=1\)` - 那么,可以得到箱组区块bins=[-1,1)、[-0.99,1.01)、[-0.98,1.02) `\(\ldots\)`[8.98,10.98)、[8.99,10.99)、[9,11),箱组数1001。 .footnote[ <sup>a</sup> 此时滚动箱组均值估计等价于谱宽 `\(h=1\)`的**矩形核函数**(rectangular kernel)估计 ] --- ### (示例)滚动箱组均值估计:数据计算表 利用箱组均值估计公式,我们可以计算得到不同箱组的均值估计:
- 箱组内因变量观测值的求和`sum_ky` `\(=\sum_{i=1}^{n} \mathbb{1}\left\{\left|X_{i}-x\right| \leq h\right\} \cdot Y_{i}\)` - 箱组的样本数`sum_k` `\(=\sum_{i=1}^{n} \mathbb{1}\left\{\left|X_{i}-x\right| \leq h\right\}\)` - 箱组的均值估计`m1` `\(=\widehat{m}(x_j),j \in (1,2,\cdots,1001)\)` --- ### (示例)滚动箱组均值估计:图形表达 - 根据箱组均值计算值,我们展示在散点图中: <img src="seminar-RDD-part01_files/figure-html/unnamed-chunk-12-1.svg" height="450" style="display: block; margin: auto;" /> --- ### (示例)滚动箱组均值估计:图形表达 - 同样地,可将**箱组均值**作为对这一箱组**条件期望函数**CEF的近似: <img src="seminar-RDD-part01_files/figure-html/unnamed-chunk-13-1.svg" height="450" style="display: block; margin: auto;" /> --- name: me-kenerl ### 1.3 核估计:回顾与思考 上述两种箱组均值估计的公式中: `$$\begin{align} \hat{m}(x)=\frac{\sum_{i=1}^{n} \mathbb{1}\left\{\left|X_{i}-x\right| \leq h\right\} \cdot Y_{i}}{\sum_{i=1}^{n} \mathbb{1}\left\{\left|X_{i}-x\right| \leq h\right\}} \end{align}$$` - 对条件期望的估计 `\(\widehat{m}(x)\)`结果都呈现一定的**锯齿**形态(jagged),也即估计结果不太**平滑**(smoothed)。 --- ### 1.3 核估计:回顾与思考 .fyi[ **思考**: - 为什么估计结果会呈现不太**平滑**的**锯齿**形态? - 能不能让估计结果更加平滑呢? ] -- .notes[ **回答**: - 问题的关键在于上面的估计公式使用了箱组**指示**函数(可视作为权重) `\(\sum_{i=1}^{n} \mathbb{1}\left\{\left|X_{i}-x\right| \leq h\right\}\)`,而这个**权重函数**本身是**跳跃的**。 - 有没有一种办法能够基于平滑的**权重函数**来计算CEF的估计值呢? ] --- ### 1.3 核估计:概念 **核估计**(kernel estimator):基于多种类型的**核函数**(kernel function)作为**权重函数**——可以是**连续的**,也可以是**跳跃的**——来估计条件期望函数 `\(\widehat{m}(x)\)`的一种估计方法。 `$$\begin{align} \hat{m}_{\mathrm{nw}}(x)=\frac{\sum_{i=1}^{n} K\left(\frac{X_{i}-x}{h}\right) Y_{i}}{\sum_{i=1}^{n} K\left(\frac{X_{i}-x}{h}\right)} \end{align}$$` - 其中: `\(K(u)\)`为核函数(kernel function) --- ### 1.3 核估计:特点 `$$\begin{align} \hat{m}_{\mathrm{nw}}(x)=\frac{\sum_{i=1}^{n} K\left(\frac{X_{i}-x}{h}\right) Y_{i}}{\sum_{i=1}^{n} K\left(\frac{X_{i}-x}{h}\right)} \end{align}$$` - 尽管这里使用了核函数来计算权重,但是本质上上述公式还是利用了**箱组估计**的方法,也即把箱组计算值作为这一箱组的**代表值**。 - 以上估计方法也被称为**局部常数估计**(local constant estimator)或**NW估计**(Nadaraya-Watson estimator) - 容易发现前述**箱组均值估计**和**滚动箱组均值估计**都是**局部常数估计**(local constant estimator)的两个特例。 --- ### 1.3 核估计:核函数定义 **定义**:若满足如下条件,则可称之为**核函数**(Kernel function) `\(K(u)\)` - `\(0 \leq K(u) \leq \bar{K}<\infty\)`, - `\(K(u)=K(-u)\)`, - `\(\int_{-\infty}^{\infty} K(u) d u=1\)`, - 对所有的正整数 `\(r\)`都有 `\(\int_{-\infty}^{\infty}|u|^{r} K(u) d u<\infty\)` **定义**:**正规化核函数**(normalized kernel function)需满足 `$$\int_{-\infty}^{\infty} u^2K(u) d u=1$$` .notes[ - **核函数**本质上是一种**边界约束的**概率密度函数(bounded pdf) - **核函数**是原点对称的,且核函数为非负数(因此可用于权重) ] --- exclude: true ## Code script: Hensen bruce prob21 fig 17-2 - 此代码仅用于生成本地图片,为提高渲染效率代码文件仅需运行一次。 - 如更新代码,请更改`eval = TRUE` --- ### 1.3 核估计:常见的正规化核函数1/2 - **矩形核函数**(Rectangular Kernel), `\(R_K= \frac{1}{2\sqrt{3}}\)` `$$\begin{aligned} &K(u)=\left\{\begin{array}{cc} \frac{1}{2 \sqrt{3}} & \text { if }|u|<\sqrt{3} \\ 0 & \text { otherwise } \end{array}\right.\\ \end{aligned}$$` - **高斯核函数**(Gaussian Kernel), `\(R_K= \frac{1}{2\sqrt{\pi}}\)` `$$\begin{aligned} &K(u)=\frac{1}{\sqrt{2 \pi}} \exp \left(-\frac{u^{2}}{2}\right)\\ \end{aligned}$$` - **叶氏核函数**(Epanechnikov Kernel), `\(R_K= \frac{3\sqrt{5}}{25}\)` `$$\begin{aligned} &K(u)=\left\{\begin{array}{cl} \frac{3}{4 \sqrt{5}}\left(1-\frac{u^{2}}{5}\right) & \text { if }|u|<\sqrt{5} \\ 0 & \text { otherwise } \end{array}\right. \end{aligned}$$` --- ### (示例)核函数类型1 <img src="seminar-RDD-part01_files/figure-html/unnamed-chunk-15-1.svg" height="480" style="display: block; margin: auto;" /> --- ### 1.3 核估计:常见的正规化核函数2/2 - **三角核函数**(Triangular Kernel), `\(R_K= \frac{\sqrt{3}}{9}\)` `$$\begin{aligned} &K(u)=\left\{\begin{array}{cl} \frac{1}{\sqrt{6}}\left(1-\frac{|u|}{\sqrt{6}}\right) & \text { if }|u|<\sqrt{6} \\ 0 & \text { otherwise } \end{array}\right. \end{aligned}$$` - **双权核函数**(Biweight Kernel), `\(R_K= \frac{5\sqrt{7}}{49}\)` `$$\begin{aligned} K(u)=\left\{\begin{array}{cl} \frac{15}{16 \sqrt{7}}\left(1-\frac{u^{2}}{7}\right) & \text { if }|u|<\sqrt{7} \\ 0 & \text { otherwise } \end{array}\right. \end{aligned}$$` --- ### (示例)核函数类型(续) <img src="seminar-RDD-part01_files/figure-html/unnamed-chunk-16-1.svg" height="480" style="display: block; margin: auto;" /> --- ### 1.3 核估计:操作过程 **NW核估计**(NW Estimator)的操作步骤如下: - 根据计算点 `\(X=x_j\)`,按照特定谱宽 `\(h\)`,划分出若干**箱组**(bins)(箱组会有重叠): `\(\{b_1,b_2,\cdots,b_q\}\)` $$b_{j}=[x_j-h,x_j +h] $$ - 根据样本数据集,以及 `\(X_i\)`的实际情况,计算特定**核函数化**的权重值: `$$K\left(\frac{X_{i}-x_j}{h}\right)$$` - 最后计算不同箱组的 `\(Y_i\)`的均值 `\(\widehat{m}(x_{b_j}),j \in (1,2,\cdots,q)\)`: `$$\begin{align} \hat{m}_{\mathrm{nw}}(x)=\frac{\sum_{i=1}^{n} K\left(\frac{X_{i}-x}{h}\right) Y_{i}}{\sum_{i=1}^{n} K\left(\frac{X_{i}-x}{h}\right)} \end{align}$$` --- ### (示例)核估计:设定箱组划分规则 下面我们分别设定如下**箱组**划分规则: - 设定箱组取值中心点 `\((X=x_i),x_i \in (0.00, 0.01, 0.02, 0.03, 0.04, 0.05, \cdots, 9.95, 9.96, 9.97, 9.98, 9.99, 10.00)\)` - 确定核函数类型为**高斯核函数**<sup>a</sup>,谱宽设定为 `\(h=1/\sqrt{3}\)`。 - 那么,可以得到箱组区块bins=[-1,1)、[-0.99,1.01)、[-0.98,1.02) `\(\ldots\)`[8.98,10.98)、[8.99,10.99)、[9,11),箱组数1001。 --- ### (示例)核估计:数据计算表 利用箱组均值核函数估计公式,我们可以计算得到不同箱组的均值估计:
--- ### (示例)核估计:图形表达 - 同样地,可将**箱组均值**作为对这一箱组**条件期望函数**CEF的近似: <img src="seminar-RDD-part01_files/figure-html/unnamed-chunk-18-1.svg" height="450" style="display: block; margin: auto;" /> --- ### (示例)CEF估计:三种估计方法的图形比较 <img src="seminar-RDD-part01_files/figure-html/unnamed-chunk-19-1.svg" height="450" style="display: block; margin: auto;" /> - **NW核估计**方法相比**箱组均值估计**和**滚动箱组均值估计**要更加平滑。 --- layout: false class: center, middle, duke-orange,hide_logo name: local-reg # 2.局部回归(Local Regression) ### [2.1 箱组线性回归(Binned Regression)](#binned-est) ### [2.2 滚动箱组回归(Rolling Regression)](#rolling-est) ### [2.3 局部线性回归(Local Linear Regression)](#ll-est) ### [2.4 局部多项式回归(Local Polynomial Regression)](#lp-est) --- layout: true <div class="my-header-h2"></div> <div class="watermark1"></div> <div class="watermark2"></div> <div class="watermark3"></div> <div class="my-footer"><span><a href="https://home.huhuaping.com/">https://home.huhuaping.com </a>    <a href="#chapter-RDD01"> RDD Part01 非参数估计 </a>                     <a href="#local-reg"> 2.局部回归(Local Regression) </a> </span></div> --- exclude: true ## Code script: Hensen bruce fig 19.1-b - 此代码仅用于生成本地图片,为提高渲染效率代码文件仅需运行一次。 - 如更新代码,请更改`eval = TRUE` --- ### 引子:从NW估计量说起 在一个局部区域(如 `\(X=x\)`的领域内),对 `\(m(x)\)`的Nadaraya-Watson (NW)估计量将表现为**常函数**(constant function)形态,此时也称为**局部常数估计量**(Local constant estimator)。 - 此时,Nadaraya-Watson (NW)估计量为一种局部近似,也即当在局部渐近取值 `\(X\simeq x\)`时, `\(m(X) \simeq m(x)\)` `$$Y=m(X)+e \simeq m(x)+e$$` 上述模型可以视作为回归方程,我们需要估计得到 `\(\widehat{m}(x)\)`,也即: `$$\begin{align} \widehat{m}_{\mathrm{nw}}(x)=\underset{m}{\operatorname{argmin}} \sum_{i=1}^{n} K\left(\frac{X_{i}-x}{h}\right)\left(Y_{i}-m\right)^{2} \end{align}$$` - 本质上,以上就是一个 `\(Y\)`对截距项的**加权回归**(Weighted regression)估计量。 --- ### 引子:CEF的回归估计方法 事实上,基于**箱组**原理的CEF估计方法,我们都可以使用回归方法来进行估计,具体包括: - 箱组线性回归(binned regression):箱组内直接使用OLS回归估计 - 滚动箱组核回归(NW regression):箱组内直接使用OLS回归估计 - 局部线性回归(Local linear regression):箱组内使用加权OLS回归估计 - 局部多项式回归(Local Polynomial regression):箱组内使用多项式的加权OLS回归估计 --- name: binned-est ### 2.1 箱组线性回归:原理 **箱组线性回归**(binned regression):箱组内直接使用OLS回归估计 `$$Y=m(X)+e \simeq m(x)+e$$` - 上述模型可以视作为回归方程,我们需要估计得到 `\(\widehat{m}(x)\)`,也即: `$$\begin{align} \widehat{m}_{\mathrm{bin}}(x)= \underset{m}{\operatorname{argmin}} \sum_{i=1}^{n} \mathbf{1} \{|X_i-x|\leq h\}\cdot \left(Y_{i}-m\right)^{2} \end{align}$$` - 因为在箱组内的权重都是一样的,因此可以直接在**箱组内**的子样本数据里使用OLS回归估计得到 `$$\widehat{m}(x)={\widehat{Y}_i} \quad \text{(OLS)}$$` --- ### 2.1 箱组线性回归:操作步骤 **箱组线性回归**(Binned regression)的操作步骤如下: - 根据计算点 `\(X=x_j\)`,按照特定谱宽 `\(h\)`,划分出若干**箱组**(bins): `\(\{b_1,b_2,\cdots,b_q\}\)` $$b_{j}=[x_j-h,x_j +h] $$ - 根据样本数据集,以及 `\(X_i\)`的实际情况,确定数据对 `\(\{X_i, Y_i\}\)`的箱组归属(子样本): `$$\mathbb{1}\left\{\left|X_{i}-x_j\right| \leq h\right\}$$` - 然后采用OLS方法计算不同箱组的 `\(Y_i\)`的拟合值 `\(\widehat{Y}_i = \widehat{m}(x) \simeq m(X)\)` `$$\hat{Y}_i = \hat{\alpha}_0 + \hat{\alpha}_1 X_i$$` - 最后,我们可以给定任意值 `\(x_i\)`得到预测的 `\(\widehat{m}(x_i)=\widehat{Y}_i|x_i\)` --- ### (示例)箱组回归估计:设定箱组划分规则 下面我们分别设定如下**箱组**划分规则: - 设定箱组取值中心点 `\((X=x_i),x_i \in (1, 3, 5, 7, 9)\)`,以及谱宽 `\(h=1\)` - 然后得到箱组区块bins=[0,2)、[2,4)、[4,6)、[6,8)、[8,10],箱组数为5。 - 箱组内的线性回归模型为 `$$\hat{Y}_i = \hat{\alpha}_0 + \hat{\alpha}_1 X_i$$` - 基于设定的 `\(x_i \in (0.00, 0.01, 0.02, 0.03, 0.04, 0.05, \cdots, 9.95, 9.96, 9.97, 9.98, 9.99, 10.00)\)`计算预测值 `\(\widehat{m}(x_i)=\widehat{Y}_i|x_i\)`,共有 `\(N=1001\)`个 --- ### (示例)箱组线性回归:数据计算表 `$$\hat{Y}_i = \hat{\alpha}_0 + \hat{\alpha}_1 X_i$$` - 对五个箱组的区块下,我们依次对子样本数据对 `\((Y_i,X_i)\)`进行线性OLS拟合,分别得到截距和斜率: <table> <thead> <tr> <th style="text-align:left;"> par </th> <th style="text-align:right;"> bd1 </th> <th style="text-align:right;"> bd2 </th> <th style="text-align:right;"> bd3 </th> <th style="text-align:right;"> bd4 </th> <th style="text-align:right;"> bd5 </th> </tr> </thead> <tbody> <tr> <td style="text-align:left;"> intercept </td> <td style="text-align:right;"> 0.75 </td> <td style="text-align:right;"> 1.45 </td> <td style="text-align:right;"> 2.32 </td> <td style="text-align:right;"> 0.74 </td> <td style="text-align:right;"> 0.72 </td> </tr> <tr> <td style="text-align:left;"> slope </td> <td style="text-align:right;"> 0.11 </td> <td style="text-align:right;"> -0.19 </td> <td style="text-align:right;"> -0.42 </td> <td style="text-align:right;"> -0.12 </td> <td style="text-align:right;"> -0.09 </td> </tr> </tbody> </table> .footnote[ **说明**:`bd1`、`bd2`等分别代表五个箱组。 ] --- ### (示例)箱组线性回归:数据计算表 利用箱组均值估计公式,我们可以计算得到不同箱组的均值估计:
- 基于设定的`x`表示 `\(x_i \in (0.00, 0.01, 0.02, 0.03, 0.04, 0.05,\cdots)\)`,共有 `\(N=1001\)`个 - 箱组回归估计`m0` `\(\widehat{m}(x_i)=\widehat{Y}_i|x_i\)` --- ### (示例)箱组线性回归:图形表达1/2 - 如前,我们设定区隔了5个箱组 <img src="seminar-RDD-part01_files/figure-html/unnamed-chunk-23-1.svg" height="450" style="display: block; margin: auto;" /> --- ### (示例)箱组线性回归:图形表达2/2 - 根据前面计算表的1001个拟合数据对 `\((x_i,\widehat{m}(x_i))\)`,我们可以得到箱组线性回归估计结果: <img src="seminar-RDD-part01_files/figure-html/unnamed-chunk-24-1.svg" height="430" style="display: block; margin: auto;" /> --- name: rolling-est ### 2.2 滚动箱组回归:原理 **滚动箱组线性回归**(rolling regression):简单地,是通过构建滚动箱组(可能出现箱组重叠状态),然后再对箱组内样本数据直接使用OLS回归估计。此外,容易发现它实际上就是等价于**矩形核函数**的**NW回归**( `\(h=1\)`)。 - 此时,Nadaraya-Watson (NW)估计量为一种局部近似,也即当在局部渐近取值 `\(X\simeq x\)`时, `\(m(X) \simeq m(x)\)` `$$Y=m(X)+e \simeq m(x)+e$$` 上述模型可以视作为回归方程,我们需要估计得到 `\(\widehat{m}(x)\)`,也即: `$$\begin{align} \widehat{m}_{\mathrm{nw}}(x)=\underset{m}{\operatorname{argmin}} \sum_{i=1}^{n} K\left(\frac{X_{i}-x}{h}\right)\left(Y_{i}-m\right)^{2} \end{align}$$` - 本质上,以上就是一个 `\(Y\)`对截距项的**加权回归**(Weighted regression)估计量。 --- ### 2.2 滚动箱组回归:操作步骤 **滚动箱组线性回归**(Rolling regression)的操作步骤如下: - 根据计算点 `\(X=x_j\)`,按照特定谱宽 `\(h\)`,划分出若干**箱组**(bins)(箱组会有重叠): `\(\{b_1,b_2,\cdots,b_q\}\)`,其中 `\(b_{j}=[x_j-h,x_j +h]\)`。 - 根据样本数据集 `\(X_i\)`的实际情况,计算**矩形核函数**的权重值( `\(h=1\)`): `$$K\left(\frac{X_{i}-x_j}{h}\right)$$` - 然后采用**加权OLS**方法计算不同箱组下仅含截距项的回归模型 的截距系数 `\(\widehat{m}(x) \simeq m(X)\)` `$$Y=m(x)+e$$` --- ### (示例)滚动箱组回归:设定箱组划分规则 下面我们分别设定如下**箱组**划分规则: - 设定箱组取值中心点 `\((X=x_i),x_i \in (0.00, 0.01, 0.02, 0.03, 0.04, 0.05, \cdots, 9.95, 9.96, 9.97, 9.98, 9.99, 10.00)\)` - 确定核函数类型为**矩形核函数**<sup>a</sup>,谱宽设定为 `\(h=1\)`。 - 那么,可以得到箱组区块bins=[-1,1)、[-0.99,1.01)、[-0.98,1.02) `\(\ldots\)`[8.98,10.98)、[8.99,10.99)、[9,11),箱组数1001。 - 采用加权OLS方法进行拟合,直接计算得到 `\(\widehat{m}(x)\)` --- ### (示例)滚动箱组回归:数据计算表 利用箱组均值估计公式,我们可以计算得到不同箱组的均值估计:
- 基于设定的`x`表示 `\(x_i \in (0.00, 0.01, 0.02, 0.03, 0.04, 0.05,\cdots)\)`,共有 `\(N=1001\)`个 - 箱组回归估计`m1` `\(=\widehat{m}(x_i)\)` --- ### (示例)滚动箱组回归:图形表达 - 根据前面计算表的1001个拟合数据对 `\((x_i,\widehat{m}(x_i))\)`,我们可以得到滚动箱组线性回归估计结果: <img src="seminar-RDD-part01_files/figure-html/unnamed-chunk-26-1.svg" height="420" style="display: block; margin: auto;" /> --- name: ll-est ### 2.3 局部线性回归:渐进近似的另一种选择 > **回顾**:Nadaraya-Watson (NW)估计量的局部近似选择为 `\(m(X)\simeq m(x)\)` **局部线性回归**(Local linear regression, LLR)方法则选择了另一种局部线性近似,也即: `$$\begin{align} m(X)\simeq m(x) + m^{\prime}(x)(X-x) \\ \end{align}$$` 因此**局部线性**(LL)模型可以表达为: `$$\begin{align} Y &= m(X) +e \\ &\simeq m(x) + m^{\prime}(x)(X-x) +e \\ & = \beta_0 + \beta_1\cdot(X-x) +e \end{align}$$` - 以上模型可以视作为**一元线性回归模型**,其中 `\(\beta_0 =m(x);\beta_1=m^{\prime}(x)\)`。 - 我们的目标是估计得到 `\(\widehat{m}(x)\)`,也即上述模型的**截距项**。 --- ### 2.3 局部线性回归:矩阵表达式 前述一元线性模型也可以表达为如下矩阵形式: `$$\begin{align} Y &\simeq Z(X, x)^{\prime} \boldsymbol{\beta}(x)+e \end{align}$$` 其中: `$$\begin{align} Z(X, x) &=\left(\begin{array}{c} 1 \\ X-x \end{array}\right)\\ \boldsymbol{\beta}(x) &=\left(m(x), m^{\prime}(x)\right)^{\prime} \end{align}$$` --- ### 2.3 局部线性回归:最小化问题求解 对于前述渐近近似回归模型,最小化问题可以表达为: `$$\begin{align} \left\{\widehat{m}_{\mathrm{LL}}(x), \widehat{m}_{\mathrm{LL}}^{\prime}(x)\right\}= \underset{\beta_0, \beta_1}{\operatorname{argmin}} \sum_{i=1}^{n} K\left(\frac{X_{i}-x}{h}\right)\left(Y_{i}-\beta_0-\beta_{1}\left(X_{i}-x\right)\right)^{2} \end{align}$$` .fyi[ 我们可以发现两个有意思的结论: - 当谱宽趋近于无穷大时 `\(h \rightarrow \infty\)`,局部线性估计量将趋近于全样本OLS估计量 `\(\widehat{m}_{\mathrm{LL}}(x) \rightarrow \widehat{\beta}_0+\widehat{\beta}_1 x\)`。因为,此时所有样本的权重都会相等。 - 局部线性估计量会同时得到在点 `\(x\)`处,条件期望函数CEF的估计量 `\(\widehat{m}(x)\)`,及其斜率 `\(\widehat{m}^{\prime}(x)\)`。 ] --- ### 2.3 局部线性回归:估计量 以**核函数**为权重,以及利用**加权最小二乘法**原理,可以证明局部线性回归的估计量为: `$$\begin{aligned} \widehat{\beta}_{\mathrm{LL}}(x) &=\left(\sum_{i=1}^{n} K\left(\frac{X_{i}-x}{h}\right) Z_{i}(x) Z_{i}(x)^{\prime}\right)^{-1} \sum_{i=1}^{n} K\left(\frac{X_{i}-x}{h}\right) Z_{i}(x) Y_{i} \\ &=\left(\boldsymbol{Z}^{\prime} \boldsymbol{K} \boldsymbol{Z}\right)^{-1} \boldsymbol{Z}^{\prime} \boldsymbol{K} \boldsymbol{Y} \end{aligned}$$` 其中: - `\(\boldsymbol{K}=\operatorname{diag}\left\{K\left(\left(X_{1}-x\right) / h\right), \ldots, K\left(\left(X_{n}-x\right) / h\right)\right\}\)` - `\(\boldsymbol{Z}\)`是 `\({Z_i}(x)^{\prime}\)`的堆栈(stacked)形态 - `\(\boldsymbol{Y}\)`是 `\({Y_i}(x)^{\prime}\)`的堆栈形态 --- ### 2.3 局部线性回归:操作步骤 **局部线性回归**(Local Linear regression)的操作步骤如下: - 根据计算点 `\(X=x_j\)`,按照特定谱宽 `\(h\)`,划分出若干**箱组**(bins)(箱组会有重叠): `\(\{b_1,b_2,\cdots,b_q\}\)`,其中 `\(b_{j}=[x_j-h,x_j +h]\)`。 - 根据样本数据集 `\(X_i\)`的实际情况,计算**高斯核函数**的权重值( `\(h=1/\sqrt{3}\)`): `$$K\left(\frac{X_{i}-x_j}{h}\right)$$` - 然后采用**加权OLS**方法计算不同箱组下一元回归模型的截距系数 `\(\widehat{m}(x) \simeq m(X)\)` `$$\begin{align} Y &= m(X) +e \\ &\simeq m(x) + m^{\prime}(x)(X-x) +e \end{align}$$` --- ### (示例)局部线性回归:数据计算表 利用局部线性回归估计公式,我们可以计算得到不同箱组的估计:
- 基于设定的`x`表示 `\(x_i \in (0.00, 0.01, 0.02, 0.03, 0.04, 0.05,\cdots)\)`,共有 `\(N=1001\)`个 - 箱组回归估计`m2` `\(=\widehat{m}(x_i)\)` --- ### (示例)局部线性回归估计:图形表达 - 根据前面计算表的1001个拟合数据对 `\((x_i,\widehat{m}(x_i))\)`,我们可以得到局部线性回归估计结果: <img src="seminar-RDD-part01_files/figure-html/unnamed-chunk-28-1.svg" height="420" style="display: block; margin: auto;" /> --- ### (示例)CEF估计:三种回归估计方法的图形比较 <img src="seminar-RDD-part01_files/figure-html/unnamed-chunk-29-1.svg" height="450" style="display: block; margin: auto;" /> - 尽管三种回归估计方法都较好地估计了真实CEF的趋势,但是**局部线性回归**LLR拟合方法要更平滑。 --- name: lp-est ### 2.4 局部多项式回归:模型表达 **局部多项式回归**(Local Polynomial regression)模型可以表达为: `$$\begin{aligned} Y &=m(X)+e \\ & \simeq m(x)+m^{\prime}(x)(X-x)+\cdots+m^{(p)}(x) \frac{(X-x)^{p}}{p !}+e \\ &=Z(X, x)^{\prime} \beta(x)+e_{i} \end{aligned}$$` - 其中: `$$\begin{align} Z(X, x)=\left(\begin{array}{c} 1 \\ X-x \\ \vdots \\ \frac{(X-x)^{p}}{p !} \end{array}\right) \quad \quad \beta(x)=\left(\begin{array}{c} m(x) \\ m^{\prime}(x) \\ \vdots \\ m^{(p)}(x) \end{array}\right) \end{align}$$` --- ### 2.4 局部多项式回归:估计量 同样地,以**核函数**为权重,可以证明**局部多项式回归**估计量的理论公式为: `$$\begin{align} \widehat{\beta}_{\mathrm{LP}}(x) &=\left(\sum_{i=1}^{n} K\left(\frac{X_{i}-x}{h}\right) Z_{i}(x) Z_{i}(x)^{\prime}\right)^{-1}\left(\sum_{i=1}^{n} K\left(\frac{Y_{i}-x}{h}\right) Z_{i}(x) Y_{i}\right) \\ &=\left(\boldsymbol{Z}^{\prime} \boldsymbol{K} \boldsymbol{Z}\right)^{-1} \boldsymbol{Z}^{\prime} \boldsymbol{K} \boldsymbol{Y} \end{align}$$` - 其中 `\(Z_i(x)=Z(X_i,x)\)` --- ### 2.4 局部多项式回归:特点 .fyi[ - 局部多项式回归LPR具有一般性: > 其中有两种特定情形:a)当 `\(p=0\)`时为Nadaraya-Watson回归估计;b)当 `\(p=1\)`时为局部线性回归估计(LL) ] .fyi[ - 估计结果需要在阶数 `\(p\)`和局部平滑谱宽 `\(h\)`之间寻求平衡。 > 一方面如果增加多项式阶数 `\(p\)`,可以改进模型渐近精度,因而倾向选择更大的谱宽 `\(h\)`。另一方面,增加多项式阶数 `\(p\)`同时也会导致估计量方差的增大,估计可靠性会降低。 ] --- layout: false class: center, middle, duke-orange,hide_logo name: performance # 3.估计效果(Performance Analysis) ### [3.1 渐近偏误(Asymptotic Bias)](#asym-bias) ### [3.2 渐近方差(Asymptotic Variance)](#asym-var) ### [3.3 渐近均方误(AMSE & AIMSE)](#amse) ### [3.4 谱宽选择(Bandwidth Selection)](#band-sel) ### [3.5 渐近分布(Asymptotic Distribution)](#asym-dist) ### [3.6 方差估计(Variance Estimation)](#var-est) ### [3.7 工资案例](#case-wage) --- layout: true <div class="my-header-h2"></div> <div class="watermark1"></div> <div class="watermark2"></div> <div class="watermark3"></div> <div class="my-footer"><span><a href="https://home.huhuaping.com/">https://home.huhuaping.com </a>    <a href="#chapter-RDD01"> RDD Part01 非参数估计 </a>                     <a href="#performance"> 3.估计效果(Performance Analysis) </a> </span></div> --- name: asym-bias ### 3.1 渐近偏误:估计量的期望 期望: `$$\begin{align} \mathbb{E}\left[\widehat{m}_{\mathrm{nw}}(x) \mid \boldsymbol{X}\right]=\frac{\sum_{i=1}^{n} K\left(\frac{X_{i}-x}{h}\right) \mathbb{E}\left[Y_{i} \mid X_{i}\right]}{\sum_{i=1}^{n} K\left(\frac{X_{i}-x}{h}\right)}=\frac{\sum_{i=1}^{n} K\left(\frac{X_{i}-x}{h}\right) m\left(X_{i}\right)}{\sum_{i=1}^{n} K\left(\frac{X_{i}-x}{h}\right)} \end{align}$$` --- ### 3.1 渐近偏误:NW和LL下估计量的期望 对于局部NW估计量有: `$$\begin{align} &\mathbb{E}\left[\hat{m}_{\mathrm{nw}}(x) \mid \boldsymbol{X}\right]=m(x)+h^{2} B_{\mathrm{nw}}(x)+o_{p}\left(h^{2}\right)+O_{p}\left(\sqrt{\frac{h}{n}}\right) \\ &\qquad B_{\mathrm{nw}}(x)=\frac{1}{2} m^{\prime \prime}(x)+f(x)^{-1} f^{\prime}(x) m^{\prime}(x) \end{align}$$` - 其中: `\(\qquad B_{\mathrm{nw}}(x)=\frac{1}{2} m^{\prime \prime}(x)+f(x)^{-1} f^{\prime}(x) m^{\prime}(x)\)` 对于局部线性LL估计量有: `$$\begin{align} \mathbb{E}\left[\hat{m}_{\mathrm{LL}}(x) \mid \boldsymbol{X}\right]=m(x)+h^{2} B_{\mathrm{LL}}(x)+o_{p}\left(h^{2}\right)+O_{p}\left(\sqrt{\frac{h}{n}}\right) \end{align}$$` - 其中 `\(\qquad B_{\mathrm{LL}}(x)=\frac{1}{2} m^{\prime \prime}(x)\)` --- ### 3.1 渐近偏误:定义 **渐近偏误**(Asymptotic Bias):我们称 `\(h^2 B_{\mathrm{nw}}(x)\)`和 `\(h^2 B_{\mathrm{LL}}(x)\)`为**渐近偏误**。 --- exclude: true ## Code script: Hensen bruce fig 19-2 - 此代码仅用于生成本地图片,为提高渲染效率代码文件仅需运行一次。 - 如更新代码,请更改`eval = TRUE` --- ### (示例)不同谱宽下渐近偏误的表现 - 根据前面计算表的1001个拟合数据对 `\((x_i,\widehat{m}(x_i))\)`,我们可以得到滚动箱组线性回归估计结果: <img src="seminar-RDD-part01_files/figure-html/unnamed-chunk-31-1.svg" height="420" style="display: block; margin: auto;" /> --- name: asym-var ### 3.2 渐近方差: 渐近方差(Asymptotic Variance) `$$\begin{align} \hat{m}_{\mathrm{nw}}(x)-\mathbb{E}\left[\hat{m}_{\mathrm{nw}}(x) \mid \boldsymbol{X}\right]=\frac{\sum_{i=1}^{n} K\left(\frac{X_{i}-x}{h}\right) e_{i}}{\sum_{i=1}^{n} K\left(\frac{X_{i}-x}{h}\right)} \end{align}$$` `$$\begin{align} \operatorname{var}\left[\hat{m}_{\mathrm{nw}}(x) \mid \boldsymbol{X}\right]=\frac{\sum_{i=1}^{n} K\left(\frac{X_{i}-x}{h}\right)^{2} \sigma^{2}\left(X_{i}\right)}{\left(\sum_{i=1}^{n} K\left(\frac{X_{i}-x}{h}\right)\right)^{2}} \end{align}$$` --- ### 3.2 渐近方差: `$$\begin{align} &\operatorname{var}\left[\hat{m}_{\mathrm{nw}}(x) \mid \boldsymbol{X}\right]=\frac{R_{K} \sigma^{2}(x)}{f(x) n h}+o_{p}\left(\frac{1}{n h}\right) \\ &\operatorname{var}\left[\hat{m}_{\mathrm{LL}}(x) \mid \boldsymbol{X}\right]=\frac{R_{K} \sigma^{2}(x)}{f(x) n h}+o_{p}\left(\frac{1}{n h}\right) \end{align}$$` - **核函数** `\(K(u)\)`的**粗糙度**(roughness)定义为: `\(R_{K}=\int_{-\infty}^{\infty} K(u)^{2} d u\)` --- name: amse ### 3.3 渐近均方误:定义 **渐近均方误**(Asymptotic MSE, AMSE):是估计量 `\(\widehat{m}(x)\)`的**平方渐近偏误**以及**渐近误差**二者之和。定义为: `$$\begin{align} \operatorname{AMSE}(x) \stackrel{\text { def }}{=} h^{4} B(x)^{2}+\frac{R_{K} \sigma^{2}(x)}{n h f(x)} \end{align}$$` --- ### 3.3 渐近均方误:连续情形下 **渐近积分均方误**(Asymptotic integrated MSE, AIMSE):类似地定义如下: `$$\begin{align} \operatorname{AIMSE} &\stackrel{\text { def }}{=} \int_{S} \operatorname{AMSE}(x) f(x) w(x) d x \\ &=\int_{S}\left(h^{4} B(x)^{2}+\frac{R_{K} \sigma^{2}(x)}{n h f(x)}\right) f(x) w(x) d x \\ &=h^{4} \bar{B}+\frac{R_{K}}{n h} \bar{\sigma}^{2} \end{align}$$` - 其中: `$$\begin{aligned} \bar{B} &=\int_{S} B(x)^{2} f(x) w(x) d x \\ \bar{\sigma}^{2} &=\int_{S} \sigma^{2}(x) w(x) d x \end{aligned}$$` --- name: band-sel ### 3.4 谱宽选择:最优谱宽 最小化**渐近积分均方误**目标下,可以得到**最优谱宽**(Optimal Bandwidth): `$$\begin{align} h_{0}=\left(\frac{R_{K} \bar{\sigma}^{2}}{4 \bar{B}}\right)^{1 / 5} n^{-1 / 5} \end{align}$$` - 此时,随着 `\(h \sim n^{-1 / 5}\)`,则有 `\(\text {AIMSE }[\hat{m}(x)]=O\left(n^{-4 / 5}\right)\)` - 因此,可以计算出上述最优谱宽下的**渐近积分均方误**: `$$\begin{align} \operatorname{AIMSE}_{0} \simeq 1.65\left(R_{K}^{4} \bar{B} \bar{\sigma}^{8}\right)^{1 / 5} n^{-4 / 5} \end{align}$$` --- ### 3.4 谱宽选择:最优谱宽 至此,我们可以得到最优谱宽的理论计算公式: `$$\begin{align} h_{0}=\left(\frac{R_{K}}{4}\right)^{1 / 5}\left(\frac{\bar{\sigma}^{2}}{n \bar{B}}\right)^{1 / 5} \simeq 0.58\left(\frac{\bar{\sigma}^{2}}{n \bar{B}}\right)^{1 / 5} \end{align}$$` `$$\begin{align} \bar{B}=\mathbb{E}\left[B(X)^{2} w(X)\right]=\mathbb{E}\left[\left(\frac{1}{2} m^{\prime \prime}(X)\right)^{2} \mathbb{1}\left\{\xi_{1} \leq X \leq \xi_{2}\right\}\right] . \end{align}$$` --- ### 3.4 谱宽选择:参考谱宽Rot Fan and Gijbels(1996)提出了一种**经验参考谱宽**(Rule of Thumb bandwidth,ROT)的计算办法: `$$\begin{align} h_{\text {rot }}=0.58\left(\frac{\widehat{\sigma}^{2}\left(\xi_{2}-\xi_{1}\right)}{n \widehat{B}}\right)^{1 / 5} \end{align}$$` - 首先构建**先验q阶多项式回归模型**,并分别估计得到 `\(\widehat{m}(x)\)`及其二阶导 `\(\widehat{m}^{\prime \prime}(x)\)` `$$\begin{align} {m}(x)&={\beta}_{0}+{\beta}_{1} x+{\beta}_{2} x^{2}+\cdots+{\beta}_{q} x^{q} + {\epsilon} \\ \end{align}$$` - 并使用 `\(\overline{B}\)`的**矩估计量**: `\(\widehat{B}=\frac{1}{n} \sum_{i=1}^{n}\left(\frac{1}{2} \hat{m}^{\prime \prime}\left(X_{i}\right)\right)^{2} \mathbb{1}\left\{\xi_{1} \leq X_{i} \leq \xi_{2}\right\}\)` - 其次,假定随机干扰项为同方差,也即 `\(\mathbb{E}\left[e^{2} \mid X\right]=\sigma^{2}\)`,从而可以使用: `$$\begin{align} \bar{\sigma}^{2}=\sigma^{2}\left(\xi_{2}-\xi_{1}\right) \approx \hat{\sigma}^{2}\left(\xi_{2}-\xi_{1}\right) \end{align}$$` --- ### 3.4 谱宽选择:参考谱宽Rot(主要步骤) 参考谱宽rot的主要机选步骤具体包括: - 步骤1:确定权重取值范围 `\(w(x)=\mathbb{1}\left\{\xi_{1} \leq x \leq \xi_{2}\right\}\)` - 步骤2:构建**先验q阶多项式回归模型**(建议为4阶),并分别估计得到 `\(\widehat{m}(x)\)`及其二阶导 `\(\widehat{m}^{\prime \prime}(x)\)` `$$\begin{align} {m}(x)&={\beta}_{0}+{\beta}_{1} x+{\beta}_{2} x^{2}+\cdots+{\beta}_{q} x^{q} + {\epsilon} \\ \widehat{m}(x)&=\hat{\beta}_{0}+\hat{\beta}_{1} x+\hat{\beta}_{2} x^{2}+\cdots+\hat{\beta}_{q} x^{q} \\ \widehat{m}^{\prime \prime}(x) &=2\hat{\beta}_{2}+6\hat{\beta}_{3} x + 12\hat{\beta}_{4} x^{2}+\cdots+q(q-1)\hat{\beta}_{q}x^{q-2} \end{align}$$` - 步骤3:利用上述估计结果计算 `$$\begin{align} \widehat{B}=\frac{1}{n} \sum_{i=1}^{n}\left(\frac{1}{2} \hat{m}^{\prime \prime}\left(X_{i}\right)\right)^{2} \mathbb{1}\left\{\xi_{1} \leq X_{i} \leq \xi_{2}\right\} \end{align}$$` --- ### 3.4 谱宽选择:参考谱宽Rot(主要步骤) 参考谱宽rot的主要机选步骤具体包括: - 步骤4:计算上述**先验q阶多项式回归模型**的回归误差方差 `\(\hat{\sigma}^2\)` `$$\hat{\sigma}^2 = \frac{\sum{\hat{\epsilon}^2}}{n-q-1}$$` - 步骤5:根据上述全部结果计算得到经验谱宽: `$$\begin{align} h_{\text {rot }}=0.58\left(\frac{\widehat{\sigma}^{2}\left(\xi_{2}-\xi_{1}\right)}{n \widehat{B}}\right)^{1 / 5} \end{align}$$` --- exclude: true ## Code script: Hensen bruce fig 19-4 - 此代码仅用于生成本地图片,为提高渲染效率代码文件仅需运行一次。 - 如更新代码,请更改`eval = TRUE` --- ### (示例)参考谱宽的计算:数据散点图 - 我们继续使用前述的蒙特卡洛模拟数据集: <img src="seminar-RDD-part01_files/figure-html/unnamed-chunk-33-1.svg" height="450" style="display: block; margin: auto;" /> --- ### (示例)参考谱宽的计算:多项式回归 下面我们按前述步骤来计算参考谱宽值 `\(h_{rot}\)`: - **步骤1**:根据案例数据集,设定权重取值范围 `\(\left\{\xi_{1} \leq x \leq \xi_{2}\right\}=\{0,10\}\)` - **步骤2**:构建多项式回归 `$$\begin{align} \begin{split} Y_i=&+\beta_{1}+\beta_{2}X_i+\beta_{3}X^2_i+\beta_{4}X^3_i+\beta_{5}X^4_i+u_i \end{split} \end{align}$$` > 直接使用OLS进行估计,得到估计方程: `$$\begin{equation} \begin{alignedat}{999} &\widehat{Y}=&&+0.4943&&+0.6986X_i&&-0.2807X^2_i&&+0.0326X^3_i&&-0.0012X^4_i\\ &(s)&&(0.1521)&&(0.2015)&&(0.0802)&&(0.0119)&&(0.0006) \end{alignedat} \end{equation}$$` > 进而得到拟合值 `\(\widehat{m}(x)\)`及其二阶导 `\(\widehat{m}^{\prime \prime}(x)\)` `$$\begin{align} \widehat{m}(x)&=+0.4943+0.6986x_i-0.2807x^2_i +0.0326x^3_i-0.0012x^4_i \\ \widehat{m}^{\prime \prime}(x)&=-2\times 0.2807+6\times0.0326x_i-12\times0.0012x^2_i \\ \end{align}$$` --- ### (示例)参考谱宽的计算:多项式回归 - 拟合值 `\(\widehat{m}(x)\)`及其二阶导 `\(\widehat{m}^{\prime \prime}(x)\)`以及残差 `\(\hat{\epsilon}\)`
--- ### (示例)参考谱宽的计算:结果 - 步骤3:利用上述估计结果计算 `$$\begin{align} \widehat{B}=\frac{1}{n} \sum_{i=1}^{n}\left(\frac{1}{2} \hat{m}^{\prime \prime}\left(X_{i}\right)\right)^{2} \mathbb{1}\left\{\xi_{1} \leq X_{i} \leq \xi_{2}\right\} = 0.0089 \end{align}$$` - 步骤4:多项式模型的回归误差方差 `$$\hat{\sigma}^2 = \frac{\sum{\hat{\epsilon}^2}}{n-q-1}=0.0687$$` - 步骤5:根据上述全部结果计算得到经验谱宽: `$$\begin{align} h_{\text {rot }}&=0.58\left(\frac{\widehat{\sigma}^{2}\left(\xi_{2}-\xi_{1}\right)}{n \widehat{B}}\right)^{1 / 5} = 0.58 \times \left(\frac{0.0687\times \left(10-0\right)}{100 \times 0.0089}\right)^{1 / 5} = 0.5508 \end{align}$$` --- ### 3.4 谱宽选择:交叉验证谱宽(目标问题) 我们期望选择谱宽,以实现估计量 `\(\widehat{m}(x, h)\)`最小化**积分均方误**(Integrated mean-squared error, IMSE),也即: `$$\begin{align} \operatorname{IMSE}_{n}(h)=\int_{S} \mathbb{E}\left[(\widehat{m}(x, h)-m(x))^{2}\right] f(x) w(x) d x \end{align}$$` - 其中 `\(f(x)\)`为 `\(X\)`的边际密度(marginal density) - `\(w(x)\)`为可积分权重函数(integrable weight function) --- ### 3.4 谱宽选择:交叉验证谱宽(可行计算) 上述最小化问题中,偏差值 `\((\widehat{m}(x, h)-m(x))\)`可以通过**留一法**下的预测误差来代替,也即: `$$\begin{align} \tilde{e}_{i}(h)=Y_{i}-\tilde{m}_{-i}\left(X_{i}, h\right) \end{align}$$` 因此,上述最小化问题 `\(IMSE_n(h)\)`的一个**可行计算方案**可以表达为如下的**交叉验证准则函数**(Cross-Validation Criterion): `$$\begin{align} \operatorname{CV}(h)=\frac{1}{n} \sum_{i=1}^{n} \widetilde{e}_{i}(h)^{2} w\left(X_{i}\right) \end{align}$$` --- ### 3.4 谱宽选择:交叉验证谱宽(准则函数) 对于**交叉验证准则函数**(Cross-Validation Criterion): `$$\begin{align} \operatorname{CV}(h)=\frac{1}{n} \sum_{i=1}^{n} \widetilde{e}_{i}(h)^{2} w\left(X_{i}\right) \end{align}$$` 我们可以证明它是去掉一个样本的积分均方误 `\(IMSE_(n-1)(h)\)`加上一个常数项的**无偏估计**,也即: `$$\mathbb{E}[\mathrm{CV}(h)]=\bar{\sigma}^{2}+\operatorname{IMSE}_{n-1}(h)$$` - 其中: `\(\bar{\sigma}^{2}=\mathbb{E}\left[e^{2} w(X)\right]\)`,它是不依赖于谱宽 `\(h\)`的。 - 显然,最小化 `\(\mathbb{E}[\mathrm{CV}(h)]\)`和最小化 `\(\operatorname{IMSE}_{n-1}(h)\)`将是等价的。 - 而且,当 `\(n\)`比较大时,最小化 `\(\operatorname{IMSE}_{n-1}(h)\)`和最小化 `\(\operatorname{IMSE}_{n}(h)\)`也将是等价的,二者各自得到的谱宽也是无偏的。 --- ### 3.4 谱宽选择:交叉验证谱宽(受约束) 另外,为了避免谱宽值选择过小,需要给定约束条件 `\(h \geq h_{\ell}\)`,此时: `$$h_{\mathrm{CV}}=\underset{h \geq h_{\ell}}{\operatorname{argmin}} \mathrm{CV}(h)$$` --- ### 3.4 谱宽选择:交叉验证谱宽(主要过程) 因此,最优交叉验证谱宽选择的主要过程如下: - 构建 `\(h\)`的序贯数值(grid value) `\([h_1, h_2, \cdots, h_J]\)`,并评估交叉效用准则函数的最小化取值 `\(CV(h_j),j\in(1,2,\cdots,J)\)`,从而得到最优谱宽: `$$h_{\mathrm{CV}}=\underset{h \in\left[h_{1}, h_{2}, \ldots, h_{J}\right]}{\operatorname{argmin}} \operatorname{CV}(h)$$` .fyi[ - 需要注意的是,以上方法得到的谱宽理论上可以是无界化的。意味着,交叉验证准则函数 `\(CV(h)\)`是单调下降的,以至于最优谱宽 `\(h = \infty\)`。 - 此时,将表明**全样本**回归估计也可能是一个最优结果。也即,使用全样本数据,NW估计方法有 `\(\widehat{m}(x)=\overline{Y}\)`;而局部线性估计方法则有 `\(\widehat{m}(x)=\hat{\beta}_0 +\hat{\beta}_1x\)` ] --- ### 3.4 谱宽选择:交叉验证谱宽(计算步骤) - 步骤1:设定**初始值**。也即选定一个**参考谱宽**作为初始值,例如前面提到的**经验谱宽** `\(h_{rot}\)`。 - 步骤2:设定**调参谱宽**(tuning bandwidth)。也即给定待评估的谱宽范围,及其待评估序贯值(grid value)。 > - 一个经验谱宽范围可供参考: `\([h_{rot}/3, 3h_{rot}]\)`。 - 范围内的待评估序贯值的个数 `\(g\)`也是需要做出尝试性选择。 --- ### (示例)交叉验证谱宽的计算:数据集 - 我们继续使用前述的蒙特卡洛模拟数据集: <img src="seminar-RDD-part01_files/figure-html/unnamed-chunk-37-1.svg" height="450" style="display: block; margin: auto;" /> --- ### (示例)交叉验证谱宽的计算:规则 - 步骤1:设定**经验谱宽** `\(h_{rot}=0.5508\)`作为**初始值**。。 - 步骤2:设定**调参谱宽**(tuning bandwidth)。 > - 一个经验谱宽范围可供参考: `\([h_{rot}/3, 3h_{rot}]=[0.1836, 1.6524]\)`。 - 给定范围内的搜寻总数为 `\(n=201\)`。则待评估序贯值为 `\(h\in (0.1836, 0.1909, 0.1983, 0.2056, 0.2130, \cdots,1.6304, 1.6377, 1.6451, 1.6524)\)`。 --- exclude: true ### R chunk: Helper function --- ### (示例)交叉验证谱宽的计算:NW方法下的预测误差平方
.footnote[ <sup>a</sup> `loo`表示交叉验证**留一法**,例如`loo_001`表示第1个样本数据不进入局部回归估计。 ] --- ### (示例)交叉验证谱宽的计算:LL方法下的预测误差平方
.footnote[ <sup>a</sup> `loo`表示交叉验证**留一法**,例如`loo_001`表示第1个样本数据不进入局部回归估计。 ] --- ### (示例)交叉验证谱宽的计算:CV计算表
--- ### (示例)交叉验证谱宽的计算:谱宽与CV变化 <img src="seminar-RDD-part01_files/figure-html/unnamed-chunk-42-1.svg" height="450" style="display: block; margin: auto;" /> --- ### (示例)最优谱宽下的估计表
--- ### (示例)真实的条件期望函数m(x) <img src="seminar-RDD-part01_files/figure-html/unnamed-chunk-44-1.svg" height="420" style="display: block; margin: auto;" /> `$$\begin{align} m(X) =\frac{sin(\frac{\pi}{4}\cdot(X_i-2))}{\frac{\pi}{4}\cdot(X_i-2)} \end{align}$$` --- ### (示例)LL方法下使用ROT谱宽估计得到的m(x) <img src="seminar-RDD-part01_files/figure-html/unnamed-chunk-45-1.svg" height="450" style="display: block; margin: auto;" /> --- ### (示例)LL方法下使用最优CV谱宽估计得到的m(x) <img src="seminar-RDD-part01_files/figure-html/unnamed-chunk-46-1.svg" height="450" style="display: block; margin: auto;" /> --- ### (示例)LL方法下使用不同谱宽估计得到的m(x):对比 <img src="seminar-RDD-part01_files/figure-html/unnamed-chunk-47-1.svg" height="450" style="display: block; margin: auto;" /> --- ### (示例)模拟数据集 为了更好地进行数据验证,我们将根据如下规则生成蒙特卡洛模拟数据集: `$$\begin{align} m(X) =\frac{sin(\frac{\pi}{4}\cdot(X_i-2))}{\frac{\pi}{4}\cdot(X_i-2)} \end{align}$$` - 此时,我们具有**上帝视角**,实际上已经知道**数据生成机制**(DGP) - 此时,我们心里面已知真实模型为**非线性的** --- name: asym-dist ### 3.5 渐近分布:渐近正态分布 渐近分布(Asymptotic Distribution) - 对于局部NW估计: `$$\begin{align} \sqrt{n h}\left(\widehat{m}_{\mathrm{nw}}(x)-m(x)-h^{2} B_{\mathrm{nw}}(x)\right) \underset{d}{\longrightarrow} \mathrm{N}\left(0, \frac{R_{K} \sigma^{2}(x)}{f(x)}\right) \end{align}$$` - 对于局部线性估计: `$$\begin{align} \sqrt{n h}\left(\widehat{m}_{\mathrm{LL}}(x)-m(x)-h^{2} B_{\mathrm{LL}}(x)\right) \underset{d}{\longrightarrow} \mathrm{N}\left(0, \frac{R_{K} \sigma^{2}(x)}{f(x)}\right) \end{align}$$` --- ### 3.5 渐近分布:超平滑情形 `$$\begin{align} \sqrt{n h}\left(\hat{m}_{\mathrm{nw}}(x)-m(x)\right) &\underset{d}{\longrightarrow} \mathrm{N}\left(0, \frac{R_{K} \sigma^{2}(x)}{f(x)}\right) \\ \sqrt{n h}\left(\widehat{m}_{\mathrm{LL}}(x)-m(x)\right) &\underset{d}{\longrightarrow} \mathrm{N}\left(0, \frac{R_{K} \sigma^{2}(x)}{f(x)}\right) \end{align}$$` --- name: var-est ### 3.6 方差估计:条件方差 `$$\begin{align} \sigma^{2}(x)=\operatorname{var}[Y \mid X=x]=\mathbb{E}\left[e^{2} \mid X=x\right] \end{align}$$` - NW方法下条件方差的一个**理想**估计为: `$$\begin{align} \bar{\sigma}^{2}(x)=\frac{\sum_{i=1}^{n} K\left(\frac{X_{i}-x}{h}\right) e_{i}^{2}}{\sum_{i=1}^{n} K\left(\frac{X_{i}-x}{h}\right)} \end{align}$$` - NW方法下条件方差的一个**可行**估计为: `$$\begin{align} \widehat{\sigma}^{2}(x)=\frac{\sum_{i=1}^{n} K\left(\frac{X_{i}-x}{h}\right) \widetilde{e}_{i}^{2}}{\sum_{i=1}^{n} K\left(\frac{X_{i}-x}{h}\right)} . \end{align}$$` > 其中 `\(\widetilde{e}_{i}=Y_{i}-\hat{m}_{-i}\left(X_{i}\right)\)`为**留一法**下的预测误差(leave-one-out prediction error) --- ### 3.6 方差估计:直接公式 如前所属,NW方法、LL方法和LP方法的CEF估计可以统一表达为: `$$\begin{align} \widehat{\beta}(x) &=\left(\boldsymbol{Z}^{\prime} \boldsymbol{K} \boldsymbol{Z}\right)^{-1}\left(\boldsymbol{Z}^{\prime} \boldsymbol{K} \boldsymbol{Y}\right) \\ &=\left(\boldsymbol{Z}^{\prime} \boldsymbol{K} \boldsymbol{Z}\right)^{-1}\left(\boldsymbol{Z}^{\prime} \boldsymbol{K} \boldsymbol{m}\right)+\left(\boldsymbol{Z}^{\prime} \boldsymbol{K} \boldsymbol{Z}\right)^{-1}\left(\boldsymbol{Z}^{\prime} \boldsymbol{K} \boldsymbol{e}\right) \end{align}$$` - 那么,我们可以直接使用下列条件方差公式: `$$\begin{align} \boldsymbol{V}_{\widehat{\beta}}(x)=\operatorname{var}[\widehat{\beta} \mid \boldsymbol{X}]=\left(\boldsymbol{Z}^{\prime} \boldsymbol{K} \boldsymbol{Z}\right)^{-1}\left(\boldsymbol{Z}^{\prime} \boldsymbol{K} \boldsymbol{D} \boldsymbol{K} \boldsymbol{Z}\right)\left(\boldsymbol{Z}^{\prime} \boldsymbol{K} \boldsymbol{Z}\right)^{-1} \end{align}$$` - 其中,我们可以直接使用**平方误差** `\(\widehat{e}_i^2\)`或者**平方预测误差** `\(\tilde{e}_i^2\)`: `$$\begin{align} \widehat{\boldsymbol{V}}_{\widehat{\beta}}(x)=\left(\boldsymbol{Z}^{\prime} \boldsymbol{K} \boldsymbol{Z}\right)^{-1}\left(\sum_{i=1}^{n} K\left(\frac{X_{i}-x}{h}\right)^{2} Z_{i}(x) Z_{i}(x)^{\prime} \widetilde{\boldsymbol{e}}_{i}^{2}\right)\left(\boldsymbol{Z}^{\prime} \boldsymbol{K} \boldsymbol{Z}\right)^{-1} \end{align}$$` --- ### 3.6 方差估计:渐近公式 更为简洁地,我们也可以使用如下条件方差的渐近公式: `$$\begin{align} \widehat{V}_{\widehat{m}(x)}=\frac{R_{K} \widehat{\sigma}^{2}(x)}{n h \widehat{f}(x)} \end{align}$$` - 其中: `$$\begin{align} \widehat{f}(x)=\frac{1}{n b} \sum_{i=1}^{n} K\left(\frac{X_{i}-x}{b}\right) \end{align}$$` `$$\begin{align} \widehat{\sigma}^{2}(x)=\frac{\sum_{i=1}^{n} K\left(\frac{X_{i}-x}{h}\right) \widetilde{e}_{i}^{2}}{\sum_{i=1}^{n} K\left(\frac{X_{i}-x}{h}\right)} . \end{align}$$` --- ### 3.6 方差估计:置信带 **逐点置信区间**(Pointwise Confidence Interval) `$$\begin{align} \widehat{m}(x) \pm z_{1-\alpha/2}(n-1) \cdot \sqrt{\widehat{V}_{\widehat{m}(x)}}\\ \widehat{m}(x) \pm 1.96 \sqrt{\widehat{V}_{\widehat{m}(x)}} \end{align}$$` --- exclude: true ## Code script: Hensen bruce fig 19-5 - 此代码运行时间较长(约3-5分钟),为提高渲染效率代码文件仅需运行一次。 - 如更新代码,请更改`eval = TRUE` --- name: case-wage ### 3.7 工资案例:背景说明 .case[ **工资案例**: - 案例基于CPS数据集,重点分析其中的**子样本数据**(黑人、男性、拥有12年受教育程度——高中毕业),样本数为n=762。 - 关注的问题:**时均工资**的对数(`Y =log(wage)`)对**职业经历**(`X=experience`)的非参数回归估计。 - 后面的分析中,我们会重点划定观测窗口为:职业经历(年数)范围为 `\([0,40]\)`。因为样本中90%以上的观测对象都落在这个范围之内。 ] --- ### (工资案例)样本数据集 .pull-left[
] .pull-right[ - 样本数据的描述性统计如下: ``` X.age Y.earnings Min. :-2 Min. :-3.1 1st Qu.:15 1st Qu.: 2.2 Median :25 Median : 2.5 Mean :25 Mean : 2.5 3rd Qu.:34 3rd Qu.: 2.8 Max. :62 Max. : 4.0 ``` ] --- ### (工资案例)样本数据散点图 <img src="seminar-RDD-part01_files/figure-html/unnamed-chunk-52-1.svg" height="450" style="display: block; margin: auto;" /> --- ### (工资案例)参考谱宽的计算:多项式回归 下面我们按前述步骤来计算参考谱宽值 `\(h_{rot}\)`: - **步骤1**:根据案例数据集,设定权重取值范围 `\(\left\{\xi_{1} \leq x \leq \xi_{2}\right\}=\{0,40\}\)` - **步骤2**:构建多项式回归 `$$\begin{align} \begin{split} Y_i=&+\beta_{1}+\beta_{2}X_i+\beta_{3}X^2_i+\beta_{4}X^3_i+\beta_{5}X^4_i+u_i \end{split} \end{align}$$` > 直接使用OLS进行估计,得到估计方程: `$$\begin{equation} \begin{alignedat}{999} &\widehat{Y}=&&+2.094395&&+0.030951X_i&&-0.000103X^2_i&&-0.000021X^3_i&&+0.000000X^4_i\\ &(s)&&(0.1371)&&(0.0284)&&(0.0019)&&(0.0000)&&(0.0000) \end{alignedat} \end{equation}$$` > 进而得到拟合值 `\(\widehat{m}(x)\)`及其二阶导 `\(\widehat{m}^{\prime \prime}(x)\)`及残差 `\(\hat{\epsilon}\)` `$$\begin{align} \widehat{m}(x) &=2.094395+0.030951x_i-0.000103x^2_i-0.000021x^3_i+0.000000x^4_i \\ \widehat{m}^{\prime \prime}(x) &=-2\times0.000103-6\times0.000021x_i+12\times0.000000x^2_i \\ \end{align}$$` --- ### (工资案例)参考谱宽的计算:结果 - 步骤3:利用上述估计结果计算 `$$\begin{align} \widehat{B}=\frac{1}{n} \sum_{i=1}^{n}\left(\frac{1}{2} \hat{m}^{\prime \prime}\left(X_{i}\right)\right)^{2} \mathbb{1}\left\{\xi_{1} \leq X_{i} \leq \xi_{2}\right\} = 0.00000025 \end{align}$$` - 步骤4:多项式模型的回归误差方差 `$$\hat{\sigma}^2 = \frac{\sum{\hat{\epsilon}^2}}{n-q-1}=0.2592$$` - 步骤5:根据上述全部结果计算得到经验谱宽: `$$\begin{align} h_{\text {rot }}&=0.58\left(\frac{\widehat{\sigma}^{2}\left(\xi_{2}-\xi_{1}\right)}{n \widehat{B}}\right)^{1 / 5} = 0.58 \times \left(\frac{0.2592\times \left(40-0\right)}{762 \times 0.000000248}\right)^{1 / 5} = 5.1442 \end{align}$$` --- ### (工资案例)交叉验证谱宽的计算:规则 - 步骤1:设定**经验谱宽** `\(h_{rot}=5.1442\)`作为**初始值**。。 - 步骤2:设定**调参谱宽**(tuning bandwidth)。 > - 一个经验谱宽范围可供参考: `\([h_{rot}/3, 3h_{rot}]=[1.7147, 15.4326]\)`。 - 给定范围内的搜寻总数为 `\(n=201\)`。则待评估序贯值为 `\(h\in (1.7147, 1.7833, 1.8519, 1.9205, 1.9891, \cdots,15.2268, 15.2954, 15.3640, 15.4326)\)`。 --- ### (工资案例)交叉验证谱宽的计算:CV计算表
--- ### (工资案例)交叉验证谱宽的计算:谱宽与CV变化 <img src="seminar-RDD-part01_files/figure-html/unnamed-chunk-56-1.svg" height="450" style="display: block; margin: auto;" /> --- ### (工资案例)最优谱宽选择下的估计表
--- ### (工资案例)LL方法下使用ROT谱宽估计得到的m(x) <img src="seminar-RDD-part01_files/figure-html/unnamed-chunk-58-1.svg" height="450" style="display: block; margin: auto;" /> --- ### (工资案例)LL方法下使用最优CV谱宽估计得到的m(x) <img src="seminar-RDD-part01_files/figure-html/unnamed-chunk-59-1.svg" height="450" style="display: block; margin: auto;" /> --- ### (工资案例)LL方法下使用不同谱宽估计得到的m(x):对比 <img src="seminar-RDD-part01_files/figure-html/unnamed-chunk-60-1.svg" height="450" style="display: block; margin: auto;" /> --- exclude: true ## Code script: Hensen bruce fig 19-6 - 此代码运行时间较长(约3-5分钟),为提高渲染效率代码文件仅需运行一次。 - 如更新代码,请更改`eval = TRUE` --- ### (工资案例)方差估计:获得预测残差 (1)我们首先可以计算得到LL估计下的预测残差的平方 `\(\tilde{e}^2_i\)` - 这一步可以直接采用前述的参考谱宽 `\(h_{rot}=5.1442\)` --- ### (工资案例)方差估计:再次获得参考谱宽 (2)然后开始计算方差估计下的最优谱宽。这里我们再进行一次参考谱宽的计算流程。 - 构建残差平方的多项式回归模型 `$$\begin{align} \begin{split} \tilde{e}^2_i=&+\gamma_{0}+\gamma_{1}X_i+\gamma_{2}X^2_i+\gamma_{3}X^3_i+\gamma_{4}X^4_i+v_i \end{split} \end{align}$$` - 利用ROT公式流程,再次获得参考谱宽 `$$\begin{align} hv_{\text {rot }}&=0.58\left(\frac{\widehat{\sigma}^{2}\left(\xi_{2}-\xi_{1}\right)}{n \widehat{B}}\right)^{1 / 5} = 0.58 \times \left(\frac{1.2384\times \left(40-0\right)}{762 \times 0.000000300}\right)^{1 / 5} = 6.7708 \end{align}$$` --- ### (工资案例)方差估计:再次获得最优交叉验证谱宽 (3)然后开始计算方差估计下的最优谱宽。这里我们再进行一次参考谱宽的计算流程。 - 步骤1:设定**经验谱宽** `\(hv_{rot}=6.7708\)`作为**初始值**。。 - 步骤2:设定**调参谱宽**(tuning bandwidth)。 > - 一个经验谱宽范围可供参考: `\([2.0, 40.0]\)`。 - 给定范围内的搜寻总数为 `\(n=202\)`。则待评估序贯值为 `\(h\in (2.0000, 2.1891, 2.3781, 2.5672, 2.7562, \cdots,39.4328, 39.6219, 39.8109, 40.0000)\)`。 - 步骤3:采用交叉验证**留一法**,分别遍历计算**NW估计**和**LL估计**下的全部CV值(见后面计算表) - 步骤4:最小CV值对应的谱宽评估值,则为最优交叉验证谱宽。当然,我们最终发现**NW估计**和**LL估计**下的结果是一样的,都选择了最大边界值 `\(hv_{CV}(NW)=hv_{CV}(LL)=40\)`。(见后面的CV比较图) --- ### (工资案例)方差估计:CV值计算表(附表)
--- ### (工资案例)方差估计:CV值与谱宽(附图) <img src="seminar-RDD-part01_files/figure-html/unnamed-chunk-64-1.svg" height="420" style="display: block; margin: auto;" /> - LL估计的CV值在局部上具有最小值,也即局部最优谱宽约为 `\(hv_{cv}(LL)\simeq 5\)` - 但是从全局来看,无论是NW估计,还是LL估计,CV函数值都表现为下降趋势。因此它们都选择了最大边界值 `\(hv_{CV}(NW)=hv_{CV}(LL)=40\)`。 --- ### (工资案例)方差估计:计算方差、标准差 (4)利用前面的平方预测误差,并使用谱宽 `\(h=5.1442\)`进行LL估计,最终得到方差和标准差估计值(见后面附表)。 `$$\begin{align} \widehat{\boldsymbol{V}}_{\widehat{\beta}}(x)=\left(\boldsymbol{Z}^{\prime} \boldsymbol{K} \boldsymbol{Z}\right)^{-1}\left(\sum_{i=1}^{n} K\left(\frac{X_{i}-x}{h}\right)^{2} Z_{i}(x) Z_{i}(x)^{\prime} \widetilde{\boldsymbol{e}}_{i}^{2}\right)\left(\boldsymbol{Z}^{\prime} \boldsymbol{K} \boldsymbol{Z}\right)^{-1} \end{align}$$` --- ### (工资案例)方差估计:计算方差估计值(附表)
--- ### (工资案例)置信区间和置信带 (5)进一步计算**逐点置信区间**(Pointwise Confidence Interval)(见前面附表),并得到置信带(见后面附图)。 `$$\begin{align} \widehat{m}(x) \pm z_{1-\alpha/2}(n-1) \cdot \sqrt{\widehat{V}_{\widehat{m}(x)}}\\ \widehat{m}(x) \pm 1.96 \sqrt{\widehat{V}_{\widehat{m}(x)}} \end{align}$$` --- ### (工资案例)方差估计:置信区间和置信带(附图) <img src="seminar-RDD-part01_files/figure-html/unnamed-chunk-66-1.svg" height="450" style="display: block; margin: auto;" /> --- layout: false class: center, middle, duke-orange,hide_logo name: cluster # 4.群组分析(Cluster observations) ### [4.1 基本原理](#basic-principal) ### [4.2 成绩案例](#case-score) --- layout: true <div class="my-header-h2"></div> <div class="watermark1"></div> <div class="watermark2"></div> <div class="watermark3"></div> <div class="my-footer"><span><a href="https://home.huhuaping.com/">https://home.huhuaping.com </a>    <a href="#chapter-RDD01"> RDD Part01 非参数估计 </a>                     <a href="#cluster"> 4.群组分析(Cluster observations) </a> </span></div> --- name: basic-principal ### 4.1 基本原理:模型表达 **群组观测**(Clustered observations)定义为:观测个体 `\(i =1,2,\cdots,n_g\)`分别处在群组 `\(g = 1,2, \cdots, G\)`内,并可观测到变量对 `\(\left(Y_{i g}, X_{i g}\right)\)`。 - 此时,可以定义模型为: `$$\begin{align} \begin{aligned} Y_{i g} &=m\left(X_{i g}\right)+e_{i g} \\ \mathbb{E}\left[e_{i g} \mid \boldsymbol{X}_{g}\right] &=0 \end{aligned} \end{align}$$` > 其中: - `\(\boldsymbol{X}_g\)`是 `\(X_{i g}\)`的堆栈形式(stacked)。 - 并假定群组之间是相互独立的。 --- ### 4.1 基本原理:矩阵表达 进一步地,我们定义如下: - `\(\boldsymbol{X}_g\)`是 `\(X_{i g}\)`的堆栈形式(stacked)。 - `\(\boldsymbol{Y}_g\)`是 `\(Y_{i g}\)`的堆栈形式(stacked)。 - `\(\boldsymbol{Z}_g(x)\)`是 `\(z_{i g}\)`的堆栈形式(stacked),其中: `$$\begin{align} Z_{i g}(x)=\left(\begin{array}{c} 1 \\ X_{i g}-x \end{array}\right) \end{align}$$` - `\(\boldsymbol{K}_{g}(x)=\operatorname{diag}\left\{K\left(\frac{X_{i g}-x}{h}\right)\right\}\)` --- ### 4.1 基本原理:LL估计 `$$\begin{align} \begin{aligned} \widehat{\beta}(x) &=\left(\sum_{g=1}^{G} \sum_{i=1}^{n_{g}} K\left(\frac{X_{i g}-x}{h}\right) Z_{i g}(x) Z_{i g}(x)^{\prime}\right)^{-1}\left(\sum_{g=1}^{G} \sum_{i=1}^{n_{g}} K\left(\frac{X_{i g}-x}{h}\right) Z_{i g}(x) Y_{i g}\right) \\ &=\left(\sum_{g=1}^{G} \boldsymbol{Z}_{g}(x)^{\prime} \boldsymbol{K}_{g}(x) \boldsymbol{Z}_{g}(x)\right)^{-1}\left(\sum_{g=1}^{G} \boldsymbol{Z}_{g}(x)^{\prime} \boldsymbol{K}_{g}(x) \boldsymbol{Y}_{g}\right) . \end{aligned} \end{align}$$` - 其中LL的估计量 `\(\hat{m}(x)=\widehat{\beta}_{1}(x)\)`等于上述估计的**截距项**。 --- ### 4.1 基本原理:删组回归法及其预测误差 为了得到**预测误差**(prediction error),我们可以采用**删组回归法**(delete cluster regression)进行遍历估计: `$$\begin{align} \widetilde{\beta}_{(-g)}(x)=\left(\sum_{j \neq g} \boldsymbol{Z}_{j}(x)^{\prime} \boldsymbol{K}_{j}(x) \boldsymbol{Z}_{j}(x)\right)^{-1}\left(\sum_{j \neq g} \boldsymbol{Z}_{j}(x)^{\prime} \boldsymbol{K}_{j}(x) \boldsymbol{Y}_{j}\right) \end{align}$$` 然后得到 `\(m(x)\)`的估计量 `\(\tilde{m}_{1}(x)=\tilde{\beta}_{1(-g)}(x)\)`,并进一步得到观测个体 `\(ig\)`的**删组预测误差**(delete-cluster prediction error): `$$\begin{align} \widetilde{e}_{i g}=Y_{i g}-\widetilde{\beta}_{1(-g)}\left(X_{i g}\right) \end{align}$$` --- ### 4.1 基本原理:条件方差 与之前类似,我们可以得到条件方差: `$$\begin{align} \boldsymbol{V}_{\widehat{\beta}}(x)=\left(\sum_{g=1}^{G} \boldsymbol{Z}_{g}(x)^{\prime} \boldsymbol{K}_{g}(x) \boldsymbol{Z}_{g}(x)\right)^{-1}\left(\sum_{g=1}^{G} \boldsymbol{Z}_{g}(x)^{\prime} \boldsymbol{K}_{g}(x) \boldsymbol{S}_{g}(x) \boldsymbol{K}_{g}(x) \boldsymbol{Z}_{g}(x)\right)\left(\sum_{g=1}^{G} \boldsymbol{Z}_{g}(x)^{\prime} \boldsymbol{K}_{g}(x) \boldsymbol{Z}_{g}(x)\right)^{-1} \end{align}$$` - 其中: `\(\boldsymbol{S}_{g}=\mathbb{E}\left[\boldsymbol{e}_{g} \boldsymbol{e}_{g}^{\prime} \mid \boldsymbol{X}_{g}\right]\)`, 因此,上述理论协方差矩阵可以通过 `\(\boldsymbol{e}_{g}\boldsymbol{e}_{g}^{\prime}\)`的估计量 `\(\widetilde{\boldsymbol{e}}_{g}\widetilde{\boldsymbol{e}}_{g}^{\prime}\)`计算得出,也即: `$$\begin{align} \widehat{\boldsymbol{V}}_{\widehat{\beta}}(x)=\left(\sum_{g=1}^{G} \boldsymbol{Z}_{g}(x)^{\prime} \boldsymbol{K}_{g}(x) \boldsymbol{Z}_{g}(x)\right)^{-1}\left(\sum_{g=1}^{G} \boldsymbol{Z}_{g}(x) \boldsymbol{K}_{g}(x) \widetilde{\boldsymbol{e}}_{g} \widetilde{\boldsymbol{e}}_{g}^{\prime} \boldsymbol{K}_{g}(x) \boldsymbol{Z}_{g}(x)\right)\left(\sum_{g=1}^{G} \boldsymbol{Z}_{g}(x) \boldsymbol{K}_{g}(x) \boldsymbol{Z}_{g}(x)\right)^{-1} . \end{align}$$` .footnote[ `\(\widehat{m}(x)\)`的方差也就是上述估计得到的协方差的第1个对角线元素。 ] --- ### 4.1 基本原理:交叉验证准则函数 与前面标准的**留一法**交叉验证略有不同,群组交叉验证准则函数可以表达为: `$$\begin{align} \mathrm{CV}(h)=\frac{1}{n} \sum_{g=1}^{G} \sum_{i=1}^{n_{g}} \widetilde{e}_{i g}^{2} \end{align}$$` 此时,最优CV谱宽将出现在上述函数的最小值处: `$$h_{\mathrm{CV}}=\underset{h \geq h_{\ell}}{\operatorname{argmin}} \mathrm{CV}(h)$$` --- exclude: true ## Code script: Hensen bruce fig 19-7a - 此代码运行时间较长(约3-5分钟),为提高渲染效率代码文件仅需运行一次。 - 如更新代码,请更改`eval = TRUE` --- name: case-score ### 4.2 (成绩案例):背景说明 .case[ **成绩案例**: - 案例基于([Duflo, Dupas, and Kremer, 2011](#bib-duflo2011))学生排位追踪(percentile tracking)对成绩(testscore)影响的数据集。其中**学生排位追踪**为分位数变量(1-100之间)。我们这里重点分析其中的**子样本数据**(女生、实施了学生排位追踪),样本数为n=1487。 - 关注的问题:**学生排位跟踪**(`X=percentile`)对**学生成绩**(`Y =testscore`)的非参数回归估计,并且我们对根据**学生所在学校**(school ID)对样本进行了分组。 - 后面的分析中,我们会重点划定观测窗口为:**成绩**范围 `\([0,40]\)`。非参数估计中我们会采用基于**高斯核函数**(Gaussian Kernel)的局部线性回归LL。 ] --- ### 4.2(成绩案例)样本数据集 .pull-left[
] .pull-right[ - 样本数据的描述性统计(未分群组)如下: ``` schoolid Y X Length:1487 Min. : 0 Min. : 0 Class :character 1st Qu.: 7 1st Qu.: 28 Mode :character Median :13 Median : 51 Mean :14 Mean : 51 3rd Qu.:21 3rd Qu.: 75 Max. :40 Max. :100 ``` ] --- ### 4.2(成绩案例)样本数据集:群组描述性统计
--- ### 4.2(成绩案例)样本数据散点图 <img src="seminar-RDD-part01_files/figure-html/unnamed-chunk-72-1.svg" height="450" style="display: block; margin: auto;" /> --- ### 4.2(成绩案例)参考谱宽的计算:多项式回归 下面我们按前述步骤来计算参考谱宽值 `\(h_{rot}\)`: - **步骤1**:根据案例数据集,设定权重取值范围 `\(\left\{\xi_{1} \leq x \leq \xi_{2}\right\}=\{0,100\}\)` - **步骤2**:构建多项式回归 `$$\begin{align} \begin{split} Y_i=&+\beta_{1}+\beta_{2}X_i+\beta_{3}X^2_i+\beta_{4}X^3_i+\beta_{5}X^4_i+u_i \end{split} \end{align}$$` > 直接使用OLS进行估计,得到估计方程: `$$\begin{equation} \begin{alignedat}{999} &\widehat{Y}=&&+6.815543&&+0.078095X_i&&+0.004483X^2_i&&-0.000100X^3_i&&+0.000001X^4_i\\ &(s)&&(1.2764)&&(0.1656)&&(0.0065)&&(0.0001)&&(0.0000) \end{alignedat} \end{equation}$$` > 进而得到拟合值 `\(\widehat{m}(x)\)`及其二阶导 `\(\widehat{m}^{\prime \prime}(x)\)`及残差 `\(\hat{\epsilon}\)` `$$\begin{align} \widehat{m}(x) &=6.815543+0.078095x_i+0.004483x^2_i-0.000100x^3_i+0.000001x^4_i \\ \widehat{m}^{\prime \prime}(x) &=2\times0.004483-6\times0.000100x^3_i+12\times0.000001x^2_i \\ \end{align}$$` --- ### 4.2(成绩案例)参考谱宽的计算:结果 - 步骤3:利用上述估计结果计算 `$$\begin{align} \widehat{B}=\frac{1}{n} \sum_{i=1}^{n}\left(\frac{1}{2} \hat{m}^{\prime \prime}\left(X_{i}\right)\right)^{2} \mathbb{1}\left\{\xi_{1} \leq X_{i} \leq \xi_{2}\right\} = 2 098.61641815 \end{align}$$` - 步骤4:多项式模型的回归误差方差 `$$\hat{\sigma}^2 = \frac{\sum{\hat{\epsilon}^2}}{n-q-1}=66.4390$$` - 步骤5:根据上述全部结果计算得到经验谱宽: `$$\begin{align} h_{\text {rot }}&=0.58\left(\frac{\widehat{\sigma}^{2}\left(\xi_{2}-\xi_{1}\right)}{n \widehat{B}}\right)^{1 / 5} = 0.58 \times \left(\frac{66.4390\times \left(1-0\right)}{1487 \times 2 098.616418148}\right)^{1 / 5} = 6.7463 \end{align}$$` ??? 这里对追踪排位X进行了尺度变换,转换为百分数 --- ### 4.2(成绩案例)交叉验证谱宽的计算:规则 - 步骤1:设定**经验谱宽** `\(h_{rot}=6.7463\)`作为**初始值**。。 - 步骤2:设定**调参谱宽**(tuning bandwidth)。 > - 人为选定谱宽范围: `\([4, 20]\)`。 - 给定范围内的搜寻总数为 `\(n=202\)`。则待评估序贯值为 `\(h\in (4.0000, 4.0796, 4.1592, 4.2388, 4.3184, \cdots,19.7612, 19.8408, 19.9204, 20.0000)\)`。 --- ### 4.2(成绩案例)交叉验证谱宽的计算:CV计算表
--- ### 4.2(成绩案例)交叉验证谱宽的计算:谱宽与CV变化(对比) <img src="seminar-RDD-part01_files/figure-html/unnamed-chunk-76-1.svg" height="450" style="display: block; margin: auto;" /> --- ### 4.2(成绩案例)交叉验证谱宽的计算:谱宽与CV变化(对比) - 为了使得两种方法具有可比较性,此处对CV值做了去最小化的尺度变换,也即 `\(cv^{\ast} = cv - min(cv)\)`。 - 常规的(未分组)局部线性LL交叉验证最优谱宽结果为12.2786,且CV函数相对比较陡峭;而群组化的(按学校分组)局部线性LL交叉验证最优谱宽结果为6.2289,并且在 `\([5,11]\)`之间CV函数值表现得比较平稳。 --- exclude: true ## Code script: Hensen bruce fig 19-7b - 此代码运行时间较长(约3-5分钟),为提高渲染效率代码文件仅需运行一次。 - 如更新代码,请更改`eval = TRUE` --- ### 4.2(成绩案例)CEF m(x)估计:基于群组CV最优谱宽(计算表)
- 我们采用群组化的(按学校分组)局部线性LL交叉验证最优谱宽结果为6.2289进行CEF估算 `\(\widehat{m}(x)\)` --- ### 4.2 (成绩案例)CEF m(x)估计:基于常规局部线性LL方法 <img src="seminar-RDD-part01_files/figure-html/unnamed-chunk-80-1.svg" height="450" style="display: block; margin: auto;" /> --- ### 4.2 (成绩案例)CEF m(x)估计:基于群组局部线性LL方法 <img src="seminar-RDD-part01_files/figure-html/unnamed-chunk-81-1.svg" height="450" style="display: block; margin: auto;" /> --- ### 4.2 (成绩案例)CEF m(x)估计:两种局部线性LL方法对比 <img src="seminar-RDD-part01_files/figure-html/unnamed-chunk-82-1.svg" height="450" style="display: block; margin: auto;" /> --- ### 4.2 (成绩案例)方差估计:计算方差、标准差 (4)直接使用CEF估计中**群组LL交叉验证**最优谱宽<sup>a</sup> `\(h=6.2289\)`进行局部线性LL估计,并利用删组法计算得到预测误差 `\(\tilde{\boldsymbol{e}}_g\)`,并最终分别得到未分组的和群组化的协方差矩阵(见下式),从而得到CEF估计值的方差和标准差(见后面附表)。 `$$\begin{align} \widehat{\boldsymbol{V}}_{\widehat{\beta}}(x)=\left(\sum_{g=1}^{G} \boldsymbol{Z}_{g}(x)^{\prime} \boldsymbol{K}_{g}(x) \boldsymbol{Z}_{g}(x)\right)^{-1}\left(\sum_{g=1}^{G} \boldsymbol{Z}_{g}(x) \boldsymbol{K}_{g}(x) \widetilde{\boldsymbol{e}}_{g} \widetilde{\boldsymbol{e}}_{g}^{\prime} \boldsymbol{K}_{g}(x) \boldsymbol{Z}_{g}(x)\right)\left(\sum_{g=1}^{G} \boldsymbol{Z}_{g}(x) \boldsymbol{K}_{g}(x) \boldsymbol{Z}_{g}(x)\right)^{-1} . \end{align}$$` .footnote[ <sup>a</sup> 这里我们没有再次评估条件方差估计中的最优谱宽,而是简单直接地使用了CEF估计的交叉验证最优谱宽。但是我们还是要注意,二者的最优谱宽可以完全不相同! ] --- ### 4.2 (成绩案例)方差估计:计算方差估计值(附表)
--- ### 4.2 (成绩案例)置信区间和置信带 (5)进一步计算**群组局部线性估计**下的**逐点置信区间**(Pointwise Confidence Interval)(见前面附表),并得到置信带(见后面附图)。 `$$\begin{align} \widehat{m}(x) \pm z_{1-\alpha/2}(n-1) \cdot \sqrt{\widehat{V}_{\widehat{m}(x)}}\\ \widehat{m}(x) \pm 1.96 \sqrt{\widehat{V}_{\widehat{m}(x)}} \end{align}$$` --- ### 4.2(成绩案例)置信区间和置信带(附图) <img src="seminar-RDD-part01_files/figure-html/unnamed-chunk-84-1.svg" height="450" style="display: block; margin: auto;" /> ??? --- ### 4.2(成绩案例)置信区间和置信带:OLS回归结果(附录) 作为比较,我们还简单使用了全样本OLS估计得到估计值(绿色) - 全局OLS回归模型为 `$$\begin{align} \begin{split} Y_i=&+\beta_{1}+\beta_{2}X_i+u_i \end{split} \end{align}$$` - 直接使用OLS进行估计,得到估计结果: `$$\begin{equation} \begin{alignedat}{999} &\widehat{Y}=&&+5.9789&&+0.1593X_i\\ &(s)&&(0.4472)&&(0.0077)\\ &(t)&&(+13.37)&&(+20.74)\\ &(over)&&n=1487&&\hat{\sigma}=8.2000\\ &(fit)&&R^2=0.2247&&\bar{R}^2=0.2242\\ &(Ftest)&&F^*=430.35&&p=0.0000 \end{alignedat} \end{equation}$$` --- layout: false class: center, middle, duke-orange,hide_logo name: reference-ch02 # 本章参考文献 --- layout: true <div class="my-header-h2"></div> <div class="watermark1"></div> <div class="watermark2"></div> <div class="watermark3"></div> <div class="my-footer"><span>huhuaping@    <a href="#chapter"> 第00章 课程说明 </a>                            <a href="#reference-ch02"> 本章参考文献 </a> </span></div> --- class: remark-slide-content.roomy ### 参考文献(References):1/2 .page-font-22[ <a name=bib-duflo2011></a>[Duflo, E., P. Dupas, and M. Kremer](#cite-duflo2011) (2011). "Peer Effects, Teacher Incentives, and the Impact of Tracking: Evidence from a Randomized Evaluation in Kenya". In: _American economic review_ 101.5, pp. 1739-74. ] --- layout:false background-image: url("../pic/thank-you-gif-funny-fan.gif") class: inverse, center # 本章结束