金融数据分析 课件 第7、8章 COPULA及其应用、面板计量模型与检验

第七章copula函数及其应用

学习目标

掌握Copula函数的定义、基本性质和相关性测度；熟悉常用的Copula函数的性质以及相关性分析的特点；了解如何对相依数据分析并进行相依风险度量。

本章导读

相依性建模是现代金融研究的重要领域之一，无论是相关性分析、风险关联分析、资产定价分析及信用风险分析，相依性都是一个非常重要的工具。Copula函数能够很好好地刻画金融序列变量间的非线性和非对称相依结构，并且能独立于边际分布而对变量间的相依结构进行建模。本章将详细介绍Copula函数的定义和基本性质，并介绍基于Copula函数的相依性测度，重点研究几类常用Copula函数的相依性特征，最后介绍Copula函数的估计方法及其在金融风险计量中的应用。通过本章内容学习，学生应具备整体思维，树立全局意识，强化基础理论以及科研素养的培养。根据给出的实际案例，计量模型的重现，提高建模能力，并树立正确的风险投资意识。7.1Copula函数的定义及性质7.2Copula函数与相关性7.3常用的Copula函数7.4Copula函数的估计方法7.5Copula函数与金融风险计量7.6专题7基于GARCH-Copula模型的绿色债券投资组合风险测度目录CONTENTSCopula函数的定义及性质7.17.1.1Copula函数的定义7.1.2Copula函数的性质7.1.2Copula函数的性质7.1.3Sklar定理7.1.3Sklar定理Copula函数与相关性7.27.2Copula函数与相关性7.2Copula函数与相关性

7.2Copula函数与相关性

7.2Copula函数与相关性常用的Copula函数7.3

7.3.1椭圆类Copula函数图6.1二元正态Copula函数的分布密度图（左）及对应的等高线图（右）

R代码>library(copula)>nc=normalCopula(0.5,dim=2)>set.seed(300)>U=rCopula(1000,copula=nc)>wireframe2(nc,FUN=dCopula,col.4=adjustcolor("black",alpha.f=0.25),col="black",shade=T,delta=0.025)>contourplot2(nc,FUN=dCopula,n.grid=42,cuts=33,lwd=1/2)

7.3.1椭圆类Copula函数图6.2二元t-Copula函数的分布密度图（左）及对应的等高线图（右）

R代码>library(copula)>tc=tCopula(0.5,dim=2)>set.seed(300)>U=rCopula(1000,copula=tc)>wireframe2(tc,FUN=dCopula,col.4=adjustcolor("black",alpha.f=0.25),col="black",shade=T,delta=0.025)>contourplot2(tc,FUN=dCopula,n.grid=42,cuts=33,lwd=1/2)

7.3.2Archimedean类Copula函数

7.3.2Archimedean类Copula函数

7.3.2Archimedean类Copula函数图6.3二元GumbelCopula函数的分布密度图（左）及对应的等高线图（右）

R代码>library(copula)>gc=gumbelCopula(5,dim=2)>set.seed(300)>U=rCopula(1000,copula=gc)>wireframe2(gc,FUN=dCopula,col.4=adjustcolor("black",alpha.f=0.25),col="black",shade=T,delta=0.025)>contourplot2(gc,FUN=dCopula,n.grid=42,cuts=33,lwd=1/2)

7.3.2Archimedean类Copula函数图6.4二元ClaytonCopula函数的分布密度图（左）及对应的等高线图（右）

R代码>library(copula)>clayc=claytonCopula(0.5,dim=2)>set.seed(300)>U=rCopula(1000,copula=clayc)>wireframe2(clayc,FUN=dCopula,col.4=adjustcolor("black",alpha.f=0.25),col="black",shade=T,delta=0.025)>contourplot2(clayc,FUN=dCopula,n.grid=42,cuts=33,lwd=1/2)

7.3.2Archimedean类Copula函数图6.5二元FrankCopula函数的分布密度图（左）及对应的等高线图（右）

R代码>library(copula)>fc=frankCopula(0.5,dim=2)>set.seed(300)>U=rCopula(1000,copula=fc)>wireframe2(fc,FUN=dCopula,col.4=adjustcolor("black",alpha.f=0.25),col="black",shade=T,delta=0.025)>contourplot2(fc,FUN=dCopula,n.grid=42,cuts=33,lwd=1/2)7.3.2衍生类Copula函数图6.6混合Copula函数的分布密度图（左）及对应的等高线图（右）

