Stata+R：分位数回归一文读懂

发布时间：2020-11-18 阅读 26541

Stata 连享会主页 || 视频 || 推文

温馨提示： 定期清理浏览器缓存，可以获得最佳浏览体验。

✌ 课程详情： https://gitee.com/arlionn/Course | lianxh.cn

✨ 课程主页： https://gitee.com/arlionn/Course

谢雁翔 (南开大学)，xyxmask1995@163.com
钟舜斌 (北京工商大学)，1315829300@qq.com

编者按: 本文部分参考了游万海老师的分位数回归讲义和陈强老师的《高级计量经济学及Stata应用》，特此致谢！

1. 引言

在此前的推文中，我们对分位数回归和面板分位数回归都做过简单介绍，参见：

本文通过几篇论文的实操对分位数回归进行更为全面的介绍，内容涉及：分位数回归的基本思想、面板分位数回归、边际效应估计及图示等。

1.1 均值回归与条件分布

从实际来看，分位数信息也十分重要，比如，在评估某项干预对受众群体的影响时，我们不但希望了解干预的「平均」影响，更希望掌握干预对位于特征分布不同位置 (分布末端或顶端) 人群的「异质性」影响。

1.2 分位数回归

针对样本数据「异质性」特征，常用做法是根据数据特征进行「分组回归」，但这样的做法会导致「样本数据的损失」。为此，Koenker 和 Bassett (1978) 提出「分位数回归 (Quantile Regression, QR)」，并使用残差绝对值的加权平均 (如， $\sum_{i = 1}^{n} | e_{i} |$ ) 作为最小化的目标函数，尽可能减小极端值的影响。

由此可知，通过设置不同的分位点，分位数回归模型可以全面的刻画解释变量与被解释变量之间的关系。此外，相比于普通的线性回归 (均值回归)，分位数回归的估计结果对「偏态、多峰和异常值数据」更为稳健。

2. 分位数回归初识

对一个随机变量 $X$ 和任意一个 $0$ 到 $1$ 之间的数 $τ$ ，如果 $X$ 的取值 $x$ 满足 $prob (X \leq x) = τ$ ，则 $x$ 是 $X$ 的 $τ$ 分位数。上述过程语言表述为，在某个样本集中，从小至大排列之后，小于某值的样本子集占总样本集的比例。

有一组数据 $R = [3, 5, 10]$ ，如何找到一个数 $n$ ，使得其和 $R$ 中的元素尽可能的接近？—— 求数据的「集中趋势: 平均值」。更一般地，设 $μ$ 为最靠近该组数据的中心，则最小化残差平方和就是样本均值。

d (μ) = \sum_{i = 1}^{k} {(y_{i} - μ)}^{2}

即：

d (μ) = k μ^{2} - 2 μ \sum_{i = 1}^{k} (y_{i}) + \sum_{i = 1}^{k} (y_{i}^{2})

对 $μ$ 求偏导，得：

\frac{d}{d μ} d (μ) = 0

可得：

\frac{d}{d μ} d (μ) = 2 k μ - 2 \sum_{i = 1}^{k} (y_{i})

μ = \frac{\sum_{i = 1}^{k} (y_{i})}{k}

将标量 $μ$ 换成自变量 $X$ 的线性方程 $μ (X, β)$ ，则有:

min_{β} \sum_{i = 1}^{k} {(y_{i} - μ (x_{i}, β))}^{2}

