范文一:面板分位数例子
data(bwd)
cre.form <- dbirwt="" ~="" smoke="" +="" dmage="" +="" agesq="" +="">->
+ novisit + pretri2 + pretri3 | momid3 | smoke + + dmage + agesq
# CRE-M type fit:
crem.fit <- rqpd(cre.form,="" panel(method="cre" ),="" data="bwd)" summary(crem.fit)="">->
Call: rqpd(formula = cre.form, panel = panel(method = "cre"), data = bwd)
taus: [1] 0.25 0.50 0.75
Coefficients:
Value Std. Error t value Pr(>|t|) (Intercept)[0.25] 2440.23361 165.18256 14.77295 0.00000 smoke[0.25] -114.21319 31.65224 -3.60838 0.00031
-5.59702 17.87334 -0.31315 0.75417 dmage[0.25]
agesq[0.25] 0.42689 0.31758 1.34420 0.17891
novisit[0.25] -380.30477 121.27964 -3.13577 0.00172
-39.35954 17.63163 -2.23232 0.02561 pretri2[0.25]
pretri3[0.25] -5.33718 39.79378 -0.13412 0.89331 m.smoke[0.25] -167.86680 35.28719 -4.75716 0.00000 m.dmage[0.25] 55.15672 20.46185 2.69559 0.00704 m.agesq[0.25] -1.18906 0.35838 -3.31785 0.00091 (Intercept)[0.5] 2943.08621 146.53085 20.08510 0.00000 smoke[0.5] -118.96828 29.78308 -3.99449 0.00007 dmage[0.5] -37.64142 18.32136 -2.05451 0.03995 agesq[0.5] 0.97716 0.30761 3.17666 0.00149 novisit[0.5] -177.32668 65.09804 -2.72399 0.00646 pretri2[0.5] -23.94772 16.59886 -1.44273 0.14912 pretri3[0.5] -18.41896 32.05553 -0.57460 0.56557 m.smoke[0.5] -144.49955 35.51163 -4.06908 0.00005 m.dmage[0.5] 71.74920 21.34456 3.36148 0.00078 m.agesq[0.5] -1.44484 0.36035 -4.00960 0.00006 (Intercept)[0.75] 3367.51099 163.17049 20.63799 0.00000 smoke[0.75] -120.63708 34.06380 -3.54150 0.00040 dmage[0.75] -29.15072 23.12975 -1.26031 0.20758 agesq[0.75] 0.83485 0.39313 2.12360 0.03372 novisit[0.75] -21.21664 105.08383 -0.20190 0.84000 pretri2[0.75] 23.61037 20.99615 1.12451 0.26082 pretri3[0.75] -39.49583 36.96971 -1.06833 0.28539 m.smoke[0.75] -138.38759 41.09140 -3.36780 0.00076 m.dmage[0.75] 54.88645 25.42125 2.15908 0.03086 m.agesq[0.75] -1.14623 0.43581 -2.63012 0.00855
crem.fit <- rqpd(cre.form,="" panel(method="cre" ,taus="c(0.1," 0.25,="" 0.5,="" 0.75,="" 0.9),="">->
tauw=rep(1/5, 5)), data=bwd)
summary(crem.fit)
Call: rqpd(formula = cre.form, panel = panel(method = "cre", taus = c(0.1,
0.25, 0.5, 0.75, 0.9), tauw = rep(1/5, 5)), data = bwd)
taus: [1] 0.10 0.25 0.50 0.75 0.90
Coefficients:
Value Std. Error t value Pr(>|t|) (Intercept)[0.1] 1971.78927 200.77800 9.82074 0.00000 smoke[0.1] -189.70646 39.39055 -4.81604 0.00000 dmage[0.1] 37.59929 27.78228 1.35336 0.17597 agesq[0.1] -0.23533 0.45634 -0.51569 0.60608
-602.65402 153.40321 -3.92856 0.00009 novisit[0.1]
pretri2[0.1] -76.52337 27.04753 -2.82922 0.00467 pretri3[0.1] -2.43651 51.46743 -0.04734 0.96224
-55.36162 51.18102 -1.08168 0.27941 m.smoke[0.1]
m.dmage[0.1] 23.35438 30.98212 0.75380 0.45098
m.agesq[0.1] -0.71575 0.51149 -1.39933 0.16174 (Intercept)[0.25] 2440.23361 165.83251 14.71505 0.00000 smoke[0.25] -114.21319 29.17432 -3.91485 0.00009 dmage[0.25] -5.59702 16.97241 -0.32977 0.74158 agesq[0.25] 0.42689 0.30057 1.42025 0.15556 novisit[0.25] -380.30477 136.15232 -2.79323 0.00523 pretri2[0.25] -39.35954 17.00846 -2.31411 0.02068 pretri3[0.25] -5.33718 42.80340 -0.12469 0.90077 m.smoke[0.25] -167.86680 33.27678 -5.04456 0.00000 m.dmage[0.25] 55.15672 20.52318 2.68753 0.00721 m.agesq[0.25] -1.18906 0.35250 -3.37321 0.00075 (Intercept)[0.5] 2943.08621 147.01143 20.01944 0.00000 smoke[0.5] -118.96828 28.98334 -4.10471 0.00004 dmage[0.5] -37.64142 19.85702 -1.89562 0.05803 agesq[0.5] 0.97716 0.34095 2.86596 0.00416 novisit[0.5] -177.32668 73.66796 -2.40711 0.01609 pretri2[0.5] -23.94772 16.09472 -1.48792 0.13679 pretri3[0.5] -18.41896 30.59545 -0.60202 0.54717 m.smoke[0.5] -144.49955 35.40814 -4.08097 0.00005 m.dmage[0.5] 71.74920 22.99920 3.11964 0.00181 m.agesq[0.5] -1.44484 0.39160 -3.68962 0.00023 (Intercept)[0.75] 3367.51099 171.13242 19.67781 0.00000 smoke[0.75] -120.63708 32.01933 -3.76763 0.00017
dmage[0.75] -29.15072 21.94964 -1.32807 0.18418 agesq[0.75] 0.83485 0.37618 2.21930 0.02648 novisit[0.75] -21.21664 106.20581 -0.19977 0.84166 pretri2[0.75] 23.61037 19.26091 1.22582 0.22029 pretri3[0.75] -39.49583 29.11850 -1.35638 0.17500 m.smoke[0.75] -138.38759 39.56134 -3.49805 0.00047 m.dmage[0.75] 54.88645 25.97660 2.11292 0.03463
m.agesq[0.75] -1.14623 0.44453 -2.57854 0.00993 (Intercept)[0.9] 3523.39409 203.46151 17.31725 0.00000 smoke[0.9] -120.42193 52.77827 -2.28166 0.02253 dmage[0.9] -51.20245 26.29222 -1.94744 0.05150 agesq[0.9] 1.26688 0.45622 2.77691 0.00550 novisit[0.9] 13.89590 77.00535 0.18045 0.85680 pretri2[0.9] 38.18016 24.91348 1.53251 0.12542 pretri3[0.9] -75.90907 50.53785 -1.50202 0.13311 m.smoke[0.9] -80.40990 61.71414 -1.30294 0.19262 m.dmage[0.9] 87.11635 28.98700 3.00536 0.00266 m.agesq[0.9] -1.75179 0.49083 -3.56903 0.00036
范文二:分位数回归 张晓峒
第15章 分位数回归模型
15.1 总体分位数和总体中位数
15.2 总体中位数的估计
15.3 分位数回归
15.4 分位数回归模型的估计
15.5 分位数回归模型的检验
15.6 分位数的计算与分位数回归的EViews操作
15.7 分位数回归的案例分析
以往介绍的回归模型实际上是研究被解释变量的条件期望。人们当然也关心解释变量与被解释变量分布的中位数,分位数呈何种关系。这就是分位数回归,它最早由Koenker和Bassett(1978)提出,是估计一组回归变量X与被解释变量Y的分位数之间线性关系的建模方法。
正如普通最小二乘OLS回归估计量的计算是基于最小化残差平方和一样,分位数回归估计量的计算也是基于一种非对称形式的绝对值残差最小化,其中,中位数回归运用的是最小绝对值离差估计(LAD,least absolute deviations estimator)。它和OLS主要区别在于回归系数的估计方法和其渐近分布的估计。在残差检验、回归系数检验、模型设定、预测等方面则基本相同。
分位数回归的优点是,(1)能够更加全面的描述被解释变量条件分布的全貌,而不是仅仅分析被解释变量的条件期望(均值),也可以分析解释变量如何影响被解释变量的中位数、分位数等。不同分位数下的回归系数估计量常常不同,即解释变量对不同水平被解释变量的影响不同。
另外,中位数回归的估计方法与最小二乘法相比,估计结果对离群值则表现的更加稳健,而
数且,分位数回归对误差项并不要求很强的假设条件,因此对于非正态分布而言,分位数回归系估计量则更加稳健。
15.1 总体分位数和总体中位数
在介绍分位数回归之前先介绍分位数和中位数概念。
对于一个连续随机变量y,其总体第τ分位数是y的定义是:y小于等于y的概率是τ,即 (τ)(τ)
τ = P( y ? y) = F(y) (τ)(τ)
其中P(,)表示概率,F(y) 表示y的累积(概率)分布函数(cdf)。 (τ)
比如y = 3,则意味着y ? 3的概率是0.25。且有 (0.25)
-1 y= F(y) (τ) (τ)
即F(y)的反函数是y。当τ=0.5时,y 是y的中位数。τ= 0.75时,y 是y的第3/4分位数,(τ)(τ)(τ)(τ)
τ= 0.25时,y是y的第1/4分位数。若y服从标准正态分布,y= 0,y =1.645,y =1.960。 (τ) (0.5) (0.95)(0.975)
另外,如果随机变量y的分布是对称的,那么其均值与中位数是相同的。当其中位数小于均值时,分布是右偏的。反之,分布是左偏的。
对于回归模型,被解释变量y对以X为条件的第τ分位数用函数y,X表示,其含义是:以t(τ)t
X为条件的y小于等于y,X的概率是τ。这里的概率是用y对X的条件分布计算的。且有 t(τ)tt
-1 y,X = F(y,X) (τ)t(τ)t
其中F(y,X) 是y在给定X条件下的累积概率分布函数(cdf)。则y,X称作被解释变量y对X(τ)tt(τ)tt''的条件分位数函数。而F (y,X)= f (y,X)则称作分位数概率密度函数。其中F(y,X)表示F(y,X)(τ)t(τ)t(τ)t(τ)t
1
对y,X求导。 (τ)t
15.2 总体中位数的估计
在介绍分位数回归之前,先来看中位数的估计和中位数回归。下面以连续变量为例介绍定理
15.1。
定理15.1
连续变量用y表示,其概率密度函数用f(y)表示,累计概率密度函数用F(y)表示,y的中位
数用y表示,则y与任一值,的离差绝对值的期望以, = y 时为最小。 E(y,,)(0.5)(0.5)
证明:
,, = E(y,,)-(y,,)f(y)dy,(y,,)f(y)dy,,-,,
,, = (15.1) -(y,,)dF(y),(y,,)dF(y),,-,,
bb,,fy(,),, 根据莱布尼兹公式,若,则有。令,则,f(y,,),y-,F(,),f(y,,)dyFdy(),,,,aa
bb,,(y-),,有。运用于式(15.1),得 F(),dy,-dy,,,a,a
,,,,,,,,(y)f(y)dy(y)f(y)dy,,,,(,)Ey,,t,-, == ,dF(y)-dF(y)-,,-,,,,,,,,
, = F,()-[1-dF(y)],F(,)-(1-F(,)),2F(,)-1,-,
,,(,)Eyt 式(15.1)求极小的一阶条件是= 0,即=0,。这意味着,等2F(,)-1F(,),0.5,,
中位数y于。 (0.5)
, = y (0.5)
与定理15.1等价的表述是以, = y(中位数)时为最小。因此,中位数回归估计y,,(0.5),
量可以通过最小绝对离差法(least absolute deviation, LAD)估计。其中X和,分别为(k,1)阶列向
量。
?, 同理,对于线性回归模型y= X ,, + u,通过求最小,估计,的中位数回归系y,Xβt tt(0.5),
???y数估计量,从而得到y的中位数回归估计量(X),Xβ。 βt(0.5)t(0.5)(0.5)
15.3 分位数回归ρ
? Koenker和Bassett(1978)证明,若用y表示y的分位数回归估计量,则对于以检查函数t(,)t
?ywy,,(check function)w为权数,ρy对任意值,的加权离差绝对值和只有在, =时,,t ,t(,)t,取得最小值。其中
TT
wy,, = ,(1,,)(y,,),,(y,,) (15.2) ,ttt,,,:,,:,,iytyii
2
,(0, 1)。据此,分位数回归可以通过加权的最小绝对离差和法(weighted least absolute deviation, ,
WLAD)进行估计。
? 根据式(15.2),对于线性回归模型y= X ,, + u, 求第,分位数回归方程系数的估计量的βt t(,)方法是求下式(目标函数)最小,
TT
?? Q,,(1,,)u,,u()(),,tt,,??,0,0uu,()t(,)t
TT??,, (15.3) (1,)(yXβ),(yXβ),,,,,,,,()()tt,,??,,,,,,::tyXtyXt,()t(,)
?表示第分位数回归方程对应的残差。,(0, 1)。第分位数的回归方程表达式是 其中,,,u()t,
??, = yXβ(,)t(,)
?其中X,,都是k,1阶列向量。称作分位数回归系数估计量,或最小绝对离差和估计量,估β(,)
计方法称作最小绝对离差和估计法。
当,=0.5时,式(15.3)变为
TTT???,,, QyXβyXβyXβ,,,,,,,0.5()0.5()0.5(0.5)(0.5)(0.5)ttt,,,??,1,,t,,:,:,tyXtyXt(0.5)t(0.5)
???,y=称作中位数回归方程,称作中位数回归系数估计量。 βXβ(0.5)t(0.5)(0.5)
? 一旦得到估计的分位数回归方程,就可以计算分位数回归的残差u。 ()t,
???, u,y,y,y- Xβ(,)tt(,)tt(,)
对一个样本,估计的分位数回归式越多,对被解释变量y条件分布的理解就越充分。以一元t
回归为例,如果用LAD法估计的中位数回归直线与用OLS法估计的均值回归直线有显著差别,则表明被解释变量y的分布是非对称的。如果散点图上侧分位数回归直线之间与下侧分位数回归t
直线之间相比,相互比较接近,则说明被解释变量y的分布是左偏倚的。反之是右偏倚的。对于t
不同分位数回归函数如果回归系数的差异很大,说明在不同分位数上解释变量对被解释变量的影响是不同的。
15.4 分位数回归模型的估计
由于目标函数(15.3)不可微,因此传统的对目标函数求导的方法不再适用。估计分位数回
?归方程参数的一种较好的方法是线性规划方法。 β(,)
基于Barrodale和Roberts (1973,以下简写为BR)提出的单纯形法(simplex algorithm),Koenker和D’Orey(1987)提出一种估计分位数回归系数的方法。EViews中应用的是上述算法的改进形式。
BR算法由于其非有效性和大样本下的一些非优良特性曾备受批评。Koenker和Hallock(2001) 以及Portnoy和Koenker(1997)通过模拟证实,与内点法(interior point method)等替代方法相比,BR算法的估计次数往往较多,大约是样本容量的平方次数。然而,改进的BR算法的估计次数在一定程度上是可以接受的,大约是样本容量的线性倍次数,在实际中是可以使用的。
3
分位数回归方程的BR算法原理略。
下面讨论分位数回归系数估计量的渐近分布。
在弱条件下,分位数回归系数渐近服从正态分布(Koenker, 2005)。回归系数的方差协方差矩阵的计算在分位数回归的系数估计中占有重要位置。其方差协方差矩阵的估计方法根据分位数密度函数是否与解释变量相关分为三种方法:
?误差项独立同分布(i.i.d.)假设下的直接估计方法。由Koenker和Bassett(1978)提出。
?误差项独立但不同分布(i.n.i.d.)条件下的直接估计方法。
?误差项独立同分布(i.i.d.)和独立但不同分布(i.n.i.d.)条件下都可使用的自举法。
(1)独立同分布假设下的参数渐近分布
Koenker和Bassett(1978)在独立同分布假设下得出分位数回归系数渐近服从正态分布,可以表述为在弱条件下:
2,1? ~ (15.5) n(,,,)N(0,,(1,,)sJ)(,)(,),()
其中
,XX,XXii (15.6) Jlim()lim(),,,n,,n,,TTi
,,1,1 (15.7) s,F(,),1/f(F(,)),()
其中s 称为稀疏函数(Sparsity function)或分位数密度函数(quantile density function)。s是分位数(τ)(τ)函数的导数,或在第τ分位数条件下概率密度函数的倒数(见Welsh,1988)。另外,模型误差项独立同分布假设意味着s与解释变量X无关,因此,分位数方程只和X在局部期间相关,即所有(τ)
2,1的条件分位数平面互相平行。事实上,式(15.5)中的就是误差项独立同分布假设(,(1,,)sJ),()
2下解释变量的回归系数估计量的渐近方差协方差矩阵表达式,而代表的是一般回归方,(1,,)s,()
程中随机误差项的方差。
误差项独立同分布假设下,分位数回归参数估计量的渐近方差协方差矩阵表达式中含有s(τ),但s(τ) 是未知分布的函数,而且必须要估计。
EViews提供了三种估计s的方法。两种是基于Siddiqui(1960)的方法分别提出的差分商方法(τ)
(Siddiqui Difference Quotient)(Koenker(1994)以及Bassett和Koenker(1982)),一种是核密度(Kernel Density)估计法。简述如下:
?Siddiqui差分商法:
差分商方法是用实际的分位数函数构造一个简单的差分商,从而求得s(τ)的估计量,表达式如下:
,1,1??F(h)F(h),,,,,nn?s() (15.8) ,,2hn
其中带宽h随着样本容量n??而趋向于0。要计算?(τ)需要做两件事,一是得到分位数函数n
,1?在两个点上的值,二是确定带宽。EViews中提供了两种Siddiqui 差分商法。 F(,)
计算分位数密度函数的第一种方法由Bassett和Koenker (1982)提出,EViews将其称之为Siddiqui (mean fitted) 方法。这种方法需要重新估计两个分位数回归模型在τ - h和τ + h上的拟nn和值,进而用不同的估计参数计算分位数函数的拟和值。最终s(τ)的估计量的数学表达式如下,
4
对任意X*有:
??(,h),(,h),,,,nn?, s(),X* (15.9) ,2hn
独立同分布假设意味着X*可以取任何值,Bassett和Koenker建议取X的均值,其优点是:估计的精度在该点达到最大;且估计的分位数函数对τ是单调的,因此对一个恰当的h,?(τ)的n值总是正的。
另一种Siddiqui 差分商法由Koenker(1994) 提出。其计算量相对较小,只需计算原分位数回归方程中残差的第τ - h和τ + h实际分位数,计算时排除在估计中设为零的k个残差,并插入新nn
值以获得分位数的分段线性形式。EViews中把这种方法叫做Siddiqui (residual) 方法。
上述两种Siddiqui方法都需要估计带宽h。EViews 提供了三种估计带宽的方法:Bofinger n
(1975) 法,Hall-Sheather (1988) 法和Chamberlain (1994)方法。
Bofinger(1975)提出的估计带宽的表达式为:
1/5,14,,,,,4.5((())),1/5,, (15.10) ,hTn,122,,[2(,(,)),1],,
可以近似最小化?(τ)的均方误差(MSE)。
另外两个带宽的表达式中含有显著性水平,因此常常用来进行假设检验。其中Hall和Sheather(1988)的表达式为:
1/3,12,,,,,1.5((())),1/32/3,, (15.11) ,hTZn,,12,,2(,(,)),1,,
-1 其中T表示样本容量,,表示正态分布的积累分布函数,,表示正态分布的密度函数,Z= ,(1-,/2),选择的显著性水平为,对应的Z值。
Chamberlain(1994)的表达式为:
,(1,,) (15.12) hZ,n,T
图1是样本容量1~300时Hall和Sheather(1988)方法在第0.1、0.3、0.5、0.7、0.9分位数下得到的带宽。图2是样本容量1~1000时三种方法在第0.5分位数下的带宽比较图 (α=0.05,MATLAB计算)。
图1 图2
从图2可以看出随着样本的增加,三种带宽都减小,并且在小样本时,减小的速度较大,在大样本情况下减小的速度较小。并且在大样本情况下,带宽的大小顺序为:Bofinger的最大,Hall和Sheather的次之,Chamberlain的最小。
?核密度法(Kernel Density):
5
,,1,1??根据(15.7)式有s= =1/f(),Falk(1988)和Welsh(1988)提出了用核密度法估计F(,)F(,)(τ)
,,1?进而得到s的方法。而Powell(1986)、Jones(1992)以及Buchinsky(1995)则通过估计F(,)(τ)
,1?1/f()来得到s。EViews中使用的方法属于后者,沿用了Powell(1984,1989)中的计算方法,F(,)(τ)
其选项名称为Kernel(residual):
T,1?? (15.13) s(,),1/[(1/T)cK(u/c)](),TiT,,1i
其中?表示分位数回归的残差;c为带宽;K表示核密度函数。EViews中可以选择的核密度函(τ)T
数有Epanechnikov核函数、均匀 (Uniform) 核函数、三角(Triangular)核函数、二权(Biweight)核函数、三权(Triweight)核函数、正态(Normal)核函数、余弦(Cosinus)核函数。
EViews中使用了Koenker(2005)提出的带宽,表达式为:
,1,1 (15.14) c,k(,(,,h),,(,,h))TTT
其中k表示Silverman(1986)的一个稳健估计量;h是Siddiqui带宽。 n
(2) 独立但不同分布假设下的参数渐近分布
?当分位数密度函数独立但不同分布即与解释变量X相关时的渐近分布服从T(,(,),,(,))Huber sandwich形式:
,1,1?N(0,,(1,,)H(,)JH(,))~ (15.15) T(,,,)(,)(,)
其中J同(15.6)式,H的表达式如下:
,H(,)lim(XXf(q(,))/T), (15.16) iiii,,,Ti
其中是个体i在第τ分位数上的条件密度函数。如果条件密度函数不依赖于观测值,式f(q(,))ii
(15.15)中的方差就退化为(15.5)式中的方差。
对于H,EViews提供了两种计算方法。第一种是Hendricks和Koenker(1992)提出的Siddiqui
差分法;另一种是Powell(1984,1989)提出的核密度法。这两种方法与在独立同分布假设时计算s(τ)的算法相同,因此在EViews选单中的名称相同,分别为Siddiqui (mean fitted)和Kernel
(residual)。
?Siddiqui差分商法
这种方法需要对每个个体估计τ - h和τ + h两个分位数回归模型,将拟和值代入下式: nn
,1,1???,,,f(q()),2h/(F(q(,h)),F(q(,h)))iiTiiTiiT (15.17) ??, ,2h/(X(,(,,h),,(,,h)))TiTT
由于分位数密度函数非同分布,因此,我们需要为每一个个体估计,这时当取f(q(,))ii
时,不能保证(15.17)式为正,因此,Hendricks和Koenker对其进行了修正: X,Xi
???,f(q(,)),max(0,2h/(X(,(,,h),,(,,h),,)) (15.18) iiTiTT
其中δ是一个很小的正数,避免上式中分母为零。将(15.18)式代入(15.16)式,得到H的估计量为
??,H(,),f(q(,))XX/T (15.19) iiii,
i
6
?核密度法
Powell(1984,1989)提出的用核密度法估计H的表达式为:
T,1??, (15.20) H(,),(1/T)cK(u/c)XX()T,iTii,,1i
其中?表示分位数回归的残差;c为带宽;K表示核密度函数;各参数含义与(15.13)式相(τ)n
同。
(3)参数渐近分布的自举法
前面的方法都是先求出分位数密度函数,然后再得到参数的渐近分布。自举法则可以省略这一步,直接得到参数的方差协方差阵。EViews中给出了四种自举方法,分别为:残差自举,XY对自举,以及两种马尔可夫链边际自举法MCMB和MBMB-A。其中前两种方法见Buchinsky (1995)。
?残差自举法(residual bootstrap)
这种方法要求解释变量与随机误差项不相关。它是对残差和解释变量分别进行有放回的再抽
*和X*(其中m可以小于原样本容量T),然后运用初始参数估样,构造样本容量为m的新序列u
***?计量构造被解释变量,即,最后用X*和Y*估计新的参数β(τ). ,y,X,u,()
如此重复K次,则参数方差协方差阵的估计量为:
Bm1????, (15.21) VT,(),()(,(,),,(,))(,(,),,(,))jj,TB,1j
其中是自举参数估计量序列的均值。EViews选单中称这种方法为Residual。 ,(,)
?XY对自举法(XY-pair or design bootstrap)
这是最常用的一种自举方法,它不要求随机误差项与解释变量相互独立。使用这种方法时,我们从原始数据中有放回的抽取K次样本容量为m的子序列(y*, X*),然后用每个子序列计算β(τ),最后运用(15.21)式计算参数方差协方差阵的估计量。EViews选单中这种方法称为XY-Pair。
?马尔可夫链边际自举法(Markov Chain Marginal Bootstrap)
以上两种自举法往往计算量过大,当方程中含有p个参数时,每次自举都需要解一个p维的线性规划问题。He和Hu(2002)提出了一种新的自举法,将一个p维的最优问题简化为求解一个含p个元素的序列的一维问题。这个序列的一维解就构成了一个马尔可夫链,其样本方差协方差阵可由(15.21)式计算,且当原序列样本容量T和自举次数K较大时具有一致性。EViews选单中把这种方法称为MCMB法。
然而,给定链长B(即自举次数),上述方法计算的参数序列之间往往存在较强的自相关从而导致参数方差协方差阵估计量的统计特性较差,有可能对任何链长B,估计量都不能收敛。Kocherginsky、He和Mu(KHM,2005)提出了一种修正的方法消除可能存在的自相关。即通过先对参数空间进行某种转换,运用MCMB算法进行估计,然后再转换回原来的空间,这种方法叫做MCMB-A。它要求独立同分布的假设条件,但它对异方差的情况表现的比较稳健。Kocherginsky、He和Mu还建议对于满足T,min(τ, 1-τ) > 5p的情况,当T ? 1000,p ? 10时,B应取在100至200之间。对于Tp在10,000到2,000,000之间的情况,建议B取在50至200之间,当然,还取决于用户的耐心。
15.5 分位数回归模型的检验
评价分位数回归函数好坏的统计量主要有3个,拟合优度、拟似然比检验和Wald检验。
(1)拟合优度(Goodness-of-Fit)
7
2Koenker和Machado(1999)提出了分位数回归的拟合优度的概念。它与一般回归分析中的R很类似。
假设分位数回归直线为
??, y,X,(,)(,)
???,,X,(1,Z)将解释变量矩阵和参数向量都分为两部分,即和,且有 ,,(,,,)(,)0(,)1(,)
?, y,,,Z,(,)0(,)1(,)
定义:
TT?????,, (15.22) Qmin[(1,)(y,Z,),(y,Z,)],,,,,,,,,,()0()101tt,,,,,()()()??,,,,,,::tyXtyXt,()t(,)
TT~?? (15.23) Q,min[,(1,,)(y,,),,(y,,)],()00tt,,,,()()??,,:,,,():,,(,)tyXtyXtt
式(15.22)和(15.23)分别表示无约束分位数回归目标函数(最小绝对离差和)和约束的分位数回归目标函数(最小绝对离差和)的极小值。无约束目标函数中的减项既包含常数项也包含所有回归因子。约束目标函数中的减项仅包含常数项,其他参数都约束为零。则Koenker和Machado拟和优度准则表达式如下:
?Q,()* (15.24) R,,1,()~Q(,)
~2*?很明显,上述统计量与传统的R非常相似。因为,所以R(τ)的值在0和1之间,解Q,Q(,)(,)
~**?R(,)R(,)释变量的作用越强,越远远小于,越接近1。反之越接近0。所以可用来QQ(,)(,)
考察解释变量对被解释变量第τ分位数回归拟和的好坏。
(2)拟似然比检验(Quasi-Likelihood Ratio Tests)
Koenker和Machado(1999)根据目标函数在施加约束条件前后得到的两个极小值构造了两个拟似然比检验统计量(QLR)。这两个拟似然比检验也称作分位数ρ检验(quantile-ρ tests)。两统计量的表达式如下:
~?2(Q,Q)(,)(,),L(), (15.25) T,(1,,)s(,)
~?QQ2,,()(),,, (15.26) ()log()T?,,,s,(1)()Q(,)
~2两个统计量都渐近服从自由度为q的χ分布,其中q是原假设目标函数中约束条件的个数。Q(,)
?和分别代表约束的和无约束目标方程的极小值。β0 Q(,)
另外,两个统计量的分母都含有稀疏项s(τ),上面给出的稀疏项s(τ)的3种计算方法都可在式(15.25)和(15.26)中使用。EViews估计的是其在备择假设下的估计量。
使用上述两统计量的前提是必须满足分位数密度函数s(τ)与解释变量X不相关。然而,尽管
8
有时并不满足独立同分布的假设,EViews在进行分位数回归的时候,不管选择何种估计参数渐近分布的方法,总会估计稀疏函数s(τ),从而构造拟似然比(QLR)检验统计量。因此,这种检验方法与下面的Wald统计量相比稳健性较差。
(3)Wald检验
给定分位数回归参数估计量的渐近方差协方差矩阵,我们就可以构造Wald形式的统计量进行各种约束形式的参数检验。
31.3.5 系列分位数回归检验
前面的分析主要集中在单个分位数回归模型的假设检验上,而有些时候也需要对一系列分位数回归的回归系数进行联合检验。比如,需要通过检验不同分位数模型的斜率是否相等来判断一个模型是否具有位移特征。同时考虑多个分位数回归式称作系列分位数回归分析(quantile process
testing)。EViews在做单方程分位数回归的同时,有专门命令执行系列分位数回归分析。
操作路径是在一个分位数回归估计结果窗口,点击View键,选Quantile Process/Process
Coefficients功能。
定义系列分位数回归系数列向量为,
, (15.27) ,(,',,',?,'),(,)(,)(,)12m
则有
?N(0,Ω)n(,,,) ~ (15.28)
其中Ω由如下形式的块矩阵Ω组成: ij(km×km)
,1,1 (15.29) Ω,[min(,,,),,,]H(,)JH(,)ijijijij
i, j=1, 2, … m. k为方程待估参数个数。其中J的表达式见(4)式。H的表达式见(15.19)或(15.20)式,取决于选择的估计方法。特别的,当误差项独立同分布的假设成立时,Ω简化为:
Ω,Ω,J (15.30) 0
其中Ω中的元素如下: 0
,,,,min(,),ijij,, (15.31) ij,1,1f(F(,))f(F(,))ij
i, j=1, 2, … k.除了以上的方法以外,Ω的估计量还可以由任何一种自举方法得到。
(1)斜率相等检验
Koenker和Bassett(1982a)提出了一种对异方差很稳健的判断不同分位数回归方程斜率是否相等的检验。零假设如下:
H:,,,,?,, 0111(,)(,)(,)12m
其中β指常数项以外的解释变量所对应的(k-1)维参数列向量。因此,零假设共含有(k-1) (m-1)个1
约束条件。接下来构造Wald形式的统计量检验零假设是否成立,它渐近服从自由度为(k-1) (m -1)
2的χ分布。
(2)对称性检验
将Newey和Powell(1987)检验最小二乘估计量对称性的方法扩展到分位数回归中。假设我们要检验的分位数回归模型有m个,m是奇数,且中间值τ是0.5,其他τ都关于0.5对称, (m+1)/2
即τ=1? τ, j=1,…,(m-1)/2。参数估计量按照τ的大小排序。则对称性检验的零假设为: jm-j+1k
9
,,,()(),,jm,j,1 (15.32) :H,,0(0.5)2
其中j=1, …, (m?1)/2。m是奇数,代表分位数回归个数。即关于0.5对称的分位数回归参数估计量的两两平均值等于中位数回归参数估计量。
我们可以构造Wald形式的统计量检验上述k(m-1)/2个约束条件是否成立。该统计量服从自
2由度为k(m?1)/2的χ分布。另外,Newey和Powell指出,如果我们已知随机误差项服从独立同分布,但不一定对称,则我们只需检验常数项的对称性。即
,,,00(),(),m,j,1j (15.33) :H,,00(0.5)2
这时约束条件减少为(m-1)/2个。θ
.6分位数的计算与分位数回归的EViews操作 15
(1)分位数的计算
对一个离散的随机变量y,取其容量为T的样本序列,计算第τ分位数的方法如下: t
首先将数据从小到大排序,标号为i,i =1, 2, …, T。然后利用下表所列的方法计算随机变量y的第τ分位数的排列序号的i;如果i为整数,则随机变量y的第τ分位数即为y,如果i不是tti整数,则随机变量y的第τ分位数为: t
y= y+ (i ? [i])( y? y) (τ)[i] [i]+1 [i]
其中[i]表示不大于i的最大整数。给定一个具体的随机变量y,对于一个容量为T的样本,t
则y的第τ分位数的序号i的计算方法如下。在大样本情况下,各方法收敛到同一值。 t
Rankit Ordinary Vander Waerden Blom Tukey Gumbel
(τ?1/2)/T τ/T τ/(T+1) (τ?3/8)/ (T+1/4) (τ?1/3)/ (T+1/3) (τ?1)/ (T?1)
计算分位数的EViews 6.0的命令为:scalar q=@quantile(y, τ, s),其中y表示求分位数的序列;τ表示要取的分位数;s取1~6依次表示上表中6种计算方法,计算所得结果存入标量q中。
例:打开6garch-03文件,在空白处键入命令:
scalar q=@quantile(DASH, .5,1)
scalar q=@quantile(DASH, .25,1)
意即对序列DASH求中位数。得结果DASH = -0.78,DASH序列的中位数是-0.78。DASH = (0.5)(0.25)-13.33,DASH序列的第0.25分位数是-13.33。
用DASH画分位数图如下。打开DASH序列窗口,点击View键选Graph功能。在打开的t
Graph Option窗口,Type选择页的Specifi选择框选Distribution,在Details的Distribution选择框中选Emprical Quantile如图。点击“确定”键,得分位数图如图。
10
160DASH120
80
40
0Quantile
-40
-80
-1200.00.20.40.60.81.0
2)分位数回归 (
主要包括3部分内容。(1)介绍怎样进行分位数回归。(2)对输出结果的分析。(3)对分位数回归相关功能键的介绍。
在EViews中进行分位数回归的路径有两个,分别是
(1)点击主选单中的Quick键,选Equation Estimation,弹出Equation Estimation窗口。 或者
(2)点击主选单中的Object键,选New Object,Equation,弹出Equation Estimation窗口。
在该窗口的Method下拉选单中,选择如图所示的选项QREG-Quantile Regression(including
LAD),EViews将打开如图所示的分位数回归对话框(Equation Estimation)。
图1
Equation Estimation(方程估计)窗口包括两个选项模块,一个是Specification(设定方程),一个是Options(选项)。
可以在Equation specification(方程设定)框中输入要估计的表达式。同一般线性回归模型一样,它可以是一行用空格隔开的被解释变量和解释变量(如图1所示),也可以是一个明确的参数为线性的表达式。
Equation Estimation(方程估计)窗口与OLS估计的Equation Estimation(方程估计)窗口相比,只多了对话框quantile to estimate的选项。在该处填入要估计的分位数。系统默认为0.5,即做中位数回归(LAD)。用户可以选择任意一个0和1之间的数(当数值接近0和1时估计会变得困难)。
激活Options(选项)模块(点击对话框上的Options(选项))。得到如图2的quantile regression
Options(分位数回归选择)选择框、Iteration control(迭代控制)选择框和Bootstrap settings(自举设定)选择框。
11
quantile regression Options对话框中的选择主要包括三部分。
图2
(1)Coefficient covariance(系数估计量方差协方差矩阵)选项框
其下拉选单中包括三个选项:Ordinary (IID),Huber-Sandwich 和Bootstrap,代表了可选的估计回归系数估计量方差协方差矩阵的方法(具体介绍见15.4节)。EViews 默认的是Huber-Sandwich方法。
(2)Weight(权数)选项框
可以输入作为权重的序列或者一个序列的表达式,从而对估计式加权。(用于WLS估计)
(3)Sparsity Estimation(稀疏函数估计)选项区
其中包括5种选择框。稀疏函数的介绍见15.4节。
? Method(方法)选项框。
当第一个选项框Coefficient covariance中选项为Ordinary (IID)或Bootstrap时,Method(方法)选项框中包括三个选项:Siddiqui (mean fitted), Kernel (residual)和Siddiqui (residual)。
当Coefficient covariance选项框中选项为Huber-Sandwich时,这里的Method选项框中只包括Siddiqui (mean fitted)和 Kernel (residual)两种选择。
? Bandwidth Method(带宽)选项框。
其下拉选单中包括四个选择,即Bofinger (1975),Hall-Sheather (1988)和Chamberlain (1994)
计算带宽方法,或者你自己给出一个特定的带宽。
? Size(置信尺度)选项框。
当选择Hall-Sheather 和Chamberlain方法时,置信度的选择默认为0.05。
? Quantile Method(分位数方法)选项框。
EViews提供了六种求解经验分位数的方法。
? Kernel(核函数)选项框。
表示核函数的选用种类。EViews中可以选择的核密度函数有Epanechnikov核函数、均匀核函数(Uniform)、三角核函数(Triangular)、二权核函数(Biweight)、三权核函数(Triweight)、正态核函数(Normal)和余弦核函数(Cosinus)。
注意,不管系数方差协方差矩阵(Coefficient covariance)是否会用到,每次进行分位数回归时,系统都会自动给出一个稀疏函数估计值。
12
Iteration control(迭代控制)选项框包括3个选项。
(1)Max(最大)。迭代的最大次数,默认为500。
(2)Starting(初始值)。表示迭代的初始值,默认为0,也可以选择其他选项,如下拉选单中的OLS,即用OLS估计量作为初始值进行迭代。
(3)Display settings(设定显示)。选择是否需要在输出结果中给出这些设置。
Bootstrap settings(自举设定)
(1)Method(方法)。代表不同的自举方法。EViews提供了四种方法,分别是Residual, XY-pair,
MCMB, MCMB-A。默认方法为XY-pair方法。
(2)Replications(循环次数)。EViews 默认为100次。用户可以自己设定次数。
(3)No. of obs(自举样本容量)。空白表示与原样本容量一致。Koenker(2005)的研究表明,选择自举样本容量小于数据样本容量时,能够获得更加准确的结果,特别是当数据样本容量较大时。
(4)output(输出)。在这里键入一个名称可以得到自举的样本矩阵。
(5)Random generator(生成随机数)和seed(种子)。本选项用于控制产生随机数。其中前者用于选择随机数产生方法,seed用于选择随机数种子,Clear(清除)按钮用于清空以往选定的随机数种子。
估计结果。
按照EViews默认设置得到的一个分位数回归估计结果如下:
输出结果上部给出的是估计设定,其中包括(按顺序)被解释变量(DASH)、方法:分位数回归(中位数)、操作日期、样本范围、样本容量(421)、标准误差和方差协方差矩阵估计方法 (Huber-Sandwich方法)、稀疏函数的估计方法(Kernel方法)、带宽方法(Hall-Sheather方法,带宽=0.12963)以及对估计结果的评价。
输出结果中部给出的是回归系数估计量、标准差、t统计量及其相应p值,这与OLS估计完全一样。可以看出,上述回归系数估计量都具有显著性。在中位数回归关系条件下,B股收益DBSH每增加一个单位,A股收益DASH平均增加3.38个单位。
输出结果下部给出的是对分位数回归估计式的评价统计量。分别为
2Pseudo R-squared:伪拟合优度(伪R),
Adjusted R-squared:调整的伪拟合优度,
S.E. of regression:分位数回归式的标准误差,
Quantile dependent var:分位数回归式中只有常数项存在的系数估计值(也即被解释变量的
13
分位数估计值)。
Objective:目标函数极小值,
Objective (const. only):分位数回归式中只有常数存在的目标函数极小值,
Sparsity:分位数密度函数(稀疏函数)估计值(本例是用核估计法计算的)。
Quasi-LR statistic:准似然比估计量的值
Prob (Quasi-LR stat):准似然比估计量的值所对应的概率值。
此外,由于这里使用的是Huber-Sandwich方法,因此稀疏函数值(Sparsity)并没有用来计算参数估计量标准差。
与上述结果类似,我们也可以通过改变估计设定,运用自举方法获得参数估计量的方差协方差矩阵。例如选择MCMB-A方法进行自举,并且将自举次数增加至500。对于稀疏函数的计算
点击OK键,得到新设定所对应的估计结果。 方法,选择Siddiqui(mean fitted),
分位数回归中的Views和Procs功能键。
分位数回归方程窗口中的大部分Views和Procs功能都与OLS回归相同,下面对一些计算细节其进行必要的补充说明。
使用上述功能时需要注意以下计算细节:
(1)这里的残差是指某一特定分位数回归函数条件下的残差,计算公式为
?,?Xβu,y -; (,)(,)tt
标准化残差指用自由度调整过的残差的标准误差。而在计算QLR统计量时则使用的是Koenker
1,??和Machado(1999)给出的目标函数极小值的平均值,即。 ,,,TQ()(,)
(2)构造Wald检验和置信椭圆时使用的是参数估计量方差协方差矩阵的稳健估计量。
(3)进行遗漏和多余变量检验(omitted and redundant tests)以及Ramsey RESET检验时,报告的都是特定约束下的QLR统计量,因此它只有在满足稀疏函数的独立同分布假设时才是有效的。
(4)在默认状态下,EViews在进行预测时会自动为样本外观测值扩展序列容量,用户可以在预测窗口中自行取消该选择框中的对勾。
Quantile Process选项
在分位数回归估计结果窗口中,点击View键,选Quantile Process功能,如图。该选项包括3个子功能,Process Coefficients(系列分位数回归),Slope Equality Test(斜率相等检验),Symmetric Quantiles Test(分位数对称检验)。
注意:由于每个功能都包含了不同分位数估计方程的比较和检验,因此运行时间会稍长。尤其是当选择自举方法计算方差协方差阵时,更是这样。
(1)Process Coefficients功能(系列分位数回归系数)
在分位数回归估计结果窗口中,点击View选Quantile Process,Process Coefficients(系列分
14
位数回归系数),可以得到Quantile Process对话框如图。
在该对话框中,可以通过Output选择区选择是按表格(Table)还是按图形(Graph)输出估计结果;选择置信度(Confidence Interval Size),EViews默认选择是输出表格形式(table)的估计结果,给出不同分位数下的系数估计量、标准差、t统计量及其相应的p值。
Quantiles后面的数字用于决定要估计的分位数个数。若用数字n表示,则输出的是第1/n,2/n,…,(n-1)/ n分位数回归估计结果。EViews默认数为10,即估计第0.1、0.2、…、0.9九个分位数回归式。通过选择分位数回归数(Quantiles)以及用图或列表形式(Output),对不同分位数回归系数估计量及其标准差进行比较与检验。
还可以选择User-specified quantiles并在下面的框中输入想计算的某个(些)分位数条件下的回归。
选取DASH对常数项和DBSH进行中位数回归,并在Quantile Process选项中选择输出95%
的折线图,如下所示: 置信水平的10个不同的分位数回归中回归系数估计值
Quantile Process Estimates (95% CI)
305.0DBSHC204.5
104.00
3.5-10
3.0-20
QuantileQuantile-302.50.00.20.40.60.81.00.00.20.40.60.81.0
点击Quantile Process对话框中的Output选项页,得到如下的对话框:
用Output选项页保存分位数向量,分位数回归系数估计量矩阵,分位数回归系数估计量的方差协方差矩阵。
15
在Quantiles Vector、Coefficient Matrix和Covariance Matrix选择框后面分别填入向量和矩阵名,EViews将把分位数序列、分位数回归系数估计量矩阵以及分位数回归系数估计量的方差协方差矩阵保存到相应向量和矩阵名中。
(2)Slope Equality Test(斜率相等检验)
在分位数回归估计结果窗口中点击View选Quantile Process,Slope equality test,可以进行Koenker和Bassett(1982a)提出的斜率参数相等检验,见如下对话框:
在该窗口的Specification(设定)选项页中,Quantiles后面的框中用户可以填入要检验的分位数个数,从而对上述分位数回归式的斜率估计量进行比较和检验。以填入数字n为例,表示输出的是第1/n,2/n,…,(n-1)/n分位数。
EViews默认的是4,即检验第1/4分位数回归式与第2/4分位数回归式之间、第2/4分位数回归式与第3/4分位数回归式之间相应回归系数是否相等。用户还可以选择User-specified
quantiles,并在下面的框中输入想要检验的分位数条件下的回归系数是否相等。
该窗口中的Output选项页的功能与Quantile Process对话框中的Output选项页的功能相同。
仍然使用上面的案例,按照系统默认的检验,得到如下的结果,即对中位数和第0.25分位数、以及中位数与第0.75分位数回归系数是否相等这2个约束条件进行检验。
输出结果中第一部分是对Wald检验的总结,这里斜率相等性检验的Wald统计量为0.53,自由度2,概率为0.7670。这表明0.25、0.5、0.75分位数回归式的斜率相等。这意味着,相同条件的不同分位数回归的被解释变量拟合值的分布是相同的。
(3)对称性检验
在分位数回归估计结果窗口中点击View选Quantile Process,Symmetric Quantiles Test,可以做以中位数回归式为中心的分位数回归式对称性检验。
16
, 依据Newey和Powell(1987)方法
,,,()(1-),, ,,(0.5)2
关于对称性检验的设置可以在Specification选项页中完成,其默认选项如图所示。如果Quantiles(分位数回归式个数)选择框中选4,则检验的约束条件只有一个,即
,,,(0.25)(0.75) ,,(0.5)2
而如果分位数回归式拟合的是0.6,则除了以上约束条件外,还会增加如下一个约束条件的检验。
,,,(0.4)(0.6) ,,(0.5)2
与上述两个检验类似,用户还可以选择User-specified quantiles,并在下面的框中输入要检验
-τ两个分位数回归式,并与中位数回归式做比较。按照的分位数τ,EViews将估计τ和1
(,,,)/2,,进行检验。 (,)(1-,)(0.5)
此外,还可以通过Intercept only或者All coefficients单选框选择是否仅对常数项还是对所有参数估计量进行对称性检验。同样,该窗口中的Output选项页的功能与Quantile Process对话框中的Output选项页的功能相同。
本例按照默认设置检验结果如下:
从输出结果中可以看出,这里我们检验的是第0.25和0.75分位数回归的系数估计量是否关于中位数对称。检验表明他们是对称的,概率p值是0.80。
17
预测功能与OLS单方程操作相同。
15.7 分位数回归的案例分析
案例分析:上证A股和B股收益关系研究(6garch-03)
上证A股、B股序列(ASH和BSH)和上证A股、B股收益序列(DASH和DBSH)如图。 tttt
1602,400ASHDASH20BSHDBSH1202,200
280802,00010
240401,800
020001,600
160-401,400-10120-80
80-20-1205010015020025030035040050100150200250300350400
图1 ASH和BSH序列 图2 DASH和DBSH序列 tttt
上证A股和B股收益序列(DASH和DBSH)的单位根检验结果如下: tt2DASH = -0.9786 DASHtt-1
(-20.0)* DW=1.99
2DBSH = -0.9599 DBSHtt-1
(-19.8)* DW=2.00
图3 上证A股收益序列(DASH)的单位根检验结果 t
图4 上证B股收益序列(DBSH)的单位根检验结果 t
上证A股和B股的差分(收益)序列(DASH和DBSH)都是平稳序列。 tt
下面分析DASH和DBSH的分布特征。打开DASH序列窗口,点击View键选Graph功能。ttt
在图形类型选择框中选Distribution(分布)。同时在选择窗右侧点击Options键,在弹出的对话窗中点击Add键,在弹出的Add对话窗中选Kernel Density。点击OK键。得直方图与核密度曲线合并图如图。分布都是单峰的,呈高峰厚尾特征。
18
.025.20HistogramHistogramDASHDBSHKernelKernel.020.16
.015.12
DensityDensity.010.08
.005.04
.000.00-150-100-50050100150200-20-16-12-8-40481216
分别求DASH和DBSH的第0.25、0.5、0.75分位数值(scalar q=@quantile(DASH, τ, 1)),τ=0.25、tt
0.5、0.75,得下表。
均值 偏度 标准差 第0.25分位数 第0.5分位数 第0.75分位数
DASH -1.40 0.77 25.7 -13.33 -0.78 9.74 t
DBSH 0.04 -0.05 4.1 -1.510 -0.06 1.73 t
上证A股和B股收益序列散点图如下:
160
120
80
40
DASH0
-40
-80
-120-20-15-10-5051015
DBSH
OLS估计结果如下:B股收益每增加一个单位,A股收益增加3.93个单位。
19
DASH对DBSH的中位数回归结果如下。 tt
160DASHDASHF120DASHF50
80
40
0
-40
-80DBSH-120-20-15-10-5051015
OLS回归、LAD(中位数)回归直线对比如图。OLS回归直线的斜率大于LAD(中位数)回归直线的斜率,说明OLS估计方法对离群值敏感,LAD估计方法比OLS估计方法稳健。
在分位数回归输出结果窗口点击View键,选Quantile Process,Process Coefficients,在弹出的对话框中点击确定键,得9个分位数回归式如下,DBSH的回归系数全部具有显著性。
20
9个分位数回归式中的回归系数值折线图如下,回归系数值在3.38和3.97之间变化。除了第0.1分位数回归式之外,其他8个分位数回归式的回归系数基本相同,说明不同分位数回归式之间,DASH和DBSH的回归关系基本相同。
Quantile Process Estimates (95% CI)
305.0DBSHC204.5
104.00
3.5-10
3.0-20
QuantileQuantile-302.50.00.20.40.60.81.00.00.20.40.60.81.0
OLS回归直线和第0.25、0.5、0.75分位数回归直线如下图。回归直线的斜率大于中位数回归线的斜率。
160DASHDASHF120DASHF25DASHF50DASHF7580
40
0
-40
-80
DBSH-120-20-15-10-5051015
斜率相等检验:
在分位数回归输出结果窗口点击View键,选Quantile Process,Slope equality test功能,在弹出的对话框中点击确定键,得输出结果如下。经检验,第0.25和0.5分位数回归直线的斜率,第0.5和0.75分位数回归直线的斜率之间无显著性差异。
原假设有两个H:β= β,β = β 01(0.25)1(0.5)1(0.5)1(0.75)
21
对称性检验
+ β = 2β ,β+ β = 2β 原假设有两个H:β00(0.25)0(0.75)0(0.5)1(0.25)1(0.75)1(0.5)
经检验,第0.25和0.75分位数回归直线的斜率以中位数斜率存在对称性。
案例分析:恩格尔定律研究。(分位数回归,) (file:6food,Murray-book 第1章)
、0.75分位数值(scalar q=@quantile(food, τ, 1)),τ=0.25、求food的分布直方图,第0.25、0.5
0.5、0.75,得下表。
25Series: FOODSample 1 199Observations 19920
Mean 580.0058Median 558.520015Maximum 1881.880Minimum 223.7500Std. Dev. 238.826810Skewness 1.568316Kurtosis 7.646177
5Jarque-Bera 260.5692Probability 0.000000
020040060080010001200140016001800
均值 标准差 偏度 0.1分位数 0.25分位数 0.5分位数 0.75分位数 0.9分位数
580.0 238.8 1.57 324.4 410.0 558.5 696.5 856.2 food
22
2,000.0025FOODHistogramKernel1,600.0020
1,200.0015
FOOD800Density.0010
400
.0005
004008001,2001,6002,0002,4002,8003,200.000002004006008001,0001,2001,4001,6001,8002,0002,200TOTAL
点击主选单中的Quick键,选Equation Estimation,弹出Equation Estimation窗口。在该窗口的Method下拉选单中,选择QREG-Quantile Regression(including LAD),EViews将打开分位数回归对话框(Equation Estimation)。点击确定键,中位数回归估计结果:
,?? = = 75.6219 + 0.5419 total ,,,totalfood0(0.5)1(0.5)t
(3.3) (18.9)
在分位数回归输出结果窗口点击View键,选Quantile Process,Process Coefficients,在弹出的对话框中点击确定键,得9个分位数回归式如下,DBSH的回归系数全部具有显著性。
23
第10、30、50、70、90分位数回归与均值回归拟合直线的比较
2,000FOODFOODFFOODF101,600FOODF30FOODF50FOODF70FOODF901,200
800
400
005001,0001,5002,0002,5003,000
TOTAL
系列分位数回归系数(Quantile Process Estimates)估计值及其置信区间(95%)。
.8150TOTALC125.7
100.6
75.550
.425
QuantileQuantile.300.00.20.40.60.81.00.00.20.40.60.81.0
用第10、20、30、40、50、60、70、80、90分位数回归得到的回归系数估计值各自连成的线(蓝线)以及相应的95%置信度的置信区间(红线)。估计结果显示,分位数回归式随着分位数的增加,食品支出与总支出的回归系数也增加,数据存在异方差。
分位数系数相等性检验结果如下。第0.25与0.5分位数回归式之间,第0.5与0.75分位数回归式之间的斜率不相等。这也从一个角度显示数据存在递增型异方差。
原假设有两个H:β= β,β = β 01(0.25)1(0.5)1(0.5)1(0.75)
24
分位数系数对称性检验如下。第0.25与0.75分位数回归式的斜率对称。第0.25与0.75分位数回归式对于中位数回归式来说,具有对称性。
原假设有两个H:β+ β = 2β ,β+ β = 2β 00(0.25)0(0.75)0(0.5)1(0.25)1(0.75)1(0.5)
下面研究对数数据情形。
中位数回归估计结果:
,??= = 0.4991 + 0.8595 Lntotal ,,,LntotalLnfood0(0.5)1(0.5)t
(1.9) (21.7)
系列分位数(第20、40、60、80、90分位数)回归结果
25
、40、60、80、90分位数回归与均值回归拟合直线的比较 第20
7.6LNFOODLNFOODFLNFOODF107.2LNFOODF30LNFOODF50LNFOODF706.8LNFOODF90
6.4
6.0
5.6
LNTOTAL5.25.66.06.46.87.27.68.0
系列分位数回归系数(Quantile Process Estimates)估计值及其置信区间(95%)。
1.6.96LNTOTALC
1.2.92
0.8.88
0.4.84
0.0.80
QuantileQuantile-0.4.760.00.20.40.60.81.00.00.20.40.60.81.0
用第10、20、30、40、50、60、70、80、90分位数回归得到的回归系数估计值各自连成的
线(蓝线)以及相应的95%置信度的置信区间(红线)。
分位数系数相等性检验。经检验,4个分位数回归式的斜率相等。
26
分位数系数对称性检验如下。经检验,两对分位数回归式的斜率对称。
27
范文三:分位数回归。
分位数回归
三部分:分位数回归简介
分位数回归的应用
R 程序实践
一、分位数回归简介
为什么要分位数回归?
传统的线性回归描述条件均值受自变量的影响,若随机误差满足经典假设,参数估计将具有无偏性、有效性等优良性质。但实际生活假设往往不满足,如存在异方差,偏态分布等会使传统线性回归不具有以上性质。
分位数回归1、随机扰动项不做分布的假定,估计具有很强稳健型
2、对所有分位数进行回归,这样对异常点具有抗耐性 一体两面的,更加精确地描述自变量对因变量变化范围的影响 3、分位数回归具有较好的弹性性质
4、对于因变量具有单调变换性
5、估计参数在大样本下具有渐进优良性
为了方便解释清楚分位数回归, 先利用一个图形来作简要说明:
上图的横坐标表示的是家庭收入,而纵坐标表示的是食物支出。这个例子稍后会用R 实现。
分位数回归原理
回归分析的基本思想就是使样本值与拟合值之间的距离最短,对于Y 的一组随机样本
,样本均值回归是使误差平方和最小,即
样本中位数回归是使误差绝对值之和最小,即
样本分位数回归是使加权误差绝对值之和最小,即
现假设因变量Y 由k 个自变量组成的矩阵X 线性表示,对于条件均值函数得到参数
其中
加权表述方式。
二、分位数回归的应用 为检查函数,等价于上述
近10多年来, 分位数回归的理论和方法在各个领域中都得到了非常迅速的发展:
在环境科学方面, 典型的有Chock,Winkler 和Chen 使用非参数分位数回归法研究了匹兹堡这座城市中日死亡率和空气污染集中度的相互关系;
在生存分析方面,Koenker 和Hallock(2001)研究了诸多因数对于新生儿出体重的影响。
Cole 和Green 以及Royston 和Altman 讨论了分位数回归在医学上的应用。
Deaton 对于分位数回归在需求分析方面上的应用做了介绍, 并分析了巴基斯坦的Engel 曲线,等等
收入不平等问题是分位数回归的另一个研究方面,Gosling 、Machin 和Meghir 研究了英国家庭的收入和财富的分布状况;
分位数回归目前在金融方面的应用主要在两个方面。金融资产组合方面:Bassett和Chen 运用分位数回归来评估共同基金的投资类型金融风险方面,Taylor 使用分位数回归的方法来估计多期收益的风险值VaR 。
三、R 程序实践
1、了解包quantreg :包括文档、代码;
文档vignette 其他文档有crq 、 rqss 文档
2、了解线性分位数函数rq ()
rq(formula, tau=.5, data, subset, weights, na.action,
method="br", model = TRUE, contrasts, ...)
method=: 此参数指定用于计算分位数回归的算法
1、默认为“br ”
2、参数设置为“fn ”
3、参数设置为“fnc ”
介绍返回值:参数tau 决定返回值得对象类型不同。参数tau 在01之间返回rq 对象;tau 在01之外返回rq.process 对象
3、提取拟合模型信息的方式
第一、根据对象类型使用函数
对于rq 对象,回归系数、残差、效应这些特定模型信息利用对应的通用函数提取,coef(rq)、resid(rq)等,而不是$号提取,这点需要注意。
对于Rq.process 对象,包括sol 和dsol 这两种结构,并且必须使用$提取这两种结构成分,成分包括分位数点、参数系数等等
第二、利用汇总函数summary() 以及其参数提取信息 信息包括估计参数的置信区间和其显著性检验。
summary(object, se = NULL, covariance=FALSE, hs = TRUE, ...)
se:这个参数可以设置为“rank ”也是默认设置。假设残差独立同分布。提供估计参数及其置信区间。
Se=“iid ”假设残差独立同分布。用KB 的方法提供参数估计及显著t 检验 Se=“nid ”按照Huber 方法得到估计及检验
Se=“ker ”采用Powell 的核估计方法
Se=“boot ”采用自助法得到估计及检验。
covariance=FALSE:是否返回估计系数的协方差矩阵,默认不返回。 hs = TRUE:关于带宽的估计方法。
分位数回归程序举例:三件事
一、作单分位点一般为中位数 多分位点的分位数回归拟合并分别提取信息 ;并制图绘出各分位点的估计参数值,形象说明解释变量是怎样影响各个分位点因变量
二、做分位数回归的拟合直线 :为直观对比 我们绘制了变量散点图;条件均值回归线; 中位数回归线;各分位数回归线
三、针对缺失值的情况
范文四:3-分位数回归
第15章 分位数回归模型
15.1 总体分位数和总体中位数 15.2 总体中位数的估计 15.3 分位数回归
15.4 分位数回归模型的估计 15.5 分位数回归模型的检验
15.6 分位数的计算与分位数回归的EViews 操作 15.7 分位数回归的案例分析
以往介绍的回归模型实际上是研究被解释变量的条件期望。人们当然也关心解释变量与被解释变量分布的中位数,分位数呈何种关系。这就是分位数回归,它最早由Koenker 和Bassett(1978)提出,是估计一组回归变量X 与被解释变量Y 的分位数之间线性关系的建模方法。
正如普通最小二乘OLS 回归估计量的计算是基于最小化残差平方和一样,分位数回归估计量的计算也是基于一种非对称形式的绝对值残差最小化,其中,中位数回归运用的是最小绝对值离差估计(LAD,least absolute deviations estimator)。它和OLS 主要区别在于回归系数的估计方法和其渐近分布的估计。在残差检验、回归系数检验、模型设定、预测等方面则基本相同。 分位数回归的优点是,(1)能够更加全面的描述被解释变量条件分布的全貌,而不是仅仅分析被解释变量的条件期望(均值),也可以分析解释变量如何影响被解释变量的中位数、分位数等。不同分位数下的回归系数估计量常常不同,即解释变量对不同水平被解释变量的影响不同。 另外,中位数回归的估计方法与最小二乘法相比,估计结果对离群值则表现的更加稳健,而且,分位数回归对误差项并不要求很强的假设条件,因此对于非正态分布而言,分位数回归系数估计量则更加稳健。
15.1 总体分位数和总体中位数
在介绍分位数回归之前先介绍分位数和中位数概念。
对于一个连续随机变量y ,其总体第τ分位数是y (τ) 的定义是:y 小于等于y (τ) 的概率是τ,即 τ = P( y ≤ y (τ) ) = F (y (τ) )
其中P(?) 表示概率,F (y (τ) ) 表示y 的累积(概率)分布函数(cdf)。 比如y (0.25) = 3,则意味着y ≤ 3的概率是0.25。且有 y (τ) = F -1(y (τ) )
即F (y (τ) ) 的反函数是y (τ) 。当τ=0.5时,y (τ) 是y 的中位数。τ= 0.75时,y (τ) 是y 的第3/4分位数,τ= 0.25时,y (τ) 是y 的第1/4分位数。若y 服从标准正态分布,y (0.5) = 0,y (0.95) =1.645,y (0.975) =1.960。
另外,如果随机变量y 的分布是对称的,那么其均值与中位数是相同的。当其中位数小于均值时,分布是右偏的。反之,分布是左偏的。
对于回归模型,被解释变量y t 对以X 为条件的第τ分位数用函数y (τ) t |X 表示,其含义是:以X 为条件的y t 小于等于y (τ) t |X 的概率是τ。这里的概率是用y t 对X 的条件分布计算的。且有 y (τ) t |X = F -1(y (τ) t |X )
其中F (y (τ) t |X ) 是y t 在给定X 条件下的累积概率分布函数(cdf)。则y (τ) t |X 称作被解释变量y t 对X 的条件分位数函数。而F ' (y (τ) t |X )= f (y (τ) t |X ) 则称作分位数概率密度函数。其中F ' (y (τ) t |X ) 表示F (y (τ) t |X )
对y (τ) t |X 求导。
15.2 总体中位数的估计
在介绍分位数回归之前,先来看中位数的估计和中位数回归。下面以连续变量为例介绍定理15.1。
定理15.1
连续变量用y 表示,其概率密度函数用f (y ) 表示,累计概率密度函数用F (y ) 表示,y 的中位数用y (0.5)表示,则y 与任一值α的离差绝对值的期望E (y -α) 以α = y (0.5) 时为最小。 证明: E (y -α) =-
?-∞
α
(y -α) f (y ) dy +(y -α) dF (y ) +
?α(y -α) f (y ) dy
?α(y -α) dF (y ) (15.1) ?a
b ∞
∞
=-
?-∞
α
根据莱布尼兹公式,若F (α) =有F '(α) =
f (y , α) dy ,则有F '(α) =
?a
b ?f (y , α)
?α
。令f (y , α) =y -α,则
?a
b ?(y -α)
?α
,得 =-dy 。运用于式(15.1)
?a
b
?E (y t -α)
?α
?
=-
?-∞
α
(y -α) f (y ) dy ?α
α
?+
?α(y -α) f (y ) dy =
?α
∞
?-∞
α
dF (y ) -
?αdF (y )
∞
= F (α) -[1-
?-∞dF (y ) ]=F (α) -(1-F (α)) =2F (α) -1
?E (y t -α)
?α
式(15.1)求极小的一阶条件是于中位数y (0.5)。
α = y (0.5)
与定理15.1等价的表述是
= 0,即2F (α) -1=0,F (α) =0.5。这意味着α等
∑y -α以α = y
(0.5)(中位数)时为最小。因此,中位数回归估计
量可以通过最小绝对离差法(least absolute deviation, LAD)估计。其中X 和β分别为(k ?1) 阶列向量。
同理,对于线性回归模型y t = X 'β + u t ,通过求
?
(0.5)最小,估计β的中位数回归系∑y t -X 'β
???(0. 5) t X ) =X β数估计量β(0.5),从而得到y t 的中位数回归估计量(y (0.5)。
15.3 分位数回归
?(τ) t 表示y t 的分位数回归估计量,则对于以检查函数 Koenker 和Bassett(1978)证明,若用y
(check function)w τ为权数,y t 对任意值α的加权离差绝对值和得最小值。其中
?(τ) t 时取∑w τy t -α只有在α =y
(1-τ)(y t -α) +∑τ(y t -α) (15.2) ∑w τy t -α = -i :∑
y <αt :y="">αt>
i
i
T T
τ∈(0, 1)。据此,分位数回归可以通过加权的最小绝对离差和法(weighted least absolute deviation,
WLAD )进行估计。
?的 根据式(15.2),对于线性回归模型y t = X 'β + u t , 求第τ分位数回归方程系数的估计量β(τ)
方法是求下式(目标函数)最小, Q =-
?(τ) t <>
∑(1-τ) u ?(τ) t +∑τu ?(τ) t
?(τ) t ≥0u
T T
=-
t :y t
?) +?) (15.3) (y t -X 'β(τ) (τ) ∑?(1-τ)(y t -X 'β∑τ
?
t :y t ≥X 'β(τ)
T T
?(τ) t 表示第τ分位数回归方程对应的残差。τ∈(0, 1)。第τ分位数的回归方程表达式是 其中u
? ?(τ) t =X 'β y (τ)
?称作分位数回归系数估计量,或最小绝对离差和估计量,估其中X ,β都是k ?1阶列向量。β(τ)
计方法称作最小绝对离差和估计法。 当τ=0.5时,式(15.3)变为 Q =-
t :y t
∑?
T
?0.5(y t -X 'β(0.5) ) +
t :y t ≥X 'β(0.5)
∑?
T
?0.5(y t -X 'β(0.5) ) =
?
(0.5) ∑0.5y t -X 'β
t =1
T
???(0.5) t =X 'βy (0.5) 称作中位数回归方程,β(0.5) 称作中位数回归系数估计量。 ?(τ) t 。 一旦得到估计的分位数回归方程,就可以计算分位数回归的残差u
? ?(τ) t =y t -y ?(τ) t =y t -X 'β u (τ)
对一个样本,估计的分位数回归式越多,对被解释变量y t 条件分布的理解就越充分。以一元
回归为例,如果用LAD 法估计的中位数回归直线与用OLS 法估计的均值回归直线有显著差别,则表明被解释变量y t 的分布是非对称的。如果散点图上侧分位数回归直线之间与下侧分位数回归直线之间相比,相互比较接近,则说明被解释变量y t 的分布是左偏倚的。反之是右偏倚的。对于不同分位数回归函数如果回归系数的差异很大,说明在不同分位数上解释变量对被解释变量的影响是不同的。
15.4 分位数回归模型的估计
由于目标函数(15.3)不可微,因此传统的对目标函数求导的方法不再适用。估计分位数回?的一种较好的方法是线性规划方法。 归方程参数β(τ)
基于Barrodale 和Roberts (1973,以下简写为BR ) 提出的单纯形法(simplex algorithm),Koenker
和D’Orey(1987)提出一种估计分位数回归系数的方法。EViews 中应用的是上述算法的改进形式。
BR 算法由于其非有效性和大样本下的一些非优良特性曾备受批评。Koenker 和Hallock(2001) 以及Portnoy 和Koenker(1997)通过模拟证实,与内点法(interior point method) 等替代方法相比,BR 算法的估计次数往往较多,大约是样本容量的平方次数。然而,改进的BR 算法的估计次数在一定程度上是可以接受的,大约是样本容量的线性倍次数,在实际中是可以使用的。
分位数回归方程的BR 算法原理略。
下面讨论分位数回归系数估计量的渐近分布。
在弱条件下,分位数回归系数渐近服从正态分布(Koenker, 2005)。回归系数的方差协方差矩阵的计算在分位数回归的系数估计中占有重要位置。其方差协方差矩阵的估计方法根据分位数密度函数是否与解释变量相关分为三种方法:
①误差项独立同分布(i.i.d. ) 假设下的直接估计方法。由Koenker 和Bassett(1978)提出。 ②误差项独立但不同分布(i.n.i.d. ) 条件下的直接估计方法。 ③误差项独立同分布(i.i.d. ) 和独立但不同分布(i.n.i.d. ) 条件下都可使用的自举法。
(1)独立同分布假设下的参数渐近分布
Koenker 和Bassett(1978)在独立同分布假设下得出分位数回归系数渐近服从正态分布,可以表述为在弱条件下: 其中
J =lim (
n →∞
?-β) ~N (0, τ(1-τ) s 2J -1) (15.5) n (β(τ) (τ) (τ)
∑
i
X i X i 'X 'X
) =lim () (15.6) T n →∞T
s (τ) =F -1(τ) =1/f (F -1(τ)) (15.7)
其中s (τ) 称为稀疏函数(Sparsity function)或分位数密度函数(quantile density function)。s (τ) 是分位数
函数的导数,或在第τ分位数条件下概率密度函数的倒数(见Welsh,1988) 。另外,模型误差项独立同分布假设意味着s (τ) 与解释变量X 无关,因此,分位数方程只和X 在局部期间相关,即所有的条件分位数平面互相平行。事实上,式(15.5)中的(τ(1-τ) s (τ) 2J -1) 就是误差项独立同分布假设下解释变量的回归系数估计量的渐近方差协方差矩阵表达式,而τ(1-τ) s (τ) 2代表的是一般回归方程中随机误差项的方差。
误差项独立同分布假设下,分位数回归参数估计量的渐近方差协方差矩阵表达式中含有s (τ) ,但s (τ) 是未知分布的函数,而且必须要估计。
EViews 提供了三种估计s (τ) 的方法。两种是基于Siddiqui(1960)的方法分别提出的差分商方法(Siddiqui Difference Quotient)(Koenker(1994)以及Bassett 和Koenker(1982)),一种是核密度(Kernel Density )估计法。简述如下:
①Siddiqui 差分商法:
差分商方法是用实际的分位数函数构造一个简单的差分商,从而求得s (τ) 的估计量,表达式如下:
?-1(τ+h ) -F ?-1(τ-h ) F n n ?(τ) = s (15.8)
2h n
'
其中带宽h n 随着样本容量n →∞而趋向于0。要计算?(τ) 需要做两件事,一是得到分位数函数
?-1(τ) 在两个点上的值,二是确定带宽。EViews 中提供了两种Siddiqui 差分商法。 F
计算分位数密度函数的第一种方法由Bassett 和Koenker (1982)提出,EViews 将其称之为Siddiqui (mean fitted) 方法。这种方法需要重新估计两个分位数回归模型在τ - h n 和τ + h n 上的拟和值,进而用不同的估计参数计算分位数函数的拟和值。最终s (τ) 的估计量的数学表达式如下,
对任意X *有:
?(τ) =X *' s
?(τ+h ) -β?(τ-h ) βn n
2h n
(15.9)
独立同分布假设意味着X *可以取任何值,Bassett 和Koenker 建议取X 的均值,其优点是:
估计的精度在该点达到最大;且估计的分位数函数对τ是单调的,因此对一个恰当的h n ,?(τ) 的值总是正的。
另一种Siddiqui 差分商法由Koenker(1994) 提出。其计算量相对较小,只需计算原分位数回归方程中残差的第τ - h n 和τ + h n 实际分位数,计算时排除在估计中设为零的k 个残差,并插入新值以获得分位数的分段线性形式。EViews 中把这种方法叫做Siddiqui (residual) 方法。
上述两种Siddiqui 方法都需要估计带宽h n 。EViews 提供了三种估计带宽的方法:Bofinger (1975) 法,Hall-Sheather (1988) 法和Chamberlain (1994)方法。
Bofinger(1975)提出的估计带宽的表达式为:
-14?-1/5 4. 5(φ(Φ(τ))) h n =T
[2(Φ-1(τ)) 2+1]2?
????
1/5
(15.10)
可以近似最小化?(τ) 的均方误差(MSE)。
另外两个带宽的表达式中含有显著性水平,因此常常用来进行假设检验。其中Hall 和Sheather(1988)的表达式为:
-12?-1/32/3 1. 5(φ(Φ(τ))) Z α h n =T
2(Φ-1(τ)) 2+1?
????
1/3
(15.11)
其中T 表示样本容量,Φ表示正态分布的积累分布函数,φ表示正态分布的密度函数,Z α= Φ-1 (1-α/2)为选择的显著性水平α对应的Z 值。
Chamberlain(1994)的表达式为: h n =Z α
τ(1-τ)
T
(15.12)
图1是样本容量1~300时Hall 和Sheather(1988)方法在第0.1、0.3、0.5、0.7、0.9分位数下得到的带宽。图2是样本容量1~1000时三种方法在第0.5分位数下的带宽比较图 (α=0.05,MATLAB 计算)。
图1 图2
从图2可以看出随着样本的增加,三种带宽都减小,并且在小样本时,减小的速度较大,在大样本情况下减小的速度较小。并且在大样本情况下,带宽的大小顺序为:Bofinger 的最大,Hall 和Sheather 的次之,Chamberlain 的最小。
②核密度法(Kernel Density):
?-1'(τ) =1/f (F ?-1(τ) ) ,Falk(1988)和Welsh(1988)提出了用核密度法估计根据(15.7)式有s (τ) = F
?-1'(τ) 进而得到s (τ) 的方法。而Powell(1986)、Jones(1992)以及Buchinsky(1995)则通过估计F
?-1(τ) ) 来得到s (τ) 。1/f (F EViews 中使用的方法属于后者,沿用了Powell(1984,1989)中的计算方法,
其选项名称为Kernel(residual): ?(τ) =1/[(1/T ) s
∑c T -1K (u ?(τ) i /c T )] (15.13)
i =1
T
其中?(τ) 表示分位数回归的残差;c T 为带宽;K 表示核密度函数。EViews 中可以选择的核密度函
数有Epanechnikov 核函数、均匀 (Uniform) 核函数、三角(Triangular)核函数、二权(Biweight)核函数、三权(Triweight)核函数、正态(Normal)核函数、余弦(Cosinus)核函数。
EViews 中使用了Koenker(2005)提出的带宽,表达式为:
c T =k (Φ-1(τ+h T ) -Φ-1(τ-h T )) (15.14)
其中k 表示Silverman(1986)的一个稳健估计量;h n 是Siddiqui 带宽。
(2) 独立但不同分布假设下的参数渐近分布
?(τ) -β(τ)) 的渐近分布服从当分位数密度函数独立但不同分布即与解释变量X 相关时(β
Huber sandwich形式:
?-β) ~N (0, τ(1-τ) H (τ) -1JH (τ) -1) (15.15) (β(τ) (τ)
其中J 同(15.6)式,H 的表达式如下:
H (τ) =lim (
X i X i 'f i (q i (τ)) /T ) (15.16) ∑T →∞
i
其中f i (q i (τ)) 是个体i 在第τ分位数上的条件密度函数。如果条件密度函数不依赖于观测值,式(15.15)中的方差就退化为(15.5)式中的方差。
对于H ,EViews 提供了两种计算方法。第一种是Hendricks 和Koenker(1992)提出的Siddiqui 差分法;另一种是Powell(1984,1989)提出的核密度法。这两种方法与在独立同分布假设时计算s (τ) 的算法相同,因此在EViews 选单中的名称相同,分别为Siddiqui (mean fitted) 和Kernel (residual)。
①Siddiqui 差分商法
这种方法需要对每个个体估计τ - h n 和τ + h n 两个分位数回归模型,将拟和值代入下式:
?(q (τ)) =2h /(F ?-1(q (τ+h )) -F ?-1(q (τ-h ))) f i i T i i T i i T
(15.17)
?? =2h /(X '(β(τ+h ) -β(τ-h )))
T
i
T
T
由于分位数密度函数非同分布,因此,我们需要为每一个个体估计f i (q i (τ)) ,这时当取
X i =时,不能保证(15.17)式为正,因此,Hendricks 和Koenker 对其进行了修正:
?(q (τ)) =m ax(0, 2h /(X '(β?(τ+h ) -β?(τ-h ) -δ)) (15.18) f i i T i T T
其中δ是一个很小的正数,避免上式中分母为零。将(15.18)式代入(15.16)式,得到H 的估计量为
?(q (τ)) X X '/T (15.19) ?(τ) =f H
∑
i
i i i i
②核密度法
Powell(1984,1989)提出的用核密度法估计H 的表达式为: ?(τ) =(1/T ) H
∑c T -1K (u ?(τ) i /c T ) X i X i ' (15.20)
i =1
T
其中?(τ) 表示分位数回归的残差;c n 为带宽;K 表示核密度函数;各参数含义与(15.13)式相
同。
(3)参数渐近分布的自举法
前面的方法都是先求出分位数密度函数,然后再得到参数的渐近分布。自举法则可以省略这一步,直接得到参数的方差协方差阵。EViews 中给出了四种自举方法,分别为:残差自举,XY 对自举,以及两种马尔可夫链边际自举法MCMB 和MBMB-A 。其中前两种方法见Buchinsky (1995)。
①残差自举法(residual bootstrap)
这种方法要求解释变量与随机误差项不相关。它是对残差和解释变量分别进行有放回的再抽样,构造样本容量为m 的新序列u *和X *(其中m 可以小于原样本容量T ) ,然后运用初始参数估
?+u *,最后用X *和Y *估计新的参数β(τ). 计量构造被解释变量,即y *=X *β(τ)
如此重复K 次,则参数方差协方差阵的估计量为: ?) =T (m ) 1?(β V
T B
∑(β?j (τ) -β(τ) ) (β?j (τ) -β(τ) ) ' (15.21)
j =1
B
其中() 是自举参数估计量序列的均值。EViews 选单中称这种方法为Residual 。
②XY 对自举法(XY-pair or design bootstrap)
这是最常用的一种自举方法,它不要求随机误差项与解释变量相互独立。使用这种方法时,我们从原始数据中有放回的抽取K 次样本容量为m 的子序列(y *, X *),然后用每个子序列计算β(τ) ,最后运用(15.21)式计算参数方差协方差阵的估计量。EViews 选单中这种方法称为XY-Pair 。
③马尔可夫链边际自举法(Markov Chain Marginal Bootstrap) 以上两种自举法往往计算量过大,当方程中含有p 个参数时,每次自举都需要解一个p 维的线性规划问题。He 和Hu(2002)提出了一种新的自举法,将一个p 维的最优问题简化为求解一个含p 个元素的序列的一维问题。这个序列的一维解就构成了一个马尔可夫链,其样本方差协方差阵可由(15.21)式计算,且当原序列样本容量T 和自举次数K 较大时具有一致性。EViews 选单中把这种方法称为MCMB 法。
然而,给定链长B (即自举次数) ,上述方法计算的参数序列之间往往存在较强的自相关从而导致参数方差协方差阵估计量的统计特性较差,有可能对任何链长B ,估计量都不能收敛。Kocherginsky 、He 和Mu(KHM,2005)提出了一种修正的方法消除可能存在的自相关。即通过先对参数空间进行某种转换,运用MCMB 算法进行估计,然后再转换回原来的空间,这种方法叫做MCMB-A 。它要求独立同分布的假设条件,但它对异方差的情况表现的比较稳健。Kocherginsky 、He 和Mu 还建议对于满足T ?min(τ, 1-τ) > 5p 的情况,当T ≤ 1000,p ≤ 10时,B 应取在100至200之间。对于Tp 在10,000到2,000,000之间的情况,建议B 取在50至200之间,当然,还取决于用户的耐心。
15.5 分位数回归模型的检验
评价分位数回归函数好坏的统计量主要有3个,拟合优度、拟似然比检验和Wald 检验。 (1)拟合优度(Goodness-of-Fit)
Koenker 和Machado(1999)提出了分位数回归的拟合优度的概念。它与一般回归分析中的R 2
很类似。
假设分位数回归直线为
? ?(τ) =X 'β y (τ)
?=(β??将解释变量矩阵和参数向量都分为两部分,即X =(1, Z ) '和β(τ) 0(τ) , β1(τ) ) ',且有 ?(τ) =β0(τ) +Z 'β1(τ) y
定义:
?=min[- Q (τ)
t :y t
∑?∑?
T
T
??(1-τ)(y t -β0(τ) -Z 'β1(τ) ) +
T
t :y t ≥X 'β(τ)
??(y t -β0(τ) -Z 'β1(τ) ) ] (15.22) ∑τ
?
T
~
Q (τ) =min[-
?(1-τ)(y t -β0(τ) ) +
t :y t
∑?τ(y t -β?0(τ) ) ] (15.23)
式(15.22)和(15.23)分别表示无约束分位数回归目标函数(最小绝对离差和)和约束的分位数回归
目标函数(最小绝对离差和)的极小值。无约束目标函数中的减项既包含常数项也包含所有回归因子。约束目标函数中的减项仅包含常数项,其他参数都约束为零。则Koenker 和Machado 拟和优度准则表达式如下:
?Q (τ)
R (τ) =1- (15.24)
Q (τ)
*
~*?≤Q 很明显,上述统计量与传统的R 2非常相似。因为Q (τ) (τ) ,所以R (τ) 的值在0和1之间,解~**?越远远小于Q 释变量的作用越强,Q (τ) ,R (τ) 越接近1。反之越接近0。所以R (τ) 可用来(τ)
考察解释变量对被解释变量第τ分位数回归拟和的好坏。
(2)拟似然比检验(Quasi-Likelihood Ratio Tests)
Koenker 和Machado(1999)根据目标函数在施加约束条件前后得到的两个极小值构造了两个拟似然比检验统计量(QLR )。这两个拟似然比检验也称作分位数ρ检验(quantile -ρ tests) 。两统计量的表达式如下:
L T (τ) =
~?) 2(Q (τ) -Q (τ)
τ(1-τ) s (τ) ?2Q (τ)
(15.25)
~
Q (τ)
ΛT (τ) =) (15.26)
τ(1-τ) s (τ) Q (τ)
~
两个统计量都渐近服从自由度为q 的χ2分布,其中q 是原假设目标函数中约束条件的个数。Q (τ)
?分别代表约束的和无约束目标方程的极小值。 和Q (τ)
另外,两个统计量的分母都含有稀疏项s (τ) ,上面给出的稀疏项s (τ) 的3种计算方法都可在式(15.25)和(15.26)中使用。EViews 估计的是其在备择假设下的估计量。
使用上述两统计量的前提是必须满足分位数密度函数s (τ) 与解释变量X 不相关。然而,尽管
有时并不满足独立同分布的假设,EViews 在进行分位数回归的时候,不管选择何种估计参数渐近分布的方法,总会估计稀疏函数s (τ) ,从而构造拟似然比(QLR )检验统计量。因此,这种检验方法与下面的Wald 统计量相比稳健性较差。
(3)Wald 检验
给定分位数回归参数估计量的渐近方差协方差矩阵,我们就可以构造Wald 形式的统计量进行各种约束形式的参数检验。
31.3.5 系列分位数回归检验
前面的分析主要集中在单个分位数回归模型的假设检验上,而有些时候也需要对一系列分位数回归的回归系数进行联合检验。比如,需要通过检验不同分位数模型的斜率是否相等来判断一个模型是否具有位移特征。同时考虑多个分位数回归式称作系列分位数回归分析(quantile process testing )。EViews 在做单方程分位数回归的同时,有专门命令执行系列分位数回归分析。
操作路径是在一个分位数回归估计结果窗口,点击View 键,选Quantile Process/Process Coefficients 功能。
定义系列分位数回归系数列向量为, β=(β(τ1) ' , β(τ2) ' , β(τm ) ' ) ' (15.27) 则有
?-β) ~N (0, Ω) (15.28) n (β
其中Ω由如下形式的块矩阵Ωij (km×km) 组成:
Ωij =[min(τi , τj ) -τi τj ]H -1(τi ) JH -1(τj ) (15.29) i , j =1, 2, … m . k 为方程待估参数个数。其中J 的表达式见(4)式。H 的表达式见(15.19)或(15.20)式,
取决于选择的估计方法。特别的,当误差项独立同分布的假设成立时,Ω简化为: Ω=Ω0?J (15.30) 其中Ω0中的元素如下: ωij =
m in(τi , τj ) -τi τj f (F -1(τi )) f (F -1(τj ))
(15.31)
i , j =1, 2, … k . 除了以上的方法以外,Ω的估计量还可以由任何一种自举方法得到。
(1)斜率相等检验
Koenker 和Bassett(1982a)提出了一种对异方差很稳健的判断不同分位数回归方程斜率是否相等的检验。零假设如下:
H 0:β1(τ
1)
=β1(τ
2)
= =β1(τ
m )
其中β1指常数项以外的解释变量所对应的(k -1) 维参数列向量。因此,零假设共含有(k -1) (m -1) 个约束条件。接下来构造Wald 形式的统计量检验零假设是否成立,它渐近服从自由度为(k -1) (m -1)的χ2分布。
(2)对称性检验
将Newey 和Powell(1987)检验最小二乘估计量对称性的方法扩展到分位数回归中。假设我们要检验的分位数回归模型有m 个,m 是奇数,且中间值τ (m +1)/2是0.5,其他τ都关于0.5对称,即τj =1? τm-j +1, j =1,…,(m -1)/2。参数估计量按照τk 的大小排序。则对称性检验的零假设为:
H 0:
β(τj ) +β(τm -j +1)
2
=β(0.5) (15.32)
其中j =1, …, (m ?1)/2。m 是奇数,代表分位数回归个数。即关于0.5对称的分位数回归参数估计量的两两平均值等于中位数回归参数估计量。
我们可以构造Wald 形式的统计量检验上述k (m -1)/2个约束条件是否成立。该统计量服从自由度为k (m ?1)/2的χ2分布。另外,Newey 和Powell 指出,如果我们已知随机误差项服从独立同分布,但不一定对称,则我们只需检验常数项的对称性。即
β0(τ) +β0(τm -j +1)
j
=β0(0.5) (15.33) H 0:
2
这时约束条件减少为(m -1)/2个。
15.6分位数的计算与分位数回归的EViews 操作
(1)分位数的计算
对一个离散的随机变量y t ,取其容量为T 的样本序列,计算第τ分位数的方法如下:
首先将数据从小到大排序,标号为i ,i =1, 2, …, T 。然后利用下表所列的方法计算随机变量y t 的第τ分位数的排列序号的i ;如果i 为整数,则随机变量y t 的第τ分位数即为y i ,如果i 不是整数,则随机变量y t 的第τ分位数为:
y (τ) = y [i ] + (i ? [i ])( y[i ]+1 ? y[i ])
其中[i ]表示不大于i 的最大整数。给定一个具体的随机变量y t ,对于一个容量为T 的样本,则y t 的第τ分位数的序号i 的计算方法如下。在大样本情况下,各方法收敛到同一值。
Rankit
Ordinary Vander Waerden
Blom
Tukey
Gumbel
(τ?1/2)/T τ/T τ/(T +1) (τ?3/8)/ (T +1/4) (τ?1/3)/ (T +1/3) (τ?1)/ (T ?1)
计算分位数的EViews 6.0的命令为:scalar q =@quantile(y , τ, s ) ,其中y 表示求分位数的序列;τ表示要取的分位数;s 取1~6依次表示上表中6种计算方法,计算所得结果存入标量q 中。
例:打开6garch-03文件,在空白处键入命令: scalar q=@quantile(DASH , .5,1) scalar q=@quantile(DASH , .25,1)
意即对序列DASH 求中位数。得结果DASH (0.5)= -0.78,DASH 序列的中位数是-0.78。DASH (0.25)= -13.33,DASH 序列的第0.25分位数是-13.33。
用DASH 画分位数图如下。打开DASH t 序列窗口,点击View 键选Graph 功能。在打开的Graph Option 窗口,Type 选择页的Specifi 选择框选Distribution ,在Details 的Distribution 选择框中选Emprical Quantile如图。点击“确定”键,得分位数图如图。
(2)分位数回归
主要包括3部分内容。(1)介绍怎样进行分位数回归。(2)对输出结果的分析。(3)对分位数回归相关功能键的介绍。
在EViews 中进行分位数回归的路径有两个,分别是
(1)点击主选单中的Quick 键,选Equation Estimation,弹出Equation Estimation窗口。 或者
(2)点击主选单中的Object 键,选New Object,Equation ,弹出Equation Estimation窗口。 在该窗口的Method 下拉选单中,选择如图所示的选项QREG-Quantile Regression(including LAD) ,EViews 将打开如图所示的分位数回归对话框(Equation Estimation)。
图1
Equation Estimation(方程估计)窗口包括两个选项模块,一个是Specification (设定方程),一个是Options (选项)。
可以在Equation specification (方程设定)框中输入要估计的表达式。同一般线性回归模型一样,它可以是一行用空格隔开的被解释变量和解释变量(如图1所示),也可以是一个明确的参数为线性的表达式。
Equation Estimation(方程估计)窗口与OLS 估计的Equation Estimation(方程估计)窗口相比,只多了对话框quantile to estimate的选项。在该处填入要估计的分位数。系统默认为0.5,即做中位数回归(LAD )。用户可以选择任意一个0和1之间的数(当数值接近0和1时估计会变得困难)。
激活Options (选项)模块(点击对话框上的Options (选项))。得到如图2的quantile regression Options (分位数回归选择)选择框、Iteration control(迭代控制)选择框和Bootstrap settings(自举设定)选择框。
quantile regression Options对话框中的选择主要包括三部分。
图2
(1)Coefficient covariance(系数估计量方差协方差矩阵)选项框
其下拉选单中包括三个选项:Ordinary (IID),Huber-Sandwich 和Bootstrap ,代表了可选的估计回归系数估计量方差协方差矩阵的方法(具体介绍见15.4节)。EViews 默认的是Huber-Sandwich 方法。
(2)Weight (权数)选项框
可以输入作为权重的序列或者一个序列的表达式,从而对估计式加权。(用于WLS 估计) (3)Sparsity Estimation(稀疏函数估计)选项区 其中包括5种选择框。稀疏函数的介绍见15.4节。
◇ Method (方法)选项框。
当第一个选项框Coefficient covariance 中选项为Ordinary (IID)或Bootstrap 时,Method (方法)选项框中包括三个选项:Siddiqui (mean fitted), Kernel (residual)和Siddiqui (residual)。
当Coefficient covariance 选项框中选项为Huber-Sandwich 时,这里的Method 选项框中只包括Siddiqui (mean fitted)和 Kernel (residual)两种选择。
◇ Bandwidth Method(带宽)选项框。
其下拉选单中包括四个选择,即Bofinger (1975),Hall-Sheather (1988)和Chamberlain (1994)计算带宽方法,或者你自己给出一个特定的带宽。
◇ Size (置信尺度)选项框。
当选择Hall-Sheather 和Chamberlain 方法时,置信度的选择默认为0.05。 ◇ Quantile Method(分位数方法)选项框。 EViews 提供了六种求解经验分位数的方法。 ◇ Kernel (核函数)选项框。
表示核函数的选用种类。EViews 中可以选择的核密度函数有Epanechnikov 核函数、均匀核函数(Uniform)、三角核函数(Triangular)、二权核函数(Biweight)、三权核函数(Triweight)、正态核函数(Normal)和余弦核函数(Cosinus)。
注意,不管系数方差协方差矩阵(Coefficient covariance)是否会用到,每次进行分位数回归时,系统都会自动给出一个稀疏函数估计值。
Iteration control(迭代控制)选项框包括3个选项。 (1)Max (最大)。迭代的最大次数,默认为500。 (2)Starting (初始值)。表示迭代的初始值,默认为0,也可以选择其他选项,如下拉选单中的OLS ,即用OLS 估计量作为初始值进行迭代。
(3)Display settings(设定显示)。选择是否需要在输出结果中给出这些设置。 Bootstrap settings(自举设定) (1)Method (方法)。代表不同的自举方法。EViews 提供了四种方法,分别是Residual, XY-pair, MCMB, MCMB-A。默认方法为XY-pair 方法。
(2)Replications (循环次数)。EViews 默认为100次。用户可以自己设定次数。 (3)No. of obs(自举样本容量)。空白表示与原样本容量一致。Koenker(2005)的研究表明,选择自举样本容量小于数据样本容量时,能够获得更加准确的结果,特别是当数据样本容量较大时。
(4)output (输出)。在这里键入一个名称可以得到自举的样本矩阵。 (5)Random generator(生成随机数)和seed (种子)。本选项用于控制产生随机数。其中前者用于选择随机数产生方法,seed 用于选择随机数种子,Clear (清除)按钮用于清空以往选定的随机数种子。
估计结果。
按照EViews 默认设置得到的一个分位数回归估计结果如下:
输出结果上部给出的是估计设定,其中包括(按顺序) 被解释变量(DASH )、方法:分位数回归(中位数)、操作日期、样本范围、样本容量(421)、标准误差和方差协方差矩阵估计方法 (Huber-Sandwich方法) 、稀疏函数的估计方法(Kernel 方法)、带宽方法(Hall-Sheather 方法,带宽=0.12963)以及对估计结果的评价。
输出结果中部给出的是回归系数估计量、标准差、t 统计量及其相应p 值,这与OLS 估计完全一样。可以看出,上述回归系数估计量都具有显著性。在中位数回归关系条件下,B 股收益DBSH 每增加一个单位,A 股收益DASH 平均增加3.38个单位。
输出结果下部给出的是对分位数回归估计式的评价统计量。分别为 Pseudo R-squared:伪拟合优度(伪R 2) , Adjusted R-squared:调整的伪拟合优度,
S.E. of regression:分位数回归式的标准误差,
Quantile dependent var :分位数回归式中只有常数项存在的系数估计值(也即被解释变量的
分位数估计值)。
Objective :目标函数极小值,
Objective (const. only):分位数回归式中只有常数存在的目标函数极小值, Sparsity :分位数密度函数(稀疏函数)估计值(本例是用核估计法计算的)。 Quasi-LR statistic:准似然比估计量的值
Prob (Quasi-LR stat):准似然比估计量的值所对应的概率值。
此外,由于这里使用的是Huber-Sandwich 方法,因此稀疏函数值(Sparsity)并没有用来计算参数估计量标准差。
与上述结果类似,我们也可以通过改变估计设定,运用自举方法获得参数估计量的方差协方差矩阵。例如选择MCMB-A 方法进行自举,并且将自举次数增加至500。对于稀疏函数的计算方法,选择Siddiqui(mean fitted),点击OK 键,得到新设定所对应的估计结果。 分位数回归中的Views 和Procs 功能键。
分位数回归方程窗口中的大部分Views 和Procs 功能都与OLS 回归相同,下面对一些计算细节其进行必要的补充说明。
使用上述功能时需要注意以下计算细节:
(1)这里的残差是指某一特定分位数回归函数条件下的残差,计算公式为
?; ?(τ) t =y t -X 'β u (τ)
标准化残差指用自由度调整过的残差的标准误差。而在计算QLR 统计量时则使用的是Koenker ?。 ?(τ) =T -1Q 和Machado(1999)给出的目标函数极小值的平均值,即σ(τ)
(2)构造Wald 检验和置信椭圆时使用的是参数估计量方差协方差矩阵的稳健估计量。
(3)进行遗漏和多余变量检验(omitted and redundant tests)以及Ramsey RESET检验时,报告的都是特定约束下的QLR 统计量,因此它只有在满足稀疏函数的独立同分布假设时才是有效的。
(4)在默认状态下,EViews 在进行预测时会自动为样本外观测值扩展序列容量,用户可以在预测窗口中自行取消该选择框中的对勾。
Quantile Process选项
在分位数回归估计结果窗口中,点击View 键,选Quantile Process功能,如图。该选项包括3个子功能,Process Coefficients(系列分位数回归),Slope Equality Test(斜率相等检验),Symmetric Quantiles Test(分位数对称检验)。
注意:由于每个功能都包含了不同分位数估计方程的比较和检验,因此运行时间会稍长。尤其是当选择自举方法计算方差协方差阵时,更是这样。
(1)Process Coefficients功能(系列分位数回归系数)
在分位数回归估计结果窗口中,点击View 选Quantile Process,Process Coefficients(系列分
位数回归系数),可以得到Quantile Process对话框如图。
在该对话框中,可以通过Output 选择区选择是按表格(Table)还是按图形(Graph)输出估计结果;选择置信度(Confidence Interval Size),EViews 默认选择是输出表格形式(table)的估计结果,给出不同分位数下的系数估计量、标准差、t 统计量及其相应的p 值。
Quantiles 后面的数字用于决定要估计的分位数个数。若用数字n 表示,则输出的是第1/n ,2/n ,…,(n -1)/ n 分位数回归估计结果。EViews 默认数为10,即估计第0.1、0.2、…、0.9九个分位数回归式。通过选择分位数回归数(Quantiles )以及用图或列表形式(Output ),对不同分位数回归系数估计量及其标准差进行比较与检验。
还可以选择User-specified quantiles 并在下面的框中输入想计算的某个(些)分位数条件下的回归。
选取DASH 对常数项和DBSH 进行中位数回归,并在Quantile Process选项中选择输出95%置信水平的10个不同的分位数回归中回归系数估计值的折线图,如下所示:
点击Quantile Process对话框中的Output 选项页,得到如下的对话框:
用Output 选项页保存分位数向量,分位数回归系数估计量矩阵,分位数回归系数估计量的方差协方差矩阵。
在Quantiles Vector、Coefficient Matrix和Covariance Matrix选择框后面分别填入向量和矩阵名,EViews 将把分位数序列、分位数回归系数估计量矩阵以及分位数回归系数估计量的方差协方差矩阵保存到相应向量和矩阵名中。
(2)Slope Equality Test(斜率相等检验)
在分位数回归估计结果窗口中点击View 选Quantile Process ,Slope equality test ,可以进行Koenker 和Bassett(1982a)提出的斜率参数相等检验,见如下对话框:
在该窗口的Specification (设定)选项页中,Quantiles 后面的框中用户可以填入要检验的分位数个数,从而对上述分位数回归式的斜率估计量进行比较和检验。以填入数字n 为例,表示输出的是第1/n ,2/n ,…,(n -1)/n 分位数。
EViews 默认的是4,即检验第1/4分位数回归式与第2/4分位数回归式之间、第2/4分位数回归式与第3/4分位数回归式之间相应回归系数是否相等。用户还可以选择User-specified quantiles ,并在下面的框中输入想要检验的分位数条件下的回归系数是否相等。
该窗口中的Output 选项页的功能与Quantile Process对话框中的Output 选项页的功能相同。 仍然使用上面的案例,按照系统默认的检验,得到如下的结果,即对中位数和第0.25分位数、以及中位数与第0.75分位数回归系数是否相等这2个约束条件进行检验。
输出结果中第一部分是对Wald 检验的总结,这里斜率相等性检验的Wald 统计量为0.53,自由度2,概率为0.7670。这表明0.25、0.5、0.75分位数回归式的斜率相等。这意味着,相同条件的不同分位数回归的被解释变量拟合值的分布是相同的。
(3)对称性检验
在分位数回归估计结果窗口中点击View 选Quantile Process,Symmetric Quantiles Test,可以做以中位数回归式为中心的分位数回归式对称性检验。
依据Newey 和Powell(1987)方法,
β(τ) +β(1-τ)
=β(0.5)
2
关于对称性检验的设置可以在Specification 选项页中完成,其默认选项如图所示。如果Quantiles (分位数回归式个数)选择框中选4,则检验的约束条件只有一个,即
β(0.25) +β(0.75)
=β(0.5)
2
而如果分位数回归式拟合的是0.6,则除了以上约束条件外,还会增加如下一个约束条件的检验。
β(0.4) +β(0.6)
=β(0.5)
2
与上述两个检验类似,用户还可以选择User-specified quantiles,并在下面的框中输入要检验的分位数τ,EViews 将估计τ和1-τ两个分位数回归式,并与中位数回归式做比较。按照(β(τ) +β(1-τ) ) /2=β(0.5) 进行检验。
此外,还可以通过Intercept only或者All coefficients单选框选择是否仅对常数项还是对所有参数估计量进行对称性检验。同样,该窗口中的Output 选项页的功能与Quantile Process对话框中的Output 选项页的功能相同。
本例按照默认设置检验结果如下:
从输出结果中可以看出,这里我们检验的是第0.25和0.75分位数回归的系数估计量是否关于中位数对称。检验表明他们是对称的,概率p 值是0.80。
预测功能与OLS 单方程操作相同。
15.7 分位数回归的案例分析
案例分析:上证A 股和B 股收益关系研究(6garch-03)
上证A 股、B 股序列(ASH t 和BSH t )和上证A 股、B 股收益序列(DASH t 和DBSH t )如图。
图1 ASH t 和BSH t 序列 图2 DASH t 和DBSH t 序列
上证A 股和B 股收益序列(DASH t 和DBSH t )的单位根检验结果如下:
D 2ASH t = -0.9786 DASH t -1
(-20.0)*
DW =1.99
D 2BSH t = -0.9599 DBSH t -1
(-19.8)* DW =2.00
图3 上证A 股收益序列(DASH t )的单位根检验结果
图4 上证B 股收益序列(DBSH t )的单位根检验结果
上证A 股和B 股的差分(收益)序列(DASH t 和DBSH t )都是平稳序列。
下面分析DASH t 和DBSH t 的分布特征。打开DASH t 序列窗口,点击View 键选Graph 功能。在图形类型选择框中选Distribution (分布)。同时在选择窗右侧点击Options 键,在弹出的对话窗中点击Add 键,在弹出的Add 对话窗中选Kernel Density。点击OK 键。得直方图与核密度曲线合并图如图。分布都是单峰的,呈高峰厚尾特征。
分别求DASH t 和DBSH t 的第0.25、0.5、0.75分位数值(scalar q =@quantile(DASH , τ, 1)),τ=0.25、0.5、0.75,得下表。
DASH t DBSH t
均值 -1.40 0.04
偏度 0.77 -0.05
标准差 25.7 4.1
第0.25分位数
-13.33 -1.510
第0.5分位数 第0.75分位数
-0.78 -0.06
9.74 1.73
上证A 股和B 股收益序列散点图如下:
OLS 估计结果如下:B 股收益每增加一个单位,A 股收益增加3.93个单位。
DASH t 对DBSH t 的中位数回归结果如下。
OLS 回归、LAD (中位数)回归直线对比如图。OLS 回归直线的斜率大于LAD (中位数)回归直线的斜率,说明OLS 估计方法对离群值敏感,LAD 估计方法比OLS 估计方法稳健。
在分位数回归输出结果窗口点击View 键,选Quantile Process,Process Coefficients,在弹出的对话框中点击确定键,得9个分位数回归式如下,DBSH 的回归系数全部具有显著性。
9个分位数回归式中的回归系数值折线图如下,回归系数值在3.38和3.97之间变化。除了第0.1分位数回归式之外,其他8个分位数回归式的回归系数基本相同,说明不同分位数回归式之间,DASH 和DBSH 的回归关系基本相同。
OLS 回归直线和第0.25、0.5、0.75分位数回归直线如下图。回归直线的斜率大于中位数回归线的斜率。
斜率相等检验:
在分位数回归输出结果窗口点击View 键,选Quantile Process,Slope equality test功能,在弹出的对话框中点击确定键,得输出结果如下。经检验,第0.25和0.5分位数回归直线的斜率,第0.5和0.75分位数回归直线的斜率之间无显著性差异。 原假设有两个H 0:β1(0.25)= β1(0.5),β1(0.5) = β1(0.75)
对称性检验
原假设有两个H 0:β0(0.25)+ β0(0.75) = 2β0(0.5) ,β1(0.25)+ β1(0.75) = 2β1(0.5)
经检验,第0.25和0.75分位数回归直线的斜率以中位数斜率存在对称性。
案例分析:恩格尔定律研究。(分位数回归,) (file:6food,Murray-book 第1章)
求food 的分布直方图,第0.25、0.5、0.75分位数值(scalar q =@quantile(food, τ, 1)),τ=0.25、0.5、0.75,得下表。
0.75分位数 696.5
0.9分位数 856.2
均值 580.0
标准差 238.8
偏度 1.57
0.1分位数 324.4
0.25分位数 410.0
0.5分位数 558.5
food
F O O D
D e n s i t y
TOTAL
点击主选单中的Quick 键,选Equation Estimation ,弹出Equation Estimation 窗口。在该窗口的Method 下拉选单中,选择QREG-Quantile Regression(including LAD) ,EViews 将打开分位数回归对话框(Equation Estimation)。点击确定键,中位数回归估计结果: ?? food t = β0(0.5) +β1(0.5) total = 75.6219 + 0.5419 total
∧
(3.3) (18.9)
在分位数回归输出结果窗口点击View 键,选Quantile Process,Process Coefficients,在弹出的对话框中点击确定键,得9个分位数回归式如下,DBSH 的回归系数全部具有显著性。
第10、30、50、70、90分位数回归与均值回归拟合直线的比较
系列分位数回归系数(Quantile Process Estimates)估计值及其置信区间(95%)。
用第10、20、30、40、50、60、70、80、90分位数回归得到的回归系数估计值各自连成的线(蓝线)以及相应的95%置信度的置信区间(红线)。估计结果显示,分位数回归式随着分位数的增加,食品支出与总支出的回归系数也增加,数据存在异方差。 分位数系数相等性检验结果如下。第0.25与0.5分位数回归式之间,第0.5与0.75分位数回归式之间的斜率不相等。这也从一个角度显示数据存在递增型异方差。 原假设有两个H 0:β1(0.25)= β1(0.5),β1(0.5) = β1(0.75)
分位数系数对称性检验如下。第0.25与0.75分位数回归式的斜率对称。第0.25与0.75分位数回归式对于中位数回归式来说,具有对称性。
原假设有两个H 0:β0(0.25)+ β0(0.75) = 2β0(0.5) ,β1(0.25)+ β1(0.75) = 2β1(0.5)
下面研究对数数据情形。 中位数回归估计结果:
??Lnfood t = β0(0.5) +β1(0.5) Lntotal = 0.4991 + 0.8595 Lntotal
∧
(1.9) (21.7)
系列分位数(第20、40、60、80、90分位数)回归结果
第20、40、60、80、90分位数回归与均值回归拟合直线的比较
系列分位数回归系数(Quantile Process Estimates)估计值及其置信区间(95%)。
用第10、20、30、40、50、60、70、80、90分位数回归得到的回归系数估计值各自连成的线(蓝线)以及相应的95%置信度的置信区间(红线)。
分位数系数相等性检验。经检验,4个分位数回归式的斜率相等。
分位数系数对称性检验如下。经检验,两对分位数回归式的斜率对称。
范文五:分位数回归
2、不同分位点拟合曲线的比较
# 散点图
attach(engel) # 打开engel数据集,直接运行其中的列名,就可以调用相应列 plot(income,foodexp,cex=0.25,type="n", # 画图,说明?
xlab="Household Income", ylab="Food Expenditure")
points(income,foodexp,cex=0.5,col="blue") # 添加点,点的大小为0.5 abline( rq(foodexp ~ income, tau=0.5), col="blue" ) # 画中位数回归的拟合直线,颜色蓝 abline( lm(foodexp ~ income), lty = 2, col="red" ) # 画普通最小二乘法拟合直线,颜色红 taus = c(0.05, 0.1, 0.25, 0.75, 0.9, 0.95)
for(i in 1:length(taus)){ # 绘制不同分位点下的拟合直线,颜色为灰色
abline( rq(foodexp ~ income, tau=taus[i]), col="gray" )
}
detach(engel)
3、穷人和富人的消费分布比较
# 比较穷人(收入在10%分位点的那个人)和富人(收入在90%分位点的那个人)的估计结果
# rq函数中,tau不在[0,1]时,表示按最细的分位点划分方式得到分位点序列 z = rq(foodexp ~ income, tau=-1)
z$sol # 这里包含了每个分位点下的系数估计结果
x.poor = quantile(income, 0.1) # 10%分位点的收入
x.rich = quantile(income, 0.9) # 90%分位点的收入
ps = z$sol[1,] # 每个分位点的tau值
qs.poor = c( c(1,x.poor) %*% z$sol[4:5,] ) # 10%分位点的收入的消费估计值 qs.rich = c( c(1,x.rich) %*% z$sol[4:5,] ) # 90%分位点的收入的消费估计值 windows(, 10,5)
par(mfrow=c(1,2)) # 把绘图区域划分为一行两列
plot(c(ps,ps),c(qs.poor,qs.rich),type="n", # type=”n”表示初始化图形区域,但不画图
xlab=expression(tau), ylab="quantile")
plot(stepfun(ps,c(qs.poor[1],qs.poor)), do.points=F,
add=T)
plot(stepfun(ps,c(qs.poor[1],qs.rich)), do.points=F,
add=T, col.hor="gray", col.vert="gray")
ps.wts = ( c(0,diff(ps)) + c(diff(ps),0) )/2
ap = akj(qs.poor, z=qs.poor, p=ps.wts)
ar = akj(qs.rich, z=qs.rich, p=ps.wts)
plot(c(qs.poor,qs.rich), c(ap$dens, ar$dens),
type="n", xlab="Food Expenditure", ylab="Density")
lines(qs.rich,ar$dens,col="gray")
lines(qs.poor,ap$dens,col="black")
legend("topright", c("poor","rich"), lty=c(1,1),
col=c("black","gray"))
上图表示收入(income)为10%分位点处(poor,穷人)和90%分位点处(rich,富人)的食品支出的比较。从左图可以发现,对于穷人而言,在不同分位点估计的食品消费差别不大。而对于富人而言,在不同分位点对食品消费的差别比较大。右图反应了穷人和富人的食品消费分布曲线。穷人的食品消费集中于400左右,比较陡峭;而富人的消费支出集中于800到1200之间,比较分散。
(四)模型比较
# 比较不同分位点下,收入对食品支出的影响机制是否相同
fit1 = rq(foodexp ~ income, tau = 0.25)
fit2 = rq(foodexp ~ income, tau = 0.5)
fit3 = rq(foodexp ~ income, tau = 0.75)
anova(fit1,fit2,fit3)
结果:
Quantile Regression Analysis of Deviance Table
Model: foodexp ~ income
Joint Test of Equality of Slopes: tau in { 0.25 0.5 0.75 }
Df Resid Df F value Pr(>F)
1 2 703 15.557 2.449e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
其中P值远小于0.05,故不同分位点下收入对食品支出的影响机制不同。 (五)残差形态的检验
也可以理解为是比较不同分位点的模型之间的关系。主要有两种模型形式:
(1)位置漂移模型:不同分位点的估计结果之间的斜率相同或近似,只是截距不同;表现为不同分位点下的拟合曲线是平行的。
(2)位置-尺度漂移模型:不同分位点的估计结果之间的斜率和截距都不同;表现为不同分位点下的拟合曲线不是平行的。
# 残差形态的检验
source("C:/Program Files/R/R-2.15.0/library/quantreg/doc/gasprice.R") x = gasprice
n = length(x)
p = 5
X = cbind(x[(p-1):(n-1)],x[(p-2):(n-2)],x[(p-3):(n-3)],x[(p-4):(n-4)])
y = x[p:n]
# 位置漂移模型的检验
T1 = KhmaladzeTest(y~X, taus = -1, nullH="location")
T2 = KhmaladzeTest(y~X, taus = 10:290/300,
nullH="location", se="ker")
# 位置尺度漂移模型的检验
T3 = KhmaladzeTest(y~X, taus = -1, nullH="location-scale") T4 = KhmaladzeTest(y~X, taus = 10:290/300,
nullH="location-scale", se="ker")
结果:运行T1,可以查看其检验结果。其中nullH表示原假设为“location”,即原假设为
位置漂移模型。Tn表示模型整体的检验,统计量为4.8。THn是对每个自变量的检验。比较
T1和T3的结果(T3的原假设为“位置尺度漂移模型”),T1的统计量大于T3的统计量,
可见相对而言,拒绝“位置漂移模型”的概率更大,故相对而言“位置尺度漂移模型”更加
合适一些。
> T1
$nullH
[1] "location"
$Tn
[1] 4.803762
$THn
X1 X2 X3 X4
1.0003199 0.5321693 0.5020834 0.8926828
attr(,"class")
[1] "KhmaladzeTest"
> T3
$nullH
[1] "location-scale"
$Tn
[1] 2.705583
$THn
X1 X2 X3 X4
1.2102899 0.6931785 0.5045163 0.8957127
attr(,"class")
[1] "KhmaladzeTest"
(六)非线性分位数回归
这里的非线性函数为Frank copula函数。
## Demo of nonlinear quantile regression model based on Frank copula vFrank <- function(x,="" df,="" delta,="" u)="" #="" 某个非线性过程,得到的是[0,1]的值="" -log(1-(1-exp(-delta))/(1+exp(-delta*pt(x,df))*((1/u)-1)))/delta="" #="">->
FrankModel <- function(x,="" delta,="" mu,sigma,="" df,="" tau)="" {="">->
z <- qt(vfrank(x,="" df,="" delta,="" u="tau)," df)="">->
mu + sigma*z
}
n <- 200="" #="" 样本量="">->
df <- 8="" #="" 自由度="">->
delta <- 8="" #="" 初始参数="">->
set.seed(1989)
x <- sort(rt(n,df))="" #="" 生成基于t分布的随机数="">->
v <- vfrank(x,="" df,="" delta,="" u="runif(n))" #="" 基于x生成理论上的非参数对应值="" y="">-><- qt(v,="" df)="" #="" v="" 对应的t分布统计量="" windows(5,5)="">->
plot(x, y, pch="o", col="blue", cex = .25) # 散点图
Dat <- data.frame(x="x," y="y)" #="" 基本数据集="">->
us <- c(.25,.5,.75)="">->
for(i in 1:length(us)){
v <- vfrank(x,="" df,="" delta,="" u="us[i])">->
lines(x, qt(v,df)) # v为概率,计算每个概率对应的T分布统计量 }
cfMat <- matrix(0,="" 3,="" length(us)+1)="" #="" 初始矩阵,用于保存结果的系数="" for(i="" in="" 1:length(us))="" {="">->
tau <- us[i]="">->
cat("tau = ", format(tau), ".. ")
fit <- nlrq(y="" ~="" frankmodel(x,="" delta,mu,sigma,="" df="8," tau="tau)," #="" 非参数模型="">->
data = Dat, tau = tau, # data表明数据集,tau分位数回归的分位点
start= list(delta=5, mu = 0, sigma = 1), # 初始值
trace = T) # 每次运行后是否把结果显示出来
lines(x, predict(fit, newdata=x), lty=2, col="red") # 绘制预测曲线
cfMat[i,1] <- tau="" #="" 保存分位点的值="">->
cfMat[i,2:4] <- coef(fit)="" #="" 保存系数到cfmat矩阵的第i行="">->
cat("\n") # 如果前面把每步的结果显示出来,则每次的结果之间添加换行符
}
colnames(cfMat) <- c("分位点",names(coef(fit)))="" #="" 给保存系数的矩阵添加列名="" cfmat="">->
结果:
拟合结果:(过程略)
> cfMat
分位点 delta mu sigma
[1,] 0.25 14.87165 -0.20530041 0.9134657
[2,] 0.50 16.25362 0.03232525 0.9638209
[3,] 0.75 12.09836 0.11998614 0.9423476 (七)半参数和非参数分位数回归
非参数分位数回归在局部多项式的框架下操作起来更加方便。可以基于以下函数。
# 2-7-1 半参数模型----
fit<-rq(y~bs(x,df=5)+z,tau=.33)>-rq(y~bs(x,df=5)+z,tau=.33)>
# 其中bs()表示按b-spline的非参数拟合
# 2-7-2 非参数方法
lprq <-function(x,y,h,m=50,tau=0.5){>-function(x,y,h,m=50,tau=0.5){>
# 这是自定义的一个非参数计算函数,在其他数据下同样可以使用
xx<-seq(min(x),max(x),length=m) #="" m个监测点="">-seq(min(x),max(x),length=m)>
fv<-xx>-xx>
dv<-xx>-xx>
for(i in 1:length(xx)){
z<-x-xx[i]>-x-xx[i]>
wx<-dnorm(z )="" #="" 核函数为正态分布,dnorm计算标准正态分布的密度值="">-dnorm(z>
r<-rq(y~z,weights=wx,tau=tau, #="" 上面计算得到的密度值为权重="">-rq(y~z,weights=wx,tau=tau,>
ci=FALSE)
fv[i]<-r$coef[1]>-r$coef[1]>
dv[i]<-r$coef[2]>-r$coef[2]>
}
list(xx=xx,fv=fv,dv=dv) # 输出结果
}
library(MASS)
data(mcycle)
attach(mcycle)
windows(5,5) # 非参数的结果一般是通过画图查看的 plot(times,accel,xlab="milliseconds",ylab="acceleration") hs<-c(1,2,3,4) #="" 选择不同窗宽进行估计="">-c(1,2,3,4)>
for(i in hs){
h=hs[i]
fit<-lprq(times,accel,h=h,tau=0.5) #="" 关键拟合函数="">-lprq(times,accel,h=h,tau=0.5)>
lines(fit$xx,fit$fv,lty=i)
}
legend(45,-70,c("h=1","h=2","h=3","h=4"),
lty=1:length(hs))
结果:
# 2-7-3 非参数回归的另一个方法----
# 考察最大的跑步速度与体重的关系
data(Mammals)
attach(Mammals)
x<-log(weight) #="" 取得自变量的值="">-log(weight)>
y<-log(speed) #="" 取得因变量的值="">-log(speed)>
plot(x,y,xlab="Weightinlog(Kg)",ylab="Speedinlog(Km/hour)",
type="n")
points(x[hoppers],y[hoppers],pch="h",col="red") points(x[specials],y[specials],pch="s",col="blue") others<-(!hoppers&!specials)>-(!hoppers&!specials)>
points(x[others],y[others],col="black",cex=0.75) fit<-rqss(y~qss(x,lambda=1),tau=0.9) #="" 关键的拟合步骤="" plot(fit,add="T)" #="" add="T表示在上图中添加这里的曲线">-rqss(y~qss(x,lambda=1),tau=0.9)>
(八)分位数回归的分解
# MM2005分位数分解的函数,改变参数可直接使用 MM2005 = function(formu,taus, data, group, pic=F){
# furmu 为方程,如foodexp~income
# taus 为不同的分位数
# data 总的数据集
# group 分组指标,是一个向量,用于按行区分data
# pic 是否画图,如果分位数比较多,建议不画图
engel1 = data[group==1,]
engel2 = data[group==2,]
# 开始进行分解
fita = summary( rq(formu, tau = taus, data = engel1 ) )
fitb = summary( rq(formu, tau = taus, data = engel2 ) )
tab = matrix(0,length(taus),4)
colnames(tab) = c("分位数","总差异","回报影响","变量影响")
rownames(tab) = rep("",dim(tab)[1])
for( i in 1:length(taus) ){
ya = cbind(1,engel1[,names(engel1)!=formu[[2]]] ) %*% fita[[i]]$coef[,1]
yb = cbind(1,engel2[,names(engel2)!=formu[[2]]] ) %*% fitb[[i]]$coef[,1]
# 这里以group==1为基准模型,用group==2的数据计算反常规模型拟合值
ystar = cbind(1,engel2[,names(engel2)!=formu[[2]]] ) %*% fita[[i]]$coef[,1]
ya = mean(ya)
yb = mean(yb)
ystar = mean(ystar)
tab[i,1] = fita[[i]]$tau
tab[i,2] = yb - ya
tab[i,3] = yb - ystar # 回报影响,数据相同,模型不同:模型机制的不同所产生的差异
tab[i,4] = ystar - ya # 变量影响,数据不同,模型相同:样本点不同产生的差异
}
# 画图
if( pic ){
attach(engel)
windows(5,5)
plot(income, foodexp, cex=0.5, type="n",main="两组分位数回归结果比较")
points(engel1, cex=0.5, col=2)
points(engel2, cex=0.5, col=3)
for( i in 1:length(taus) ){
abline( fita[[i]], col=2 )
abline( fitb[[i]], col=3 )
}
detach(engel)
}
# 输出结果
tab
}
# 下面是用一个样本数据做的测试
data(engel)
group = c(rep(1,100),rep(2,135)) # 取前100个为第一组,后135个第二组 taus = c(0.05,0.25,0.5,0.75,0.95) # 需要考察的不同分位点
MM2005(foodexp~income, taus, data = engel, group=group, pic=F) # 参数说明,见? 说明:? 自编函数MM2005的参数说明:
函数调用形式MM2005 (formu, taus, data, group, pic=F),其中 # furmu 为回归方程,如foodexp~income;
# taus 为不同的分位数,如taus=c(0.05,0.5,0.95);
# data 总的数据集,如上例中的engel数据集;
# group 分组指标,是一个向量,用于按行区分data,第一组为1,第二组为2;目前仅能分两组;
# pic 逻辑参数:是否画图。如果分位数比较多,建议不画图;默认不画图,pic=F;如果想画图,则添加参数pic=T。
? 最终结果:
> MM2005(foodexp~income, taus, data = engel, group=group, pic=F)
分位数 总差异 回报影响 变量影响
0.05 -30.452061 -72.35939 41.90733
0.25 -2.017317 -46.20125 44.18394
0.50 30.941212 -23.24042 54.18163
0.75 43.729025 -15.76283 59.49185
0.95 52.778985 -11.29932 64.07830
第三节 一个案例
这里使用某年份流动人口中个人消费和个人收入的数据作为案例,进行整体性分析。
(一)初始化程序
基础数据是dta格式(stata软件的数据格式),所以要用foreign包里的read.dta函数读
取。这部分的主要内容包括:
1、修改工作目录至程序和数据所在文件夹,注意不能使用“\”,必须使用“/”或“\\”;
2、清除其他变量、清理内存;
3、加载分位数回归所需的包quantreg和导入数据所需的foreign包;
4、加载4-functions.R文件中的R程序,这里有三个有用的函数;
5、读取案例数据20121218.dta;
6、初始化处理:去掉含有缺失值的行,检查样本个数,检查列名,整理列名,进行对
数变换并生成新的数据集dat2。
# 一个案例
setwd("F:/实践2/2012_分位数回归等两个任务/分位数回归方法夹") # 设置工作目录 setwd(””)
rm(list=ls())
gc()
library(quantreg)
library(foreign)
source("4-functions.R") # 将4-functions.R中的函数放入内存中 dat = read.dta("20121218.dta")
dat = na.omit(dat) # 去掉有缺失值的行
dim(dat)
# 这个数据包含10000个样本,3个指标
names(dat)
# [1] "urtype" "q413" "q410a"
# q413 每月家庭总收入
# q410a 家庭食品支出
# urtype 样本类型,1农村或2城镇
# 这里考察上个月的收入与食品支出的关系,并分性别进行讨论 # 把变量的名称改一下
names(dat) = c("urtype","income","foodexp")
# 进行log变换
dat2 = dat[dat$income>0 & dat$foodexp>0,]
dat2[,"lnincome"] = log( dat2[,"income"] )
dat2[,"lnfoodexp"] = log( dat2[,"foodexp"] )
dat2 = na.omit(dat2) # 去掉log变换以后的缺失值行
dim(dat2)
(二)画散点图
画散点图的主要目的是检查数据的分布是否需要进行对数变换,另外查看一下使用什么样的模型更合适。本案例中进行对数变换是要合适一些。
# 散点图,观察是否需要log变换
attach(dat)
windows(5,5)
# 不变换
plot(income,foodexp,cex=0.5,main="散点图")
detach(dat)
# 变换
attach(dat2)
windows(5,5)
plot(lnincome,lnfoodexp,cexp=0.5,main="散点图(对数变换以后)")
detach(dat2)
结果:
(三)分位数回归的参数估计结果
调用rq函数进行分位数回归估计。
1、这里选取的分位点包括0.05、0.25、0.5、0.75和0.9。分别对整体数据、城镇和农村分别建模。如果数据量较大,建议method参数为“fn”或“pfn”,运算效率会高些。
2、用summary()函数可以得到更加详细的统计结果。
3、用笔者编写的tabs()函数和tabcoef()函数可以把summary()处理过的结果进行整合。注意,笔者的这两个函数在4-funcitons.R中,只有加载了这个文件(source(“4-functions.R”))才能使用这两个函数,
# 建立基础的分位数回归模型
taus = c(0.05,0.25,0.5,0.75,0.95)
fit1 = rq( lnfoodexp ~ lnincome, data=dat2, tau=taus, method="fn" ) fit2 = rq( lnfoodexp ~ lnincome, data=dat2[dat2$urtype==1,],
tau=taus, method="fn" )
fit3 = rq( lnfoodexp ~ lnincome, data=dat2[dat2$urtype==2,],
tau=taus, method="fn" )
# 所有数据放在一起是参数结果
s1 = summary(fit1)
# 男性
s2 = summary(fit2)
# 女性
s3 = summary(fit3)
# 参数结果的直接比较
tabs(s1) # 这里的tabs函数是笔者自己写的,方便显示结果 tabs(s2)
tabs(s3)
# 具体的系数估计值及假设检验
tabcoef(s1) # 这里的tabcoef函数是笔者自己写的,方便显示结果 tabcoef(s2)
tabcoef(s3)
结果:
表3-1-a 整体数据在不同分位点的系数(tabs(s1))
Quantiles (Intercept) lnincome
0.05 4.26 0.50
0.25 4.55 0.52
0.50 5.01 0.50
0.75 5.43 0.49
0.95 5.93 0.52
表3-1-b 农村在不同分位点的系数(tabs(s2))
Quantiles (Intercept) lnincome
0.05 3.69 0.58
0.25 4.31 0.55
0.50 4.69 0.55
0.75 5.01 0.55
0.95 5.59 0.58
表3-1-c 城镇在不同分位点的系数(tabs(s3))
Quantiles (Intercept) lnincome
0.05 4.72 0.42
0.25 5.33 0.40
0.50 5.92 0.35
0.75 6.70 0.28
0.95 7.08 0.32
表3-1-d 整体在不同分位点的系数及其检验结果(tabcoef(s1))
Quantiles Value Std. Error t value Pr(>|t|)
0.05 4.26 0.093 45.869 0.000
0.05 0.50 0.014 36.819 0.000
0.25 4.55 0.057 79.293 0.000
0.25 0.52 0.008 63.960 0.000
0.50 5.01 0.060 82.969 0.000
0.50 0.50 0.009 58.144 0.000
0.75 5.43 0.076 71.414 0.000
0.75 0.49 0.011 42.916 0.000
0.95 5.93 0.185 32.005 0.000
0.95 0.52 0.026 19.734 0.000
表3-1-e 农村在不同分位点的系数及其检验结果(tabcoef(s2))
Quantiles Value Std. Error t value Pr(>|t|)
0.05 3.69 0.082 44.741 0.000
0.05 0.58 0.012 47.113 0.000
0.25 4.31 0.080 54.037 0.000
0.25 0.55 0.011 50.725 0.000
0.50 4.69 0.076 61.495 0.000
0.50 0.55 0.011 49.483 0.000
0.75 5.01 0.096 52.016 0.000
0.75 0.55 0.014 40.039 0.000
0.95 5.59 0.198 28.181 0.000
0.95 0.58 0.028 20.534 0.000
表3-1-f 城镇在不同分位点的系数及其检验结果(tabcoef(s3))
Quantiles Value Std. Error t value Pr(>|t|)
0.05 4.72 0.147 32.141 0.000
0.05 0.42 0.022 19.667 0.000
0.25 5.33 0.098 54.526 0.000
0.25 0.40 0.015 27.160 0.000
0.50 5.92 0.097 61.130 0.000
0.50 0.35 0.014 24.700 0.000
0.75 6.70 0.103 65.335 0.000
0.75 0.28 0.015 19.047 0.000
0.95 7.08 0.290 24.422 0.000
0.95 0.32 0.044 7.212 0.000
(四)分位数回归模型形式的检验
这里主要是检验是位置漂移模型还是位置尺度漂移模型。
# 参数结果的检验
attach(dat2)
T1 = KhmaladzeTest(lnfoodexp ~ lnincome, taus=seq(.1,.9,by = .1),
nullH="location", se="ker")
T2 = KhmaladzeTest(lnfoodexp ~ lnincome, taus=seq(.1,.9,by = .1),
nullH="location-scale", se="ker")
T1
T2
detach(dat2)
结果:根据以下结果,以位置漂移模型为原假设的统计量为2.84,以位置漂移模型为原假设
的统计量为3.29,。相对而言,更容易拒绝位置漂移模型,故按照位置漂移模型相对合适一
些,即不同分位点的斜率基本上可以认为是相同的,而截距不同。
> T1
$nullH
[1] "location"
$Tn
[1] 2.840216
$THn
[1] 2.840216
attr(,"class")
[1] "KhmaladzeTest"
> T2
$nullH
[1] "location-scale"
$Tn
[1] 3.2899
$THn
[1] 3.2899
attr(,"class")
[1] "KhmaladzeTest" (五)分位数回归的图形结果 # 画图
# 参数比较
windows(5,5)
plot(s1)
windows(5,5)
plot(s2)
windows(5,5)
plot(s3)
# 散点图
# - 总体
windows(5,5)
attach(dat2)
plot(lnincome,lnfoodexp,cex=0.25,type="n",
xlab="Personal Income", ylab="Food Expenditure",
main="Income and Food Expenditure")
points(lnincome,lnfoodexp,pch=".",col="blue")
abline( rq(lnfoodexp ~ lnincome, tau=0.5, method="fn"), col="blue" ) abline( lm(lnfoodexp ~ lnincome), lty = 2, col="red" )
taus = c(0.05, 0.1, 0.25, 0.75, 0.9, 0.95)
for(i in 1:length(taus)){
abline( rq(lnfoodexp ~ lnincome, tau=taus[i], method="fn"), col="gray" ) }
legend("bottomright",c("QR","OLS"),col=c("blue","red"),lty=1:2,bty="n") detach(dat2)
# - 城镇
windows(5,5)
dat2a = dat2[dat2$urtype==1,]
attach(dat2a)
plot(lnincome,lnfoodexp,cex=0.25,type="n",
xlab="Personal Income", ylab="Food Expenditure",
main="Income and Food Expenditure")
points(lnincome,lnfoodexp,pch=".",col="blue")
abline( rq(lnfoodexp ~ lnincome, tau=0.5, method="fn"), col="blue" ) abline( lm(lnfoodexp ~ lnincome), lty = 2, col="red" )
taus = c(0.05, 0.1, 0.25, 0.75, 0.9, 0.95)
for(i in 1:length(taus)){
abline( rq(lnfoodexp ~ lnincome, tau=taus[i], method="fn"), col="gray" ) }
legend("bottomright",c("QR","OLS"),col=c("blue","red"),lty=1:2,bty="n") detach(dat2a)
# - 农村
windows(5,5)
dat2b = dat2[dat2$urtype==2,]
attach(dat2b)
plot(lnincome,lnfoodexp,cex=0.25,type="n",
xlab="Personal Income", ylab="Food Expenditure",
main="Income and Food Expenditure")
points(lnincome,lnfoodexp,pch=".",col="blue")
abline( rq(lnfoodexp ~ lnincome, tau=0.5, method="fn"), col="blue" ) abline( lm(lnfoodexp ~ lnincome), lty = 2, col="red" )
taus = c(0.05, 0.1, 0.25, 0.75, 0.9, 0.95)
for(i in 1:length(taus)){
abline( rq(lnfoodexp ~ lnincome, tau=taus[i], method="fn"), col="gray" ) }
legend("bottomright",c("QR","OLS"),col=c("blue","red"),lty=1:2,bty="n")
detach(dat2b)
结果:
六)分位数分解
按照MM2005方法进行分解。注意的是:
1、group参数必须是包含1和2的向量,其中反事实模型的拟合值是采用1类的模型参数和2类的数据计算得到。相应地,在本案例中作为参照的拟合结果是假定农村流动人口按照城镇的决定模型,会得到什么样的食品支出。
2、data参数中只能包含模型中出现的变量。比如这里使用的变量只包括lnfoodexp和lnincome,所以只能保留这两个变量在data中。
# 按照城乡进行分解
group = as.character( dat2$urtype ) # 必须保证group是包含1和2的向量 group[group=="1"] = 1
group[group=="2"] = 2
taus = c(0.05,0.25,0.5,0.75,0.95)
data = dat2[,c("lnfoodexp","lnincome")]
res = MM2005(lnfoodexp ~ lnincome,taus, data, group, pic=F, m.rq="fn") # 分位数 总差异 回报影响 变量影响
# 0.05 -0.051089859 0.06520406 -0.11629392
# 0.25 -0.056506090 0.04647729 -0.10298338
# 0.50 -0.014632774 0.08108505 -0.09571782
# 0.75 0.009629691 0.10826533 -0.09863564
# 0.95 0.028791246 0.10291667 -0.07412543
write.csv(res,"tab-MM2005.csv",row.names=F) # 将结果保存为csv表格,方便整理 结果:
表3-2 分位数分解:按城乡的分解
分位数 总差异 回报影响 变量影响
0.05 -0.22 -0.04 -0.18
0.25 -0.18 -0.01 -0.17
0.50 -0.23 -0.06 -0.17
0.75 -0.27 -0.10 -0.17
0.95 -0.41 -0.23 -0.17
以上结果显示,无论对于高收入家庭还是低收入家庭,个体对食品支出的影响都差不多。但对于高收入家庭(分位数较高),总差异和回报影响更大。说明城乡高收入人群的食品消费差异较大,其中主要的差异来自于斜率的不同,即食品消费关于收入的敏感性不同带来的差异。