R代码>library(copula)>nc=normalCopula(0.5,dim=2))>tc=tCopula(0.5,dim=2))>clayc=claytonCopula(0.5,dim=2))>gc=gumbelCopula(5,dim=2))>fc=frankCopula(0.5,dim=2))>weights=c(0.2,0.2,0.2,0.2,0.2)>mcp=mixCopula(list(nc,tc,clayc,gc,fc),w=weights)>wireframe2(mcp,FUN=dCopula,col.4=adjustcolor("black",alpha.f=0.25),col="black",shade=T,delta=0.025)>contourplot2(mcp,FUN=dCopula,n.grid=42,cuts=33,lwd=1/2)7.3.2衍生类Copula函数7.3.2衍生类Copula函数7.3.2衍生类Copula函数7.3.2衍生类Copula函数7.3.2衍生类Copula函数图6.9五维R藤的树结构7.3.2衍生类Copula函数【例7.1】为刻画银行间的高维相依关系，我们选取了平安银行、宁波银行、浦发银行、华夏银行和民生银行（分别对应数字1,2,3,4,5）进行分析，数据区间为2011年1月4日至2021年12月31日，共计2675个样本。对五家银行的收益率采用GARCH(1,1)模型进行拟合，得到残差序列后通过ecdf函数进行概率积分变换，获得拟合Copula函数的分布数据coupladata。接下来，由VineCopula包中的RVineStructrueSelect函数寻找最优藤结构。通过Summary函数可直接输出R藤的全部结果，由AIC准则选取为R藤Copula，表6.2列示了估计结果。7.3.2衍生类Copula函数6.3.2衍生类Copula函数由表7.2可看出R藤Copula函数中每个节点间的Copula函数，以及对应的估计参数、Kendall秩相关系数和上下尾相关系数；其次，考虑条件藤Copula结构，分析其条件相关系数，也就是考虑两家银行的间接相依性；最后，我们可看出由第三层至第四层的高维Copula结构的Kendall秩相关系数呈下降趋势。

R代码#载入R包>library(xts)>library(rugarch)>library(VineCopula)>library(copula)#载入数据>data=read.csv("E://jrjl/Chapter5/vinecopula.csv")>DATE=data[,1]>date=as.Date(DATE)>data=xts(data[,-1],as.Date(date,format="yyyy%mm%dd"))>PA=data$平安银行_ret>NB=data$宁波银行_ret>PF=data$浦发银行_ret>HX=data$华夏银行_ret>MS=data$民生银行_ret#边际分布拟合>spec_PA=ugarchspec(mean.model=list(armaOrder=c(1,0)),variance.model=list(model="sGARCH",garchOrder=c(1,1)),distribution.model="sstd")>garch_PA=ugarchfit(spec=spec_PA,data=PA)>spec_NB=ugarchspec(mean.model=list(armaOrder=c(1,0)),variance.model=list(model="sGARCH",garchOrder=c(1,1)),distribution.model="sstd")>garch_NB=ugarchfit(spec=spec_NB,data=NB)>spec_PF=ugarchspec(mean.model=list(armaOrder=c(1,0)),variance.model=list(model="sGARCH",garchOrder=c(1,1)),distribution.model="sstd")>garch_PF=ugarchfit(spec=spec_PF,data=PF)>spec_HX=ugarchspec(mean.model=list(armaOrder=c(1,0)),variance.model=list(model="sGARCH",garchOrder=c(1,1)),distribution.model="sstd")

R代码>garch_PF=ugarchfit(spec=spec_PF,data=PF)>spec_HX=ugarchspec(mean.model=list(armaOrder=c(1,0)),variance.model=list(model="sGARCH",garchOrder=c(1,1)),distribution.model="sstd")>garch_HX=ugarchfit(spec=spec_HX,data=HX)>spec_MS=ugarchspec(mean.model=list(armaOrder=c(1,0)),variance.model=list(model="sGARCH",garchOrder=c(1,1)),distribution.model="sstd")>garch_MS=ugarchfit(spec=spec_MS,data=MS)#数据转换>sigma_matrix=matrix(data=c(garch_PA@fit$sigma,garch_NB@fit$sigma,garch_PF@fit$sigma,garch_HX@fit$sigma,garch_MS@fit$sigma),nrow=length(garch_PA@fit$sigma),ncol=5,byrow=FALSE)>residual_matrix=matrix(data=c(garch_PA@fit$residuals,garch_NB@fit$residuals,garch_PF@fit$residuals,garch_HX@fit$residuals,garch_MS@fit$residuals),nrow=length(garch_PA@fit$residuals),ncol=5,byrow=FALSE)>std_sigma_matrix=matrix(nrow=2675,ncol=5)>copuladata=matrix(nrow=2675,ncol=5)>for(iinc(1:5)){std_sigma_matrix[,i]=residual_matrix[,i]/sigma_matrix[,i]f=ecdf(as.numeric(std_sigma_matrix[,i]))copuladata[,i]=f(std_sigma_matrix[,i])}#寻找最优的藤结构>Rst=RVineStructureSelect(copuladata,family=c(1:6),progress=TRUE,se=TRUE,method='itau',rotations=TRUE)>summary(Rst)Copula函数的估计方法7.47.4.1经验Copula7.4.1经验Copula7.4.1经验Copula【例7.2】运用经验Copula进行非参数估计。我们仍以ClaytonCopula函数为例，设定一个样本容量为n的2维样本数据，并设定ClaytonCopula函数的参数。我们运用R语言copula包中的claytonCopula函数生成所需数据，并采用copula包中的C.n()函数对样本数据进行经验估计。