令 $μ (x_{i}, β) = β^{'} x_{i}$ ，上式则变换为：

min_{β} \sum_{i = 1}^{k} {(y_{i} - β^{'} x_{i})}^{2}

对应参数 $β$ 的解即为线性回归模型 $y_{i} = β^{'} x_{i} + ε_{i}$ 的估计值。

(1) 如果损失函数定义为二次函数 $ρ (μ) = μ^{2}$ ，那么 $\hat{β} = min \sum_{i = 1}^{k} ρ (y_{i} - β^{'} x_{i})$ ，即 $\hat{β} = min \sum_{i = 1}^{k} {(y_{i} - β^{'} x_{i})}^{2}$ 对应的解则为最小二乘估计的解。

(2) 如果损失函数定义为 $ρ (μ) = | μ |$ ，那么 $\hat{β} = min \sum_{i = 1}^{k} ρ (y_{i} - β^{'} x_{i})$ ，即 $\hat{β} = min \sum_{i = 1}^{n} | y_{i} - β^{'} x_{i} |$ 对应的解则为最小一乘估计的解。

(3) 如果损失函数定义为:

ρ_{τ} (u) = {\begin{cases} τ u, & u \geq 0 \\ (τ - 1) u, & u < 0 \end{cases}

则：

\hat{β} (τ) = \arg min_{β \in R_{p}} \sum_{i = 1}^{n} ρ_{τ} (y_{i} - β^{'} (τ) x_{i})

即为对应下式的解：

min_{β_{τ}} \sum_{y_{i} \geq β^{'} (τ) x_{i}} τ (y_{i} - β^{'} (τ) x_{i}) + \sum_{y_{i} < β^{'} (τ) x_{i}} (τ - 1) (y_{i} - β^{'} (τ) x_{i})

损失函数就是上述例子中模型所表现的误差，最小化损失函数的过程其实是通过损失函数反过来优化模型参数。设分位点 $τ$ = 0.9，则损失函数为:

τ (d 2) + (1 - τ) (| d 1 + d 3 |) = 0.9 * 0.4 + 0.1 * (| - 1.3 + - 0.4 |) = 0.53

当 $τ$ =0.5 时候就是最小一乘回归。根据损失函数的形式，我们知道最小二乘回归和最小一乘回归的损失函数图像是关于 $Y$ 轴对称的，但是其他分位数的损失函数不一定对称，即对不同的分位点，高估和低估其预测值会有不同的损失权重。例如，对于0.25 分位数，高估 $Y$ 的损失权重为 0.75, 而低估的权重为 0.25；对于 0.75 分位数，则情况正好相反。可以直观理解为，越高的分位数估计，越倾向于高估预测值，越低的分位数则倾向于低估预测值。

在实际运用中，选择不对称的分位数损失函数往往可以反映因变量的分位数分布，不同分位数的分布可能由于各种原因存在差异。目前，分位回归模型的估计算法主要有单纯形法、内点算法和平滑算法。

3. 分位数回归模型与 Stata 实现

3.1 生成随机模拟数据

clear
set seed 1  //设置种子值，保证结果可重复
qui set obs 100  //设置100观测值
egen x = seq(), from(1) to(100)  //生成1到100序列
gen sig  = 0.1 + 0.05*x //sig代表方差，方差随着x在变化，不是同方差，异方差不对称 
scalar b_0 = 6                                
scalar b_1 = 0.1                              
gen e = rnormal(0,sig) //生成误差项                      
gen y = b_0 + b_1*x + e                  
twoway scatter y x, plotregion(fcolor(white)) graphregion(fcolor(white))  
*散点图呈现喇叭状，非均匀分布考虑分位数回归，若均匀分布不考虑分位数回归

从上图可以看出，随着 $x$ 的增大， $y$ 的波动性越大，这违反了线性回归模型 (OLS) 的同方差假定，此时普通线性回归模型不再适用。

3.2 分位数模型估计及 Stata 实现

分位数回归模型的基本步骤：

绘制散点图，判断数据分布 (对称/非对称) 是否适合于用分位数回归模型；
建模估计 qreg，sqreg，bsqreg，qreg2(聚类稳健标准误)；
系数差异性检验 wald 检验；
系数可视化 grqreg。

在 Stata 的进行分位数回归估计的命令中，sqreg 较为常用，详见 help sqreg。

*计算单个分位点
sysuse auto, clear
qreg price weight length foreign, quantile(.25)
est store m0

*多个分位点，从0.3-0.9，可用qreg的循环语句也可使用sqreg命令
forvalues i = 0.3(0.3)0.9{
  qreg price weight length foreign, quantile(`i')
  local j = `i' * 10
  local j = int(`j')
  est store m`j'
}   

local m "0.25RQ 0.3RQ 0.6RQ 0.9RQ"
esttab m0 m3 m6 m9, mtitle(`m') nogaps star(* 0.1 ** 0.05 *** 0.01)

----------------------------------------------------------------------------
                      (1)             (2)             (3)             (4)   
                   0.25RQ           0.3RQ           0.6RQ           0.9RQ   
----------------------------------------------------------------------------
weight              1.832***        1.705**         5.975***        9.666***
                   (2.89)          (2.59)          (4.00)          (4.53)   
length              2.846           5.090          -100.1*         -248.8***
                   (0.13)          (0.23)         (-1.96)         (-3.40)   
foreign            2209.9***       2148.6***       3731.0***       4164.4***
                   (5.24)          (4.90)          (3.75)          (2.93)   
_cons             -1879.8         -1843.9          5901.2         25191.0***
                  (-0.76)         (-0.72)          (1.01)          (3.02)   
----------------------------------------------------------------------------
N                      74              74              74              74   
----------------------------------------------------------------------------
t statistics in parentheses
* p<0.1, ** p<0.05, *** p<0.01

3.3 Wald 检验

不同分位数下系数差异检验用 Wald 检验，具体如下：

假设 $b$ 为系数， $V$ 为方差协方差矩阵。现在需要做如下假设检验：

R b = r

相应的 wald 检验统计量为

W = (R b - r)^{'} {({R V R}^{'})}^{- 1} (R b - r)

这个统计量服从卡方分布

W \sim χ_{q}^{2}

其中：q 为待检验的线性假设个数。

set seed 1   //设置种子值
sysuse auto, clear
*sqreg不用循环语句，默认情况下reps为20，bootstrap100次
sqreg price weight length foreign, quantile(.25 .5 .75) reps(100)
*wald检验
test [q25]weight=[q75]weight //检验0.25分位点处与0.75分位点处是否存在显著差异,拒绝原假设则存在差异

 (1)  [q25]weight - [q75]weight = 0
       F( 1, 70) = 11.62
         Prob > F = 0.0011

3.4 系数可视化

*ssc install grqreg,all replace  //安装外部命令
sysuse auto, clear
qreg price mpg headroom
grqreg, cons ci ols olsci title(Fig.1a Fig.1b Fig.1c)  //cons表示常数项绘图，ci表示置信区间，ols表示绘出普通ols，olsci表示ols回归置信区间

Median regression                                   Number of obs =         74
  Raw sum of deviations  71102.5 (about 4934)
  Min sum of deviations 64549.83                    Pseudo R2     =     0.0922

------------------------------------------------------------------------------
       price |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         mpg |  -127.3333   60.78234    -2.09   0.040    -248.5299   -6.136791
    headroom |  -130.6667   415.6721    -0.31   0.754    -959.4933      698.16
       _cons |   8272.333   2158.205     3.83   0.000     3968.995    12575.67
------------------------------------------------------------------------------

上图分别为 cons (常数项) 、mpg 与 headroom 的分位数估计结果，横轴的虚线为 ols 估计结果及其置信区间。但也存在问题，比如纵轴坐标重叠和背景非白色。

4. 面板分位数回归

面板分位数回归，顾名思义即将截面分位数回归模型扩展至面板数据框架下，考虑了研究单元之间的个体异质性。特别注意的是，面板分位数回归为非线性关系，对于个体效应的估计不能通过 i.id 虚拟变量的形式控制。

4.1 Koenker (2004)

惩罚分位 (Penalized quantile regression with fixed effects)，应用了LASSO的思想:

min_{(α, β)} \sum_{k = 1}^{q} \sum_{j = 1}^{n} \sum_{i = 1}^{m_{i}} w_{k} ρ_{τ_{k}} (y_{i j} - α_{i} - x_{i j}^{⊤} β (τ_{k})) + λ \sum_{i = 1}^{n} | α_{i} |

缺点: (1) 当 $n$ 比较大时，计算费时；(2) $λ$ 的选择会影响估计结果。

实现：R 中的 rqpd 包进行估计，rqpd不支持 R 3.6.3 及以后版本，建议直接在 R 中运行，更多内容见「rqpd-package: Regression quantiles for panel data」，也可以使用 Stata 外部命令rcall 运行，详细请参考「Rcall：Stata 与 R 的无缝对接」。

温馨提示: 文中链接在微信中无法打开，请点击「阅读原文」。

## rcall 可以从Stata中运行R命令，为简便运算可以直接复制粘贴进R运算
install.packages("rqpd", repos="http://R-Forge.R-project.org")
install.packages("plm", repos="http://R-Forge.R-project.org")
library(rqpd)   ##加载面板分位数回归包
library(plm)    ##加载面板数据回归包

data(Produc,package="plm")   ##调用plm包自带数据集进行面板估计
mydata = Produc
m <- 17    ##17年
n <- 48    ##美国48州
s <- as.factor(rep(1:n,rep(m,n)))   ##声明面板数据结构变量
x <-cbind(mydata$pcap,mydata$hwy)   ##数据中取出解释变量x
y <- mydata$unemp                   ##数据中取出被解释变量y

##rqpd直接运算公式，个体效应惩罚系数lambda为1,系统默认分位点是0.25,0.5,0.75
fit <- rqpd(y ~x | s, panel(lambda = 1))  
sfit <- summary(fit)
sfit
sfit$coefficients  ##能一并取出估计的个体效应系数

##依次定义九个分位点同时定义tau和tauw
fit1 <- rqpd(y ~x | s, panel(taus=seq(0.1,0.9,by=0.1),tauw=rep(1/9,9),lambda = 1))   
sfit1 <- summary(fit1)
sfit1

4.2 Canay (2011)

两阶段估计法 (Two-step estimator)，分位数回归为非线性运算，面板数据不能通过组内去心的方法解决。

Y_{i t} = X_{i t}^{'} θ (τ) + α_{i} + e_{i t} (τ), Q_{e_{i t} (τ)} (τ | X_{i}) = 0

思想：先消除个体效应，再进行估计。与 Koenker (2004) 的相同点：将 $α_{i}$ 设置为 location，即不随分位点改变。实现：作者主页 R 程序，详见「程序」。

温馨提示: 文中链接在微信中无法打开，请点击「阅读原文」。

4.3 Machado and Silva (2019)

Q_{Y} (τ | X_{i t}) = (α_{i} + δ_{i} q (τ)) + X_{i t}^{'} β + Z_{i t}^{'} γ q (τ)

优点：(1) 个体效应允许随着分位点在变动 (location-scale)；(2) 可以纳入其他内生解释变量；(3) 计算速度较快。

缺点: 假定变量系数与分位点之间关系是单调的（要不单调上升，要不单调下降），即不存在分位交叉 (No Quantile-Crossing)，画图无意义。

实现： Stata 中 xtqreg 进行估计，详细参见 help xtqreg。

关于 xtqreg 系数可视化，请参考以下代码：

clear    //清空数据 
set obs 100    //设置100观测值
egen country_id = seq(),b(10)    //每10个重复构造整数序列，生成变量country_id
bysort country_id: gen time = _n    //根据分组country_id生成时间变量

set seed 12345   //设置种子值确保可重复
gen y = uniform()   //生成被解释变量y服从正态分布
forv i = 1/8{
gen z`i' = uniform()
}   //循环语句生成8个解释变量

local x1 "z1 z2 z3 z4 z5 z6 z7 z8"   //局部宏定义所有解释变量
*list  `x1'   //列出所有被解释变量

local v_num : word count `x1'   //利用扩展函数 word count string 得到x1中变量的个数
di `v_num'   //v_num表示最大的循环个数,8个
est clear
mat result1 = J(`v_num',1,.) //生成空矩阵result1,8行1列
matrix rownames result1 = `x1'   //矩阵行命名为对应的被解释变量
mat list result1   //列出空矩阵

tempname a b b1   //设定临时矩阵变量名
forvalues quantile = 0.1(0.1)0.9 {   //从0.1-0.9没隔0.1开始循环 
    local wanted : di %2.1f `quantile'   //列出循环分位点长度设置为2，保留一位小数，定义为`wanted'
 local wanted = `wanted'*10    //`wanted'为每一分位点重复十次
    xtqreg y `x1', i(country_id) q(`quantile')   //面板分位数回归
    matrix `a' = r(table)   //对应各个分位点估计的`x1'系数
 mat list `a'	//列出面板分位数回归`a'矩阵
    matrix `b' = `a'[1,1..`v_num'] \ `a'[2,1..`v_num'] //仅取出面板分位数回归`a'矩阵的前两行构成`b'矩阵,第一行是系数。第二行是标准误，第四行是p值
 mat list `b'   //列出`b'矩阵
 mat `b1' = `b''   
 matrix colnames `b1' = b_q`wanted' p_q`wanted'   //`b1'为`b'矩阵的转置
 mat list `b1'   //列出`b1'矩阵
 matrix result1 = (result1, `b1')    //`b1'数据导入列出空矩阵result1
}

mat list result1
matselrc result1 cresult, c(2/19) //19列 = 2*分位个数+1列无用项（第一列），保留有用数据列
mat list cresult

local newnames : colfullnames cresult
di "`newnames'"
local cn : word count `newnames'
di `cn'
svmat cresult, names(reg)   //矩阵数据转换为变量
keep reg*

format reg* %4.3f   //设定格式为保留三位小数
drop if reg1==. //删除缺失值，行数为8,8个解释变量：变量为18，对应不同分位点的系数和标准误；

forv i = 1/9 {  /*9：分位个数*/
local k1 = 2*`i' - 1   //奇数列为系数
local k2 = 2*`i'       //偶数列为标准误
gen up_`i' = reg`k1' + 1.96*reg`k2'   //上界
gen low_`i' = reg`k1' - 1.96*reg`k2'  //下界
}

preserve 
keep reg1 reg3 reg5 reg7 reg9 reg11 reg13 reg15 reg17   //保留系数列
rename (reg1 reg3 reg5 reg7 reg9 reg11 reg13 reg15 reg17) (reg1 reg2 reg3 reg4 reg5 reg6 reg7 reg8 reg9)   //重命名
xpose,clear  //转置数据
renvars v*,  prefix(coef_)   //重命名加前缀coef_
save coef.dta,replace  //储存系数数据
restore   //preserve+restore运行不保存

keep up* low*   //保留上下界数据
xpose,clear   //转置数据
preserve
keep if mod(_n,2)==1   //保留奇数行作为上界数据
renvars v*,  prefix(up_)     //重命名加前缀up_
save up.dta,replace    //储存上界数据
restore

preserve   
keep if mod(_n,2)==0   //保留偶数行作为下界数据
renvars v*,  prefix(low_)   //重命名加前缀low_
save low.dta,replace   //储存下界数据
restore 

use up.dta,clear    //打开上界数据
merge 1:1 _n using low.dta   //一对一合并数据
keep if _merge==3   //合并匹配成功为_merge==3，保留合并成功行
drop _merge  //删除合并信息
merge 1:1 _n using coef.dta   //一对一合并数据
keep if _merge==3
drop _merge  //删除合并信息

matrix input myrvec = (0.10,0.20,0.30,0.40,0.50,0.60,0.70,0.80,0.90)
mat list myrvec
mat qq = myrvec'  //qq矩阵为myrvec的转置
mat list qq
svmat qq, names(qqp)  //矩阵数据转换为变量

forv i = 1/8{
twoway (line up_v`i' qqp,lc(black*1.4) lpattern(dash) xlabel(0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90)) ///
(line low_v`i' qqp,lc(black*1.4) lpattern(dash) xlabel(0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90))  ///
(line coef_v`i' qqp, lc(black*2)  lw(thick)) ,xtitle("Quantile") ytitle("Coefficient(z`i')") plotregion(margin(zero)) graphregion(color(white)) legend(off)
graph save g_`i'
}

graph combine g_1.gph g_2.gph g_3.gph g_4.gph g_5.gph g_6.gph g_7.gph g_8.gph

4.4 Powell (2015)

qregpd 估计是实现广义分位数估计 genqreg 的特例，其解决了由其他固定效应分位数估计引起的一个基本问题：包含个体固定效应会改变对处理变量的估计系数的解释。由于标准差的计算有时可能收极值的影响，可以使用 MCMC 方法或网格搜索方法来估计广义分位数回归。关于该命令更多详细介绍，请参考 help qregpd。相关代码如下：

*ssc install qregpd, replace
*ssc install moremata, replace

sysuse nlswork, clear
qregpd ln_wage tenure union, id(idcode) fix(south) 
estimates store QREG1
matrix list e(gamma)      //显示south固定效应估计值

qregpd ln_wage tenure union, id(idcode) fix(year)
estimates store QREG2
matrix list e(gamma)      //显示year时点固定效应估计值

*带工具变量的情形--采用MCMC优化
*tenure为内生变量，union为外生变量做自己的工具变量，ttl_exp wks_work为工具变量  
qregpd ln_wage tenure union, id(idcode) fix(year) optimize(mcmc) draws(1000) burn(100) arate(.5) noisy instruments(ttl_exp wks_work union)
mat list e(gamma)

*估计1/4分位数的结果--默认是0.5
qregpd ln_wage tenure union, q(0.25) id(idcode) fix(year) optimize(mcmc) draws(1000) burn(100) arate(.5) noisy instruments(ttl_exp wks_work union)
mat list e(gamma)

5. 更多参考资料

[1] 连享会直播课程：分位数回归 (游万海) -Link-
连享会直播：分位数回归
连享会推文：Stata：分位数回归简介
连享会推文：Stata: 面板分位数回归
[2] 分位数回归及Stata实现 -Link-
[3] You, Wan-Hai, et al. Democracy, Financial Openness, and Global Carbon Dioxide Emissions: Heterogeneity Across Existing Emission Levels. World Development, vol. 66, 2015, pp. 189–207. -Link-
[4] You, Wanhai, et al. Oil Price Shocks, Economic Policy Uncertainty and Industry Stock Returns in China: Asymmetric Effects with Quantile Regression. Energy Economics, vol. 68, 2017, pp. 1–18. -Link-
[5] You, Wanhai, et al. Twitter’s Daily Happiness Sentiment and the Predictability of Stock Returns. Finance Research Letters, vol. 23, 2017, pp. 58–64. -Link-
[6] Chuang, C.C., Kuan, C.M., Lin, H.Y., 2009. Causality in quantiles and dynamic stock return volume relations.Journal of Banking & Finance. 33, 1351–1360. -Link-
[7] Xiao, Z. 2009. Quantile cointegrating regression. Journal of Econometrics, 150(2), 248-260. -Link-
[8] White, H., Kim, T. H., & Manganelli, S. 2015. VAR for VaR: Measuring tail dependence using multivariate regression quantiles. Journal of Econometrics, 187(1), 169-188. -Link-
[9] Cho, J. S., Kim, T. H., & Shin, Y. 2015. Quantile cointegration in the autoregressive distributed-lag modeling framework. Journal of Econometrics, 188(1), 281-300. -Link-
[10] Koenker, R. (2004). Quantile Regression for Longitudinal Data, Journal of Multivariate Analysis, 91, 74–89. -Link-
[11] Canay, I.A. (2011). A Simple Approach to Quantile Regression for Panel Data, The Econometrics Journal, 14, 368–386. -Link-
[12] Machado, J.A.F. and Santos Silva, J.M.C. (2019), Quantiles via Moments, Journal of Econometrics, 213(1), pp. 145-173. -Link-
[13] Powell, David. 2015. Quantile Regression with Nonadditive Fixed Effects. RAND Labor and Population Working Paper. -Link-
[14] Powell, David. 2016. Quantile Treatment Effects in the Presence of Covariates. RAND Labor and Population Working Paper. -Link-

专题	嘉宾	直播/回看视频
⭐ 最新专题		因果推断, 空间计量，寒暑假班等
⭕ 数据清洗系列	游万海	直播, 88 元，已上线
研究设计	连玉君	我的特斯拉-实证研究设计，-幻灯片-
面板模型	连玉君	动态面板模型，-幻灯片-
面板模型	连玉君	直击面板数据模型 [免费公开课，2小时]