R代码>library(copula)>d=2>cc=claytonCopula(3,dim=d)>n=10000>set.seed(123)>U=rCopula(n,copula=cc)>v=matrix(runif(n*d),nrow=n,ncol=d)>ec=C.n(v,X=U)>True=pCopula(v,copula=cc)>error=round(mean(abs(True-ec)/True)*100,2)>error[1]0.26

7.4.2参数估计法

7.4.2参数估计法

7.4.2参数估计法

7.4.2参数估计法

7.4.2参数估计法

R代码>cc=claytonCopula(3,dim=2)>mcc=mvdc(cc,margins=c("norm","norm"),paramMargins=list(list(mean=0,sd=1),ist(mean=0,sd=2)))>set.seed(123)>n=1000>X=rMvdc(n,mvdc=mcc)>mle=fitMvdc(X,mvdc=mcc,start=c(0,1,0,2,2))>summary(mle)

7.4.2参数估计法

7.4.2参数估计法

【例6.4】运用两阶段极大似然估计拟合Copula函数。在这里我们仍然使用例6.3中的数据，首先估计两个序列的边际分布，得到相应的参数值，再通过fitCopula()函数选ClaytonCopula，其中method=“ml”，进而得到Copula函数的参数估计值。Copula函数与金融风险计量7.57.5Copula函数与金融风险计量7.5Copula函数与金融风险计量7.5Copula函数与金融风险计量

7.5Copula函数与金融风险计量R代码##R包加载##>library(rugarch);library(mistr);library(VineCopula);library(copula)##数据导入##>data=read.csv("E://jrjl/Chapter5/shuju.csv")>DATE=data[,1]>date=as.Date(DATE)>data=xts(data[,-1],as.Date(date,format="yyyy%mm%dd"))>PA=data$PA_ret>ZS=data$ZS_ret##边际分布拟合##>garchspec_PA=ugarchspec(mean.model=list(armaOrder=c(0,0)),variance.model=list(model="sGARCH",garchOrder=c(1,1)),distribution.model="sstd")>garchfit_PA=ugarchfit(data=PA,spec=garchspec_PA)>garchspec_ZS=ugarchspec(mean.model=list(armaOrder=c(0,0)),variance.model=list(model="sGARCH",garchOrder=c(1,1)),distribution.model="sstd")>garchfit_ZS=ugarchfit(data=ZS,spec=garchspec_ZS)#提取标准化残差>standardize_residual_PA=residuals(garchfit_PA)/sigma(garchfit_PA)>standardize_residual_ZS=residuals(garchfit_ZS)/sigma(garchfit_ZS)#数据转换>PIT_PA=pdist("sstd",standardize_residual_PA,mu=0,sigma=1,skew=+coef(garchfit_PA)["skew"],shape=coef(garchfit_PA)["shape"])>PIT_ZS=pdist("sstd",standardize_residual_ZS,mu=0,sigma=1,skew=+coef(garchfit_ZS)["skew"],shape=coef(garchfit_ZS)["shape"])##选择最优Copula函数##>cop_select=BiCopSelect(PIT_PA,PIT_ZS,familyset=0:10)>summary(cop_select)7.5Copula函数与金融风险计量

R代码#抽取样本>T_model=tCopula(coef(fit_tStudent_PA_ZS)[1],dim=2,df=coef(fit_tStudent_PA_ZS)[2])>set.seed(123)>PA_ZS_tCopula_est=rCopula(2675,copula=T_model)#模拟平安银行的收益率>inverse_PA=qdist("sstd",PA_ZS_tCopula_est[,1],mu=0,sigma=1,skew=coef(garchfit_PA)["skew"],shape=coef(garchfit_PA)["shape"])>SR_PA_T=xts(x=inverse_PA,order.by=index(standardize_residual_PA))>simulate_PA=inverse_PA*coredata(sigma(garchfit_PA))>simulate_log_return_PA=simulate_PA+fitted(garchfit_PA)>plot(simulate_log_return_PA)#模拟招商银行的收益率>inverse_ZS=qdist("sstd",PA_ZS_tCopula_est[,2],mu=0,sigma=1,skew=coef(garchfit_ZS)["skew"],shape=coef(garchfit_ZS)["shape"])>SR_ZS_T=xts(x=inverse_ZS,order.by=index(standardize_residual_ZS))>simulate_ZS=inverse_ZS*coredata(sigma(garchfit_ZS))

R代码>simulate_log_return_ZS=simulate_ZS+fitted(garchfit_ZS)#构建投资组合计算VaR和ES>PA_ZS_tCopula_est=cbind(simulate_log_return_PA,simulate_log_return_ZS)>Rpa=PA_ZS_tCopula_est[,1]>Rzs=PA_ZS_tCopula_est[,2]>weight=c(0.7,0.3)>port_return_tCopula=weight[1]*Rpa+weight[2]*Rzs>VaR.95.est=quantile(port_return_tCopula,0.95)>VaR.99.est=quantile(port_return_tCopula,0.99)>ES.95.est=mean(port_return_tCopula[port_return_tCopula>VaR.95.est])>ES.99.est=mean(port_return_tCopula[port_return_tCopula>VaR.99.est])>ES.95.est[1]4.682247>ES.99.est[1]8.4976227.5Copula函数与金融风险计量表6.6投资组合的VaR值和ES值95%VaR95%ES99%VaR99%ES风险值2.78074.68225.68758.4976

从tCopula中进行2675次抽样，并考虑一个权重为0.7和0.3的投资组合，分别投资于平安银行和招商银行，模拟得到投资组合收益率，并通过式（6.25）和（6.26）计算VaR和ES。表6.6给出了该投资组合的不同置信水平的VaR值和ES值。专题7基于GARCH-Copula模型的绿色债券投资组合风险测度

7.6基于GARCH-Copula模型的绿色债券投资组合风险测度党的二十大报告对“推动绿色发展，促进人与自然和谐共生”作出战略部署，提出“必须牢固树立和践行绿水青山就是金山银山的理念，站在人与自然和谐共生的高度谋划发展”。习近平总书记指出：“推动经济社会发展绿色化、低碳化，推动经济实现质的有效提升和量的合理增长”。2023年10月召开的中央金融工作会议，也明确提出做好绿色金融这篇文章。发展绿色金融是推动实现绿色发展的必然要求，也是推动经济实现质的有效提升和量的合理增长的关键。作为绿色金融的重要融资渠道，绿色债券具备债券和绿色发展的良好属性，吸引了许多传统债券市场中关注社会责任的投资者。为了避免绿色债券可能存在的不确定性风险，一些投资者选择将绿色债券与传统债券进行组合投资。投资有风险，但应尽可能地将风险最小化。因此选择合适的分析工具来探讨金融资产间相依结构以及有效测度资产组合的风险水平，已成为亟需解决的重要问题。为了对金融风险进行预测并准确刻画资产组合的风险水平，本专题我们选择使用Copula模型和VaR模型来分析资产组合联合分布的相关问题。另外，金融资产收益率往往具有“尖峰厚尾”特征，我们通过GARCH(1,1)-skewt模型拟合各收益率的边际分布。

7.6基于GARCH-Copula模型的绿色债券投资组合风险测度1.数据来源本专题将研究中国绿色债券与传统债券投资组合风险这一问题。数据选取绿色债券、企业债和公司债三部分，并分别选取“中债-中国绿色债券财富（总值）指数”、“中债-企业债财富（总值）指数”和“中债-公司债财富（总值）指数”作为代表。研究数据区间为2013年1月4日至2023年6月30日，共计2625个观测值，数据来源于Wind数据库。为了更好的呈现收益率的波动特征，这里将收益率序列放大100倍进行分析，图6-10给出了绿色债券、企业债和公司债三个收益率波动图。

7.6基于GARCH-Copula模型的绿色债券投资组合风险测度图6-10三种债券收益率时序图7.6基于GARCH-Copula模型的绿色债券投资组合风险测度7.6基于GARCH-Copula模型的绿色债券投资组合风险测度7.6基于GARCH-Copula模型的绿色债券投资组合风险测度一方面，在投资组合和置信水平不变的情况下，绿色债券与公司债组合所计算的VaR和ES值在绝大多数情况下略小于绿色债券与企业债组合。另一方面，通过对比同一类型不同投资权重组合下的VaR和ES结果发现组合1_1和2_1的VaR和ES最小，组合1_3和组合2_3的VaR和ES最大。7.6基于GARCH-Copula模型的绿色债券投资组合风险测度本专题通过GARCH(1,1)-skewt模型拟合边缘分布，并结合Copula模型得到联合分布，对绿色债券与企业债和公司债的投资组合风险进行了定量研究。在考虑资产相依结构的前提下，模拟出资产不同权重以及不同置信度下的投资组合风险值。实证结果表明，相较于企业债而言，公司债更有助于降低绿色债券投资组合的风险，为投资者进行投资决策提供了经验。习题Theending第八章面板数据计量模型与检验学习目标掌握基础几个面板数据模型的概念、内容区分混合模型、随机效应模型和固定效应模型的使用范围与适用条件了解面板数据模型的计量程序与实际应用了解我国数字金融发展现状，掌握数字金融对区域经济发展的影响和重要性。8.1面板数据的基本界定8.2面板数据的设定和加载8.3面板回归模型8.4面板数据模型的检验8.5动态面板数据与广义矩GMM估计8.6

专题8：数字金融对地区经济发展的影响目录CONTENTS面板数据的基本界定

8.1面板数据的定义面板数据（PanelData），与时间序列数据与截面数据所区别的是在时间序列的基础上取多个截面，在这些截面上同时选取样本观测值所构成的样本数据集。或者说从数据形式上来说面板数据集是一个m*n的数据矩阵，因此具有（m，n）二维的数据性质，记载的是n个时间节点上，m个对象的某一数据指标。如表8-1所示，面板数据的格式是每个样本不同年份一个接一个叠起来的。8.1面板数据的基本界定表8-1数据型态序号年份y_pricex1_per-gdpx2_populationx3_income12014456026868326452821201543612852932951146120164456307123335601512017486633589338615371201856573601433970034120195857420613417528512020606542852330779901202162824031333280792220144811100653280562462201543601019422835950022016456310411828663987220174851908602886970622018538276218289753182201965246855529081889220206351743952718339222021618387416272849238.1面板数据的基本界定

面板数据为何有用呢？在我们进行回归分析中，有一部分受到一些观测不到的因素所影响，举个例子，地区贷款量和地区贷款政策的关系。地区贷款政策就是观察不到但是对地方贷款量求有影响，且如果对这些因素忽略，会导致遗漏变量偏差使得估计有误，但是这些因素因为观察不到，并且可能无法进行数据收集整理纳入回归。这些因子在横截面的N不同，但不会随着时间而变动，进而我们可以采用面板数据捕捉控制这些观察不到的因素对被解释变量的影响。假设纯截面数据、时间序列数据回归（地区贷款数和地区gdp）,其中yi和yt为贷款规模，xi和xt为地区gdp，同时也有一个观察不到的变量zi和zt，如我国不同地区的地方性信贷政策。8.1面板数据的基本界定如下两个方程式的回归：（8.1)式是典型横截面数据回归式，（8.2）式则是时间序列数据回归式。这两种形态的线性模型，最小二乘法的参数估计式是无偏且一致。但是，如果z是观察不到的变量，实证上就缺了这一个变量，存在遗漏变量。式(8.1)和式(8.2)使用工具变量估计法才会是无偏且一致的。但是，观察不到变量的两种情况，在使用面板数据的架构均可以解决。假设纯截面数据、时间序列数据回归（地区贷款数和地区gdp）,其中yi和yt为贷款规模，xi和xt为地区gdp，同时也有一个观察不到的变量zi和zt，如我国不同地区的地方性信贷政策。8.1面板数据的基本界定第一种情况。如果我们将(8.1)式的数据进行延展成面板数据。此时假设zit不随时间变动，故zit=zi，可得回归式(8.3）：

(8.3）将（8.3）式同步滞后一期：

（8.4）（8.3）式和（8.4）式相减：

（8.5）通过两式的差分我们可以发现无法估计的变量z在（8.5）式中抵消了，由此通过估计Δyit=β1Δxit+eit将我们想知道的β1进行无偏估计得出。8.1面板数据的基本界定第二种情况。如果我们将式（8.2）延展成面板数据。假设zit不随i变动，故zit=zt，可以使用均数移除法：移除每一个时间点的横截面平均，可得回归式（8.6）：（8.6）再进行均数处理得：（8.7）是无偏且一致。（8.6）式减去（8.7）式可得：（8.7）（8.8）面板数据的设定和加载

8.28.2

面板数据的设定和加载首先学习如何将数据导入进R。这里将数据类型分成两种情况介绍。平衡面板：对于所有的样本N，其覆盖的时间区间T都一样。非平衡面板：所有样本N中至少有一个样本的时间区间T与其他样本的时间区间T不同。我们采用收集于我国各个省或地级市的统计年鉴整理成的数据表（见下表8-2）：8.2

面板数据的设定和加载表8-2观察表头前六个样本数据

注：BalancedPanel:n=100,T=17,N=1700接下来，我们加载非平衡面板。这里我们调用r语言中自带的数据包“Hedonic”,该数据集为调查波士顿地区的自住房数量影响因素，其中变量含义如书中所示idyearprizepoplutionincomerjgdpsup120051906.178322.88139809465.008.56787120061924.866322.001598610000.989.85881120072418.383321.001977911910.4710.4429120082575.568321.002248711882.5912.8611120093237.848321.002551112968.1547.6307120104041.577316.902870814987.24342.90378.2

面板数据的设定和加载表8-3观察数据包“Hedonic”样本变量

mvcrimzninduschasnoxrmage1-110.090.01182.31no28.9443.2365.202-19.980.0307.07no22.0041.2378.902-210.450.0307.07no22.0051.6261.103-110.420.0302.18no20.9848.9745.803-210.500.0702.18no20.9851.0854.203-310.260.0302.18no20.9841.3458.70

disradtaxptratioblackslstattownidtime1-11.410.0029615.300.40-3.00112-11.600.6924217.800.40-2.39212-21.600.6924217.800.39-3.21223-11.801.1022218.700.39-3.53313-21.801.1022218.700.40-2.93323-31.801.1022218.700.39-2.96338.2

面板数据的设定和加载pdim(mydata2)#查看面板数据信息UnbalancedPanel:n=92,T=1-30,N=506进一步，将介绍如何在R中对面板数据进行描述性统计分析。Aggregate(Data,by=list(),FUN=)上面语法中有三个设定参数取第3~7列的数据为列，以N的维度当作群组因子，分组计算平均数，结果如下：aggregate(mydata1[3:7],by=list(mydata1[,"id"]),FUN="mean")8.2

面板数据的设定和加载表8-4分组计算的平均数组一prizepoplutionincomerjgdpsup14172.507326.4806244581.1924612.55534.6248724207.116270.0729453316.5679423.90148.9742534098.8351112.0154543114.8618177.231082.7583544662.147162.5679543878.6734629.50194.83064522676.1391994.0764799855.2881321.51590.1485663377.539566.5499042580.1229914.58589.1732277822.694468.8623065449.6674841.60936.130838.2

面板数据的设定和加载进一步，通过设定不同参数，我们可以得到更多的数据分析结果。aggregate(mydata1[,-c(1,2)],by=list(mydata1[,"year"]),FUN="mean")上面的代码是将数据集中的前两列数据去掉，然后依照时间T维度(year),计算平均数，这样算出来的，就是横截面面平均时间序列数据:每年都是100个城市的平均值。8.2

面板数据的设定和加载表8-5各变量每年100个城市的平均值组一yearprizepoplutionincomerjgdpsup120052508.752568.229318861.1022458.75129.3624220062847.050574.498821246.0825565.03139.0088320073526.182582.765025069.9928959.94149.6675420083896.435591.004828787.9431939.92172.2874520094555.806601.163231617.9834643.64250.3389620105523.836621.682935786.3439042.02548.75638.2

面板数据的设定和加载如果需要更多的统计衡量函数，如偏度和峰度等，可以通过加载包fBasics来实现:library(fBasics)#加载fBasics包比如：aggregate(mydata1[3:5],by=list(mydata1[,"id"]),FUN="skewness")上面我们通过fBasics实现了对第三到第五列变量的偏度8.2

面板数据的设定和加载表8-6相关变量的偏度组一prizepoplutionincome1-0.191660.531760.235622-0.17408-0.362200.0263830.37827-1.699410.245114-0.066711.138150.2949550.51973-0.735470.3760860.004350.363510.2735270.911630.589760.28759面板回归模型8.38.3.1

面板回归模型的一般形式一般来说，面板数据的回归方程式可以表示如下：用矩阵可以表示为：8.3.2

面板数据回归模型的分类在上述面板数据一般回归模型的基础上，附加上相应不同的限制性假设，使其成为不同类型的面板数据回归模型。我们可以将其分为以下几种（1）混合效应回归模型混合效应（PooledRegressionModels）估计模型：其中，8.3.2

面板数据回归模型的分类从理论上我们就可以发现，混合效应回归模型假设了解释变量对被解释变量的影响与个体以及时间无关。实际上，混合效应回归模型假设了解释变量对被解释变量的影响与个体无关。

在许多问题的研究中，由于过强的假设前提，混合效应模型有时并不适用。8.3.2

面板数据回归模型的分类（2）单因素效应模型单因素效应模型是对模型（8.10）施加了如下假设：

（8.14）在个体单因素效应模型(individualeffectregressionmodel)（8.15）

8.3.2

面板数据回归模型的分类其中，ξi与uit相互独立，并且ξi反应了个体i的非时变异质性，被称为个体效应，也可以写作下面的矩阵形式：

（8.16）其中，T是是N阶单位矩阵IN和T阶列向量lT=(1,1,1,⋯,1)'的克罗内克积。接下来考虑时间单因素效应模型(timeeffectsregressionmodel)：

（8.17）其中λt与uit相互独立，并且λt反应了个体i的时变同质性，被称为时间效应，也可以写作下面的矩阵形式：

（8.18）8.3.2

面板数据回归模型的分类同样地,如果劳动力市场、资本市场和商品市场是有效的,那么在各地区(个体)技术效率相同的假设下,根据各地区的面板数据利用柯布道格拉斯生产函数估计劳动和资本对产出的贡献时,可以将模型设定为时间单因素效应模型,前提条件允许全球技术进步是时变。因此，可以很明显得分析得出相对于混合回归模型,时间单因素效应模型更合意,它控制了不可观测的时变同质性对模型参数估计的影响。8.3.2

面板数据回归模型的分类（3）双因素效应模型所谓双因素效应模型,实际上就是在模型中既考虑了不可观测的非时变异质性效应,又考虑了不可观测时变同质性效应的面板数据线性回归模型，直观地讲就是两种单因素的结合。

（8.19）

其矩阵表示为

（8.20）8.3.2

面板数据回归模型的分类（4）固定效应模型基于单因素效应模型和双因素效应模型,如果Xit与λt是相互独立的,与ξi相关,或者Xit与λt相关,与ξi独立以及与λt、ξi都相关，这时的单因素效应模型和双因素模型都属于固定效应模型，更具体地可以分为:个体固定效应模型时间固定效应模型时间个体固定效应模型。8.3.2

面板数据回归模型的分类1）个体固定效应模型个体固定效应模型是对于不同的纵剖面时间序列（个体）只有截距项不同的模型

（8.21）或者表示为矩阵形式

（8.22）其中I_N⊗l_T是N阶单位矩阵I_N和T阶列向量l_t=(1,1,1,⋯,1)的克罗内克乘8.3.2

面板数据回归模型的分类2）时间固定效应模型时间固定效应模型就是对于不同的截面（时点）有不同截距的模型。如果确知对于不同的截面，模型的截距显著不同，但是对于不同的时间序列（个体）截距是相同的，那么应该建立时间固定效应模型（8.23）其矩阵表示为：（8.24）其中是N阶单位矩阵IN和T阶列向量lt=(1,1,1,⋯,1)的克罗内克乘，8.3.2

面板数据回归模型的分类3）时间个体双固定效应模型时间个体固定效应模型就是对于不同的截面（时点）、不同的时间序列（个体）都有不同截距的模型。表示如下:

（8.25）其矩阵表示为

（8.26）8.3.2

面板数据回归模型的分类4）固定效应模型的估计常用的固定效应的估计方法有最小二乘虚拟变量法(LeastSquareDummyVariable,LSDV)和广义最小二乘法(GeneralizedLeastSquares,GLS)两种。

LSDV法的是将每个个体的观测值表示为一个拟合常数项和一个个体固定效应的和。可以写成如下方程式：（8.27）（8.28）上式中的残差如果未知，则使用可行广义最小二乘法(FeasibleGLS。8.3.2

面板数据回归模型的分类（5）随机效应模型1）随机效应模型个体随机效应模型和时间随机效应模型，可表示为：

（8.29）

（8.30）双因素随机效应模型模型则为：

（8.31）8.3.2

面板数据回归模型的分类2）随机效应的估计在随机效应假设之下，用GLS和MLE(最大似然法)皆可以。随机效应之下的GLS和前面的不同，主要差异在随机效应需要进行基本分布假设。已知一个面板数据回归：yit=α+βxif+(μi+εit)，随机效应GLS有如下假设：E[εi]=0；E[μi]=0；E[εij

μj]=0（8.32）E[εij2]=σ_e2；E[μij2]=σμ2

（8.33）E[εijεj]=0，s≠t（8.34）E[μi

μj]=0，i≠j（8.35）8.3.2

面板数据回归模型的分类2）随机效应的估计故如同一般GLS的观念，就是一个内插逆矩阵的做法，结果如下：

（8.36）上式为：

（8.37）

且：

（8.38）8.3.2

面板数据回归模型的分类如果假设同质变异，则：

（8.39）对照Q转换，随机效应则是Ω-1转换。假设其为正态分布时，其概似函数如下：

（8.40）

（8.41）8.3.2

面板数据回归模型的分类

图1面板回归模型简易分类图面板数据模型的检验8.48.4.1

固定效应模型检验如果估计的模型设定是“固定效应”，我们就要检验用扩张的虚拟变量矩阵（LSDV）方法所估计出的个体效应，在统计是否不显著。原假设如下：

（8.42）上面的原假设中也隐含了横截面N的异质性是否在统计上不显著的问题，但它可以被忽略。标准的检验方法如F检验，概念类似于ANOVA，建立在残差平方和（RSS）的基础上：截面F=（8.43）8.4.1

固定效应模型检验除了F检验，另一个方法就是似然比（likelihoodratio，LR）检验，LR统计量在渐近上是卡方分布的，所以一般也称为卡方检验：截面（8.44）同时似然比和卡方检验结果一般可以通过观测频数和期望频数之间的差异方式进行观察8.4.2随机效应模型检验随机检验和固定效应模型检验一样，但是检验对象是由随机效应估计后的模型。如果是单维模型，原假设为：（8.45）如果是双维模型，则原假设为：（8.46）8.4.3随机效应和固定效应的选择在涉及随机效应和固定效应模型的选择问题上，随机效应较好还是固定效应较好，是一个需要检验的问题。在计量上我们使用豪斯曼检验。原假设如下：H0：(E(ui，t|Xi，t)=0（8.47）此原假设的统计量为豪斯曼统计量:

（8.48）

上式中

标符号RE代表随机效应，FE代表固定效应。根据书中结果，最终，接受原假设，选择随机效应，拒绝原假设，选择固定效应。

8.4.3随机效应和固定效应的选择【案例8.1】中国百城房价的影响因素前文为大家介绍面板数据的导入方法以及描述统计，进一步我们将介绍面板数据的混合、固定、随机效应回归以及检验的实操。通过此案例研究影响我国商品房房价的因素。在前文所展示的数据集中适配我国100个地级市的商品房均价，其中包括北京、上海、成都、重庆等一线大城市，将关键数据进一步拟合为4个变量，分别为城市常住人口、城镇居民人均年收入、实际人均gdp以及供应的住宅商品房面积，通过这四个变量来构建模型，根据不同结果得出影响房价的结论：8.4.3随机效应和固定效应的选择表8-7混合回归模型的估计结果

估计系数标准误t检验p值截距项-3.1438e+032.1625e+02-14.53742.2e-16***x12.6497e+002.1943e-0112.07542.2e-16***x21.0074e-014.4487e-032.64562.2e-16***x37.6842e-024.9929e-0315.39022.2e-16***x4-6.4606e-017.6496e-02-8.44572.2e-16***8.4.3随机效应和固定效应的选择同样如果我们假设存在影响房价的地区或时间因素，且与x相关，这是我们对方程进行固定效应模型估计。>gsp_fe0=plm(y~x1+x2+x3+x4,data=mydata1,model="within",effect="individual")#估计共同截距pool模型。利用函数内的model="within"设定固定效应model,并将估计结果存入对象gsp_fe0>summary(gsp_fe0)#对象gsp_fe0内的估计结果8.4.3随机效应和固定效应的选择表8-8固定效应模型的估计结果

估计系数标准误t检验p值x116.428020.7736421.23462.2e-16***x20.082750.0041719.85632.2e-16***x30.027440.006923.96417.691e-05***x4-0.104790.05279-1.98530.04728**8.4.3随机效应和固定效应的选择表8-9固定效应模型的估计结果对比

估计系数标准误t检验p值x116.428020.7736521.234642.4624e-88***x20.082750.0041719.856331.3760e-78***x30.027440.006923.964147.6915e-05***x4-0.104790.05279-1.985290.04728**比较两个模型，我们发现得出的估计值还是有差别，进一步介绍R的plm()函数，提供6种模型处理个体效应设定，通过其中的model参数的设置来实现model=c(“within”“random”,...)8.4.3随机效应和固定效应的选择>write.csv(summary(gsp_fe)$coef,file="table1.csv")#使用函数write,csv(),将估计系数，输出成.csv格式前文中有讲解到还可以通过可行广义最小二乘法(FeasibleGLS）来估计混合最小二乘法(pooledOLS)和固定效应。如范例程序如下。>gsp_poolFGLS=pggls(myFormula,data=mydata1,model="pooling")>summary(gsp_poolFGLS)8.4.3随机效应和固定效应的选择表8-10广义最小二乘法混合模型的估计结果

估计系数标准误z检验p值截距项-1.0428e+031.9211e+02-5.42805.700e-08***x11.3933e+002.5445e-015.47574.357e-08***x28.1421e-024.1192e-0319.76632.2e-16***x34.8315e-023.9585e-0312.20542.2e-16***x4-3.6106e-023.8751e-02-0.93170.3528.4.3随机效应和固定效应的选择>gsp_feFGLS=pggls(myFormula,data=mydata1,model="within")>summary(gsp_feFGLS)表8-11广义最小二乘法固定模型的估计结果

估计系数标准误z检验p值x18.378220.7333411.42482.2e-16***x20.073520.0038219.23792.2e-16***x30.026980.004466.04541.491e-09***X40.034960.018171.92360.054*8.4.3随机效应和固定效应的选择进一步如果我们假设存在影响房价的地区或时间因素，且不与x相关，这是我们对方程进行随机效应模型估计：>gsp_re=plm(myFormula,data=mydata1,model="random",random.method="walhus")#执行随机效应回归估计>summary(gsp_re)表8-12随机效应模型walhus的估计结果

估计系数标准误z检验p值截距项-4.4312e+034.0022e+02-11.07202.2e-16***x17.1347e+004.9506e-0114.41172.2e-16***x29.5359e-024.0822e-0323.35962.2e-16***x34.1270e-026.7478e-036.11619.591e-10***x4-1.6067e-015.6371e-02-2.85030.00437**8.4.3随机效应和固定效应的选择估计随机效应，R提供4个GLS用的权重矩阵"swar"(默认)"walhus""amemiya"和"nerlove"。"swar"就是Swamy-Arora估计式；"walhus"是Wallace-Hussain估计式，后面两个是工具变量估计随机效应模型所使用的。同时我们还可以运用极大似然法MLE对随机效应模型进行估计的具体操作如下所示。>library(nlme)>gsp_reMLE=nlme::lme(myFormula,data=house_prize,random=~1|id)#执行MLE的随机效应回归估计>summary(gsp_reMLE)8.4.3随机效应和固定效应的选择表8-13MLE随机效应模型的估计结果

估计系数标准误t检验p值截距项-7083.206658.9443-10.7493260.000***x112.4840.663918.8037700.000***x20.0890.004121.7677880.000***x30.0310.00684.5784810.000***x4-0.1190.0531-2.2435230.025**8.4.3随机效应和固定效应的选择上面的程序是对于MLE的估计，进一步来说明lme这个函数：

lme(myFormula,data=house_p