笨人工科转金融，零 stata 经验纯菜鸟，感觉陈强老师的《计量经济学及 Stata 应用》非常适合初学者，分享书中代码笔记，一起学习一起进步！
手工码字，部分 OCR，大部分都添加了注释，个别章节附有 smcl 日志文件，有错的地方直接喊我修改即可！
我用的是 Stata/MP 18.0 for Windows [64-bit x86-64]
陈强老师主页：
陈强教授的计量经济学及 Stata 主页 (econometrics-stata.com)
可以下载数据集、课件以及勘误表
我的 B 站链接 ：https://space.bilibili.com/108575181欢迎私信~
pwd //查看的工作路径，stata左下角也是显示的当前工作路径
//可直接把陈强老师的数据集复制到默认工作路径
cd “新的工作路径” //改变工作路径，但每次重启stata还会变回去
//也可以每次启动stata修改一次工作路径
sysdir // 查看各类路径 system directories
2.stata 入门
2.3.1 导入数据
use E:\grilic_small.dta,clear
2.3.3 审视数据
describe // 查看变量属性和名称 type format label
list x1 x2 //查看变量x1和x2的数据
list //查看所有变量的数据
set more off //连续滚屏显示
set more on //分页显示
list s lnw in 1/5 //只看s与lnw的前5个数据
list s lnw in 11/15 //只看s与lnw的第11一15个观测值
list s lnw if s>=16 //列出所有满足条件“s≥16”（教育年限为16年及以上）的数据
// “>=”表示“大于等于
// “==”（等于）
// “>”（大于）
// “<”（小于）
// “<="（小于等于）
// “~=”（不等于，也可用“！="表示）
// 一个等号“=”表示“赋值”，两个等号“==”表示“等于”
drop if s>=16 //删除满足“s≥16”条件的观测值
drop ln* //将所有以“ln”开头的变量都删除
keep if s>=16 //保留满足“s≥16”条件的观测值
sort s //将数据按照变量s升序排列
gsort -s //将数据按照变量s降序排列
2.3.4 画图
histogram s, width(1) frequency
// “histogram”表示直方图
// 选择项“width（1）”表示将组宽设为1（否则将使用Stata根据样本容量计算的默认分组数）
// 选择项“frequency”表示将纵坐标定为频数（默认使用密度）
// 横轴标签为变量s的label
help histogram //查看该命令的帮助文档，更多有关命令histogram的选项与用法
kdensity q //连续的经验分布图（核密度图）
scatter lnw s //绘制散点图，横轴纵轴分别为s和lnw的label
gen n =_n //定义（生成一列）变量n，表示第n个观测值
scatter lnw s, mlabel(n) mlabpos(6)
//选择项"mlabel(n)"表示:以变量n作为标签(mark label)
//选择项“mlabpos(6)”(mark label position)表示将此标签放在散点正下方(6 点钟的位置)，默认位置为散点的右边(3 点钟)
twoway (scatter tc q)(lfit tc q) //在散点图上同时画出回归直线,“lfit”表示“linear fit”(线性拟合)
graph save scatter1 //将此散点图存为文件名为“scatter1”的图像文件
twoway (scatter tc q)(qfit tc q) //在散点图上同时画出二次回归曲线,“qfit”表示“quadratic fit”(二次拟合)
graph save scatter2 //将此散点图存为文件名为“scatter2”的图像文件
graph combine scatter1.gph scatter2.gph //将上述两个图并列排放在一张图上
2.3.5 统计分析
summarize // 描述性统计（不指明变量，将显示数据集中所有变量的统计指标） Variable | Obs Mean Std. dev. Min Max
summarize s // 指定变量的描述性统计
su q if q>=10000 //子样本的统计指标
su q,detail //新增：百分位数(percentiles)，方差(variance)，偏度(skewness)与峰度(kurtosis)
tabulate s //显示变量s的经验累积分布函数(empirical cumulative distribution function)
//"Freg"表示频数,"Percent"表示百分比,"cum."表示累积百分比
pwcorr lnw s expr,sig star(.05) //显示工资对数、教育年限与工龄之间的相关系数
//"pwcorr"表示"pairwise correlation"(两两相关)
//选择项"sig"表示显示相关系数的显著性水平(即p值,列在相关系数的下方)
//选择项"star(.05)"表示给所有显著性水平小于或等于5%的相关系数打上星号
2.3.6 生成新变量
generate lns = log(s) //定义一个新变量"教育年限的对数"
gen s2 = s^2 //定义s的平方项
gen exprs = s*expr //定义s与expr的互动项(交乘项)(interaction term)
gen w = exp(lnw) //根据工资对数lnw计算工资
g colleg = (s>=16) //"虚拟变量"(dummy variable,也称"哑变量"),"s≥16"为"受过高等教育"
rename colleg college //对变量进行重命名
ren larg large //对变量进行重命名
replace college = (s>=15) //将原变量(s≥16)替换为新变量(s≥15)
2.3.7 stata 的计算器功能
display log(2) //计算ln2的值
di 2^0.5 //计算根号2的值
di normal(1.96) //计算标准正态变量小于 1.96 的概率,“normal”表示标准正态的累积分布函数
2.3.8 调用命令与终止命令
键盘"PgUp"键调用上一条命令,而"PgDn"键调用下一条命令
在历史命令窗口单击旧命令,可将旧命令调人命令窗口,然后进行编辑;双击旧命令,则将再次执行此旧命令
中途停止命令执行,可点击快捷键 Break 图标，或键盘上同时按"Ctrl+Break"
2.3.9 stata 的日志
//如果希望在每次使用Stata时,储存其运行结果,可点击菜单"File"→"Log"→"Begin"来定义"日志文件"(log file),
//Stata日志文件的扩展名为smcl
log using today //在当前路径生成一个名为"today.smcl"的日志文件Stata Markup and Control Language
log off //如果要暂时关闭日志(不再记录输出结果)
log on //恢复使用日志
log close //彻底退出日志
//查看日志文件的内容:菜单"File"→"Log"→"View'
2.4 stata、命令库的更新
update all //更新所有命令
ssc install newcommand //从SSC下载Stata命令 https://ideas.repec.org/s/boc/bocode.html
sysdir // 查看当前工作路径 system directories
//手工安装命名：将下载的新命令文件复制到PLUS所指示的那个文件夹
search keyword //搜索Stata帮助文件、Stata常见问题、Stata案例、Stata Journal,Stata Technical Bulletin等
findit keyword //搜索范围比命令search更广,还包括Stata的网络资源
//命令search的搜索结果较少,直接在Stata结果窗口显示
//命令findit的搜索结果较多,将打开另一页面显示
2.5 进一步学习 stata 的资源
Stata 常见问题 http://www.stata.com/support/faqs/
Stata 案例 http://www.stata.com/links/examples-and-datasets/
Stata Journal http://www.stata.com/bookstore/stata-journal/
Stata Technical Bulletin http://www.stata.com/products/stb/
加州大学洛杉矶分校：大量 Stata 的资源及实例 https://stats.oarc.ucla.edu/stata/
3.数学回顾
3.5 随机变量的数字特征
use grilic.dta,clear
describe
sum lnw, detail //选择项"detai1"：查看lnw的更多统计指标,比如偏度、峰度
hist lnw, width(0.1) //画lnw的直方图来看其(无条件)分布
kdensity lnw,normal normop(lpattern(dash)) //概率密度函数的连续估计
//"kdensity"表示核密度估计(kernel density estimation)
//选择项"normal"表示画正态分布的密度函数作为对比
//选择项"normop(lpattern(dash))"将正态密度用虚线(dash)来画(normop表示normal options;
//lpattern表示line pattern).如果有多条曲线，需要设置不同线型，括号用内空格隔开，如lpattern(dash dash dot)
gen wage =exp(lnw)
kdensity wage //wage的概率密度函数的连续估计
kdensity lnw if s==16 //条件分布：给定教育年限为16年(大学毕业),工资对数的条件密度
//twoway 绘制数值型变量x和y之间的关系图; scatter 散点图
twoway kdensity lnw || kdensity lnw if s ==16, lpattern(dash) //将lnw的无条件密度与条件密度画在一起
twoway (kdensity lnw)(kdensity lnw if s ==16, lpattern(dash)) //分隔符"||"(separator)可以通过两个括号"()()"来等价实现
sum lnw
sum lnw if s==16
twoway (kdensity lnw if s==12)(kdensity lnw if s ==16,lpattern(dash)) //比较在s=12(中学毕业)与s=16(大学毕业)情况下,lnw的条件密度
3.6 迭代期望定律
use grilic.dta,clear
sum lnw if rns ==0 //北方居民有554位,其条件期望E(lnw | rns=0)=5.725644
dis 5.72544*(554/(554+204))+5.581083*(204/(554+204))
sum lnw //计算无条件期望
3.8.1 正态分布(normal distribution)
//函数normalden(x)与normal(x)分别表示标准正态的密度函数中φ(x)与累积分布函数Φ(x)
dis normal(1.96) //计算标准正态变量小于1.96的概率
twoway function y=normalden(x), range(-5 5) xline(0) ytitle(概率密度) //画出标准正态的密度函数
//选择项"range(-5 5)"表示在横轴区间(-5,5)上画此图;默认为"range(0 1)",即在(0,1)区间画图
//选择项"x1ine(0)"表示在横轴x=0处画一条直线
//选择项"ytitle(概率密度)"表示将纵轴的标签设为"概率密度"
//正态分布N(m,s^2)的密度函数可用normalden(x,m,s)来表示,其中m与s分别为期望与标准差
twoway function y=normalden(x), range(-5 10) || function z=normalden(x,1,2), range(-5 10) lpattern(dash) ytitle(概率密度)
//将N(0,1)与N(1,4)的密度函数画在一起
//选择项"lpattern(dash)"表示使用虚线画图
3.8.2 卡方分布(Chi-square)
//函数chi2den(k,x)与chi2(k,x)分别表示自由度为k的卡方分布的概率密度与累积分布函数
twoway function chi3=chi2den(3,x), range(0 20) || function chi5=chi2den(5,x), range(0 20) lpattern(dash) ytitle(概率密度)
//将X^2(3)与X^2(5)的密度函数画在一起
3.8.3 t 分布
//函数tden(k,t)与t(k,t)分别表示自由度为k的t分布的概率密度与累积分布函数
twoway function t1=tden(1,x), range(-5 5) || function t5=tden(5,x), range(-5 5) lpattern(dash) ytitle(概率密度)
//将t(1)与t(5)的密度函数画在一起
//函数ttail(k,t)表示自由度为k的t分布的右侧尾部概率,即P(T>t)
dis ttail(1,95)
3.8.4 F 分布
//函数Fden(k1,k2,x)与F(k1,k2,x)分别表示自由度为(k1,k2)的F分布的概率密度与累积分布函数
twoway function F20=Fden(10,20,x),range(0 5) || function F5=Fden(10,5,x), range(0 5) lpattern(dash) ytitle(概率密度)
//将F(10,20)与F(10,5)的密度函数画在一起
//更多有关概率分布的Stata函数,参见"help density function"
4.一元线性回归
4.1 一元线性回归模型
use grilic.dta,clear
list s lnw in 1/10
twoway scatter lnw s || lfit lnw s
//scatter 散点图，x1为横坐标，x2为纵坐标
//lfit：linear fit 线性拟合
4.7 一元回归的 stata 实例
regress y x,noconstant
//"y"为被解释变量,"x"为解释变量,选择项"noconstant"表示无常数项(默认有常数项)
use grilic.dta,clear
reg lnw s
//Coef."表示回归系数(Coefficient),"_cons"表示常数项(constant)
//TSS(Tota1)
//ESS(Mode1)可解释部分
//RSS(Residua1)不可解释部分
//R2(R-squared)
reg lnw s, noc //无常数项回归，选择项noconstant
vce //显示估计系数的协方差矩阵variance covariance matrix estimated
reg lntc lnq lnpl lnpk lnpf if q>=6000 //只对“大企业”这个子样本进行回归
reg lntc lnq lnpl lnpk lnpf if large //只对虚拟变量“large=1”的子样本进行回归
reg lntc lnq lnpl lnpk lnpf if large==0 //对“小企业”(除了“大企业”以外的所有企业)进行回归
reg lntc lnq lnpl lnpk lnpf if ~large //同上，“~”表示逻辑的“否”(not)运算
4.8 Stata 命令运行结果的存储与调用
//e-类命令(e-class commands):e-类命令为"估计命令"(estimation commands),比如"regress"
//e-类命令的运行结果都存储在"e()",可以通过输入命令"ereturn list"来显示
//r-类命令(r-class commands):而所有其他命令为r-类命令,比如,"summarize".
//r-类命令的运行结果都存储在"r()",可以通过输入命令"return list"来显示
use grilic.dta,clear
sum s
return list
display r(sd)/r(mean) //计算"变异系数"(coefficient of variation,即标准差除以平均值)
reg lnw s
ereturn list
//标量(scalars)、宏(macros)、矩阵(matrices,即系数矩阵e(b)与协方差矩阵e(V)),函数(functions)
//①"宏"是Stata编程使用的一种缩写方式,它以一个简洁的字符串来代指另一个更为复杂的字符串.
//②有时,我们可能只用数据集中的一个子样本进行估计.这里的函数"e(sample)"为虚拟变量,即如果观测值在样本中,则取值为1;反之,则取值为0.
Stored results
ereturn post stores the following in e():
Scalars
e(N) number of observations
e(df_r) degrees of freedom, if specified
Macros
e(wtype) weight type
e(wexp) weight expression
e(properties) estimation properties; typically b V
Matrices
e(b) coefficient vector
e(Cns) constraints matrix
e(V) variance-covariance matrix of the estimators
Functions
e(sample) marks estimation sample
ereturn repost stores the following in e():
Macros
e(wtype) weight type
e(wexp) weight expression
e(properties) estimation properties; typically b V
Matrices
e(b) coefficient vector
e(Cns) constraints matrix
e(V) variance-covariance matrix of the estimators
Functions
e(sample) marks estimation sample
With ereturn post, all previously stored estimation results -- e() items -- are removed. ereturn repost, however, does not remove previously stored estimation results. ereturn clear removes the current e() results.
ereturn display stores the following in r():
Scalars
r(level) confidence level of confidence intervals
Macros
r(label#) label on the #th coefficient, such as (base), (omitted), (empty), or (constrained)
r(table) information from the coefficient table (see below)
r(table) contains the following information for each coefficient:
b coefficient value
se standard error
t/z test statistic for coefficient
pvalue observed significance level for t/z
ll lower limit of confidence interval
ul upper limit of confidence interval
df degrees of freedom associated with coefficient
crit critical value associated with t/z
eform indicator for exponentiated coefficients
4.9 总体回归函数与样本回归函数:蒙特卡罗模拟
clear //(删除内存中已有数据)
set obs 30 //(确定随机抽样的样本容量为30)
set seed 10101 //(指定随机抽样的"种子"为10101)
gen x=rnormal(3,2) //(得到服从N(3,2)分布的随机样本,记为x)
gen e=rnormal(0,3) //(得到服从N(0,3)分布的随机样本,记为e)
gen y= 1 + 2*x + e //(计算被解释变量y)
reg y x //(把y对x进行OLS回归)
twoway function PRF=1 +2*x, range(-5 15) || scatter y x || lfit y x, lpattern(dash)
//选择项"range(-515)"用于指定画图的横轴范围介于-5与15之间;默认为0与1之间,即"range(0 1)"
//选择项"lpattern(dash)"表示画虚线,默认画实线.
附录 A4.2 随机数的产生
set seed 10101 //(确定"种子"为10101)
set obs 30 //(确定样本容量为30)
gen x=runiform() //(产生均匀分布的随机变量x)
gen x=rnormal() //(产生标准正态分布N(0,1)的随机样本)
gen x=rnormal(m,s) //(产生正态分布N(m,s^2)的随机样本)》
gen x=rt(m) //(产生自由度为m的t分布随机样本)
gen x=rchi2(m) //(产生自由度为m的卡方分布随机样本)
5.多元线性回归
5.1 二元线性回归
use cobb_douglas.dta,clear
list
reg lny lnk lnl
predict lny1 //将被解释变量lny的拟合值记为"lny1"
predict e,residual //计算残差,并记为e.选择项"residual"表示计算残差(如果省略此选择项,则默认为计算拟合值)
list lny lny1 e //列出lny及其拟合值和残差
line lny lny1 year,lpattern(solid dash)
//将产出对数及其拟合值画在一起，纵坐标为lny和lny1，横坐标为year
//lpattern(solid dash):第一个(lny)的线型为solid，第二个(lny1)的线型为dash
5.12 多元回归的 stata 实例
//如果变量的边际效应不是常数,可考虑在回归方程中加入平方项
use grilic.dta,clear
reg lnw s expr tenure smsa rns
vce //显示回归系数的协方差矩阵"variance covariance matrix estimated".
reg lnw s expr tenure smsa rns,noc //无常数项回归
reg lnw s expr tenure smsa if rns //使用虚拟变量rns(south = 1),只对南方居民的子样本进行回归
reg lnw s expr tenure smsa if ~rns //只对北方居民的子样本进行回归，"~"表示逻辑的"否"(not)运算
reg lnw s expr tenure smsa if s >= 12 & rns //只对中学以上(S>=12)且在南方居住的子样本进行回归
quietly reg lnw s expr tenure smsa rns //前缀"quietly"表示不汇报回归结果
predict lnw1 //将lnw的拟合值记为"lnw1"
predict e,residual //计算残差,并记为e
test s =0.1 //检验原假设:变量s的系数等于0.1
//对于单个系数的（双边）检验，t分布的平方=F分布
scalar df_r=e(df_r) //save the value of df
sca list df_r //list the value of the scalar
sca coef_s=_b[s] //save the coefficient of -s-
matrix SE=e(V) //save variance-covariance matrix of the estimators
matrix list SE //list the SE matrix
sca se_s=sqrt(SE[1,1]) //if you place -s- the 2nd after the dependent variable, then change SE[1,1] to SE[2,2]
sca list coef_s //here list the coefficient of -s-
sca list se_s //here list the SE of -s-
//t=(估计量-假想值)/估计量的标准误
sca t_s=(coef_s-0.1)/se_s //手工计算t统计量
dis ttail(df_r,t_s)*2 //手工计算t统计量对应的P值,双边检验所以最后要乘2
dis ttail(e(df_r),(_b[s]-0.1)/sqrt(e(V)[1,1]))*2 //上面一堆等价于这条命令,dis可以直接调用return中的结果
test expr=tenure //检验expr与tenure的系数是否相等
test expr+tenure=s //检验工龄回报率与现单位年限回报率之和是否等于教育回报率
6.大样本 OLS
6.1 为何需要大样本理论
use grilic.dta,clear
gen wage =exp(lnw)
twoway kdensity wage, xaxis(1) yaxis(1) xvarlab(wage) || kdensity lnw, xaxis(2) yaxis(2) xvarlab(ln(wage)) lpattern(dash)
//选择项"xaxis(1) yaxis(1)"与"xaxis(2) yaxis(2)"指定对于变量wage与lnw分别使用不同的x轴与y轴
//选择项"xvarlab(wage)"与"xvarlab(ln(wage))"将变量wage与lnw核密度图的横轴标签分别指定为"wage"与"ln(wage)".
gen lns =log(s)
twoway kdensity s,xaxis(1) yaxis(1) xvarlab(s) || kdensity lns,xaxis(2) yaxis(2) xvarlab(lns) lpattern(dash)
6.2.4 依分布收敛
twoway function N=normal(x),range(-5 5)||function t1 =t(1,x),range(-5 5) lpattern(dash)||function t5 =t(5,x),range(-5 5) lpattern(shortdash) ytitle(累积分布函数)
//t分布依分布收敛于标准正态分布
twoway function N=normalden(x),range(-5 5)||function t1 =tden(1,x),range(-5 5) lpattern(dash)||function t5 =tden(5,x),range(-5 5) lpattern(shortdash) ytitle(累积分布函数)
//t分布的概率密度函数依分布收敛于标准正态分布
6.4 使用蒙特卡罗法模拟中心极限定理
program onesample,rclass //(定义程序onesample,并以r()形式储存结果)
drop _all //(删去内存中已有数据)
set obs 30 //(确定随机抽样的样本容量为30)
gen x=runiform() //(得到在(0,1)上均匀分布的随机样本)
sum x //(使用命令sum计算样本均值)
return scalar mean_sample=r(mean) //(将样本均值记为mean_sample)
end //(程序onesample结束)
set more off //(指定Stata输出结果连续翻页)
simulate xbar = r(mean_sample),seed(101) reps(10000) nodots: onesample
//选择项"reps(10000)"表示：命令simulate将运行"onesample"程序10000遍,并生成变量xbar来记录这10000个样本均值
//选择项"seed(101)"用来确定随机数的初始值
//选择项"nodots"表示不显示表示模拟过程的点(默认以一个点表示抽取一个样本)
hist xbar,normal
//选择项"norma1"表示画出相应的正态分布
program drop onesample //删除onesample程序
gen x=rchi2(10)
//在Stata中轻松运用program编写命令https://www.163.com/dy/article/GQPJV6A20538FCPL.html
6.5.1 一致估计量
twoway function y1=normalden(x,4,.2), range(3 8)||function y2=normalden(x,4.5,.3), range(3 8)||function y3=normalden(x,5,.5),range(3 8)||function y4=normalden(x,5.5,.7),range(3 8) xline(4)
6.10 大样本 0LS 的 Stata 实例
use nerlove.dta,clear
reg lntc lnq lnpl lnpk lnpf //普通标准误
display 1/_b[lnq] //估计规模报酬，_b[lnq]表示lnq的OLS系数估计值
test lnq=1 //检验规模报酬不变的原假设
reg lntc lnq lnpl lnpk lnpf,r //稳健标准误
test lnq=1
6.11 大样本理论的蒙特卡罗模拟
program chi2data,rclass //(定义程序chi2data,以r()形式储存结果)
drop _all //(删去内存中已有数据)
set obs 20 //(确定随机抽样的样本容量为20)
gen x=rchi2(1) //(生成服从X(1)分布的解释变量)
gen y=1+2*x+rchi2(10)-10 //(生成被解释变量)
reg y x //(线性回归)
return scalar b=_b[x] //(存储B的估计值)
end //(程序chi2data结束)
set more off //(指定Stata输出结果连续翻页)
simulate bhat =r(b),reps(10000) seed(10101) nodots:chi2data
//选择项"reps(10000)"表示通过命令simulate将程序"chi2data""模拟10000次.
sum bhat //计算均值与标准差
hist bhat,normal
set obs 100
set obs 1000
7.异方差
7.5.1 画残差图
use nerlove.dta,clear
reg lntc lnq lnpl lnpk lnpf
rvfplot //(residual-versus-fitted plot)画出残差与拟合值的散点图
rvpplot lnq //(residual-versus-predictor plot)画出残差与解释变量lnq的散点图
7.5.2 BP 检验
quietly reg lntc lnq lnpl lnpk lnpf
//前缀(prefix)"quietly"表示执行此命令,但不在Stata的结果窗口显示运行结果
estat hettest,iid //使用拟合值y_hat进行BP检验
estat hettest,iid rhs //使用所有解释变量进行BP检验
//"estat": post-estimation statistics(估计后统计量),即在完成估计后所计算的后续统计量
//"hettest": heteroskedasticity test
//选择项"iid"表示仅假定数据为iid,而无需正态假定
//选择项"rhs"表示,使用方程右边的全部解释变量进行辅助回归,默认使用拟合值y_hat进行辅助回归
//如果想指定使用某些解释变量进行辅助回归,可使用如下命令:
//estat hettest [varlist],iid
//其中,"[varlist]"为指定的变量清单;"[]"表示其中的内容可出现在命令中,也可不出现.
estat hettest lnq,iid //使用变量lnq进行BP检验
7.5.3 怀特检验
estat imtest, white //"imtest"指information matrix test(信息矩阵检验)
7.5.4 WLS
quietly reg lntc lnq lnpl lnpk lnpf
predict e1, residual //计算残差,并记为e1,residual可简写为r
gen e2=e1^2 //生成残差的平方,并记为e2
gen lne2=log(e2) //将残差平方取对数
reg lne2 lnq //假设lnσ^2为变量lnq的线性函数,进行辅助回归
//常数项不显著
reg lne2 lnq, noc //去掉常数项,重新进行辅助回归
predict lne2f //计算以上辅助回归的拟合值,并记为lne2f
gen e2f =exp(lne2f) //去掉对数后,即得到方差的估计值,并记为e2f
reg lntc lnq lnpl lnpk lnpf [aw=1/e2f] //使用方差估计值的倒数作为权重,进行WLS回归
//如果担心对条件方差函数的设定不准确,导致加权变换后的新扰动项仍有一定的异方差,可使用稳健标准误进行WLS估计
reg lntc lnq lnpl lnpk lnpf [aw=1/e2f], r
7.6 Stata 命令的批处理（日志）
* WLS for Nerlove(1963) //"*"表示不执行其后的命令,常用来作为注释
capture log close //表示如有已打开的日志文件,先将其关闭(如有打开的日志文件,无法定义新的日志文件)
log using wls_nerlove, replace //在当前路径创建名为"wls_nerlove.smcl"的日志文件(选择项replace表示可覆盖此文件的原有内容),并将Stata运行结果记录于此日志文件
set more off
use nerlove.dta,clear
reg lntc lnq lnpl lnpk lnpf
predict e1, r
gen e2=e1^2
gen lne2=log(e2)
reg lne2 lnq, noc
predict lne2f
gen e2f =exp(lne2f)
* Weighted least square regression
reg lntc lnq lnpl lnpk lnpf [aw=1/e2f]
reg lntc lnq lnpl lnpk lnpf [aw=1/e2f], r
log close
exit
//可在Stata中点击菜单File→Do,寻找并执行此文件
8.自相关
8.5.1 时间序列算子 time-series operator
tsset year //”tsset"表示"time series set".它告诉Stata,该数据集为时间序列,且时间变量为year.
//常用的时间序列算子包括滞后(lag)与差分(difference),分别以"L."与"D."来表示(可以小写)
//一阶滞后算子为"L.",即L.x(t)=x(t-1);二阶滞后算子为"L2.",即L2.x(t)=x(t-2),以此类推
reg y L.x L2.x L3.x L4.x //同时表示一阶至四阶滞后
reg y L(1/4).x //同时表示一阶至四阶滞后的简写
//一阶差分算子为"D.",即D.x(t)=△x(t)=x(t)-x(t-1);
//二阶差分算子为"D2.",即D2.x(t)=△(△x(t))=△(x(t)-x(t-1))=)=(x(t)-x(t-1))-(x(t-1)-x(t-2))=x(t)-2x(t-1)+x(t-2)(二阶差分为一阶差分的差分)
//以上时间序列算子可以混合使用.比如,"LD."表示一阶差分的滞后值,"DL."表示滞后值的一阶差分,二者实际上是等价的.因为LD.xt=L.(xt-xt-1)=xt-1-xt-2=D.xt-1=DL.xt
//有关时间序列算子的更多说明,参见"help tsvarlist"
8.5.2 画残差图
//假设在作完回归后,将残差记为e1
//predict e1, r esidual //计算残差,并记为e1,residual可简写为r
scatter e1 L.e1 //画出残差与其滞后的散点图
ac e1 //画出残差自相关图(即各阶自相关系数),"ac"表示autocorrelation(自相关)
8.5.3 BG 检验
//作完OLS回归后,可使用如下命令进行BG检验:
estat bgo dfrey, l ags(p) nom iss0
//选择项"lags(p)"用来指定BG检验的滞后阶数p,默认"lags(1)",即p=1;
//选择项"nomiss0"表示进行不添加0的BG检验,默认以0代替缺失值,即Davidson-MacKinnon的方法.
确定滞后阶数 p：
方法 1：使用 Sata 命令 ac 画自相关图时,所有落在 95% 的置信区域(以阴影表示)以外的自相关系数均显著地不等于 0.
方法 2：设定一个较大的 p 值,作回归：
e_{t}=\gamma_{1} e_{t-1}+\cdots+\gamma_{p} e_{t-p}+\delta_{2} x_{t 2}+\cdots+\delta_{K} x_{t K}+v_{t} \quad(t=p+1, \cdots, n)
然后看最后一个系数\gamma_{p}的显著性；如果\gamma_{p}不显著，则考虑滞后（p-1）期，以此类推，直至显著为止。
8.5.4 Q 检验
//假设将0LS残差记为e1，则可使用如下命令进行Q检验:
wntestq e1, lags(p)
//“wntestq”指white noise test O，因为白噪声没有自相关。
//选择项“lags（p）”用来指定滞后阶数，默认滞后阶数为min{loor(n/2)-2,40}.
//进行Q检验的另一命令是:
corrgram e1, lags(p)
//“corrgram”表示correlogram，即画自相关图。
//选择项“lags(p)”用来指定滞后阶数，而默认滞后阶数也是min{1oor(n/2)-2,40}.
8.5.5 DW 检验
estat dwatson //（作完OLS回归后）显示DW统计量
8.5.6 HAC 稳健标准误
newey y x1 x2 x3,lag(p) //进行OLS估计,但提供Newey-West标准误
//必选项"lag(p)"用来指定截断参数p,即用于计算HAC标准误的最高滞后阶数.
8.5.7 处理一阶自相关的 FGLS
prais y x1 x2 x3, corc //使用准差分法处理自相关
//选择项"corc"表示使用CO估计法,默认使用PW估计法.
use icecream.dta,clear
tsset time
twoway connect consumption time, msymbol(circle) yaxis(1) || connect temp time, msymbol(triangle) yaxis(2) //画出冰淇淋的消费量与气温的时间趋势图
//"connect"表示将观测点用线连接起来
//选择项"msymbol(circle)"与"msymbol(triangle)"分别表示点的"图标"(marker symbol)分别为圆圈与三角形
//选择项"yaxis(1)"与"yaxis(2)"指定使用不同的纵坐标,因为冰淇淋消费量与气温的取值范围很不相同
reg consumption temp price income
predict e1,r //计算残差(记为e1)
twoway scatter e1 L.e1 || lfit e1 L.e1 //画残差与其滞后值(L.e1)的散点图.
//"lfit"表示linear fit(线性拟合),即画出e1与1.e1的拟合回归线
twoway scatter e1 L2.e1 || lfit e1 L2.e1 //画残差与其二阶滞后的散点图
//残差似乎不存在二阶自相关(散点分布没有规律,且样本回归线的斜率接近于0).
ac e1 //画残差的自相关图
//横轴为滞后阶数,纵轴为残差的自相关系数,而阴影部分为置信度为95%的置信区间(区域).
//各阶自相关系数的取值均在95%的置信区间之内,故可接受各阶自相关系数为0的原假设.
//然而,一阶自相关系数已很接近置信区间的边界,故仍怀疑存在一阶自相关,而更高阶自相关则可大致忽略.
estat bgodfrey //BG检验,考察是否存在一阶自相关
estat bgodfrey,nomiss0 //不以0取代缺失值
wntestq e1 //Q检验, Prob > chi2(13) = 0.0160表明默认的滞后阶数为13阶,且可在5%水平上拒绝"无自相关"的原假设
corrgram e1 //使用命令corrgram进行Q检验,可以同时计算各阶Q统计量
//汇报了从1-13阶的自相关系数(AC),Q统计量(Q)及其相应p值(Prob>Q).其中,第13阶Q统计量及其p值与命令wntestq的结果完全相同
estat dwatson //计算DW统计量
newey consumption temp price income,lag(3) //异方差自相关稳健的HAC标准误.Newey-West估计量的滞后阶数为p=3: n^(1/4)=30^(1/4)≈2.34
newey consumption temp price income,lag(6) //将滞后阶数增大一倍，考察Newey-West标准误是否对于截断参数敏感
prais consumption temp price income,corc //FGLS:CO估计法
prais consumption temp price income,nolog //PW估计法,选择项"nolog"表示不显示迭代过程
reg consumption temp L.temp price income //自相关的存在可能是由于模型设定不正确,在解释变量中加入气温(temp)的滞后值,然后进行OLS回归
estat bgodfrey //使用BG检验判断重新设定的模型是否存在自相关
estat dwatson //计算DW统计量
本节较为复杂，提供 icecream.dta 的相关记录（鼠标挪到文件上，点击右边出现的图标可以下载）：
9.模型设定与数据问题
一篇专业水准的实证论文几乎总是需要说明,它是如何在存在遗漏变量的情况下避免遗漏变量偏差的.
9.4 解释变量个数的选择
use icecream.dta,clear
quietly reg consumption temp price income
estat ic //作完回归后,计算信息准则"ic"表示information criterion(信息准则)
qui reg consumption temp L.temp price income //加入气温的一阶滞后项(L.temp),重新进行估计
estat ic
qui reg consumption temp L.temp L2.temp price income //引入气温的二阶滯后项(L2.temp)
estat ic
reg consumption temp L.temp L2.temp price income //假设p_max=2，使用序贯t规则
//结果显示，L2.temp的系数高度不显著
reg consumption temp L.temp price income //令p=1，重新进行估计
//L.temp的系数在1%水平上显著(p值为0.006),故最终选择p=1.此结果与信息准则的结果相同.
9.5 对函数形式的检验
"Ramsey's RESET 检验"(Regression Equation Specification Error Test)：如果怀疑非线性项被遗漏,那么就把非线性项引入方程，并检验其系数是否显著。
use grilic.dta,clear
qui reg lnw s expr tenure smsa rns
estat ovtest //作完回归,进行RESET检验,(默认)使用拟合值的高次项进行RESET检验
estat ovtest,rhs //"ovtest""表示omitted variable test,遗漏高次项的后果类似于遗漏解释变量.
//选择项"xhs"表示使用解释变量的幂为非线性项,直接使用解释变量的高次项进行RESET检验;默认为使用拟合值的高次项进行RESET检验
gen expr2 = expr^2 //工资对数与工龄(expr)的关系可能存在非线性.引入工龄的平方项,记为expr2,再进行回归.
reg lnw s expr expr2 tenure smsa rns
estat ovtest,rhs
//工龄平方(expr2)在1%显著为正,但工龄本身(expr)却变得很不显著:因为二者存在多重共线性
//在5%水平上,接受"无遗漏变量"的原假设.事实上,在本例中,最重要的模型设定误差是遗漏了对个人能力的度量
9.6 多重共线性
"严格多重共线性"(strict multicollinearity)
"多重共线性"(multicollinearity)或"共线性"：近似(非严格)的多重共线性
多重共线性的主要表现是,如果将第 k 个解释变量x_{k}对其余的解释变量进行回归,所得可决系数(记为R_{k}^{2})较高。
解释变量x_{k}的"方差膨胀因子"(Variance Inflation Factor,VIF)为\mathrm{VIF}_{k} \equiv \frac{1}{1-R_{k}^{2}}
方差膨胀因子\mathrm{VIF}_{k}越大,则说明x_{k}的多重共线性问题越严重。
twoway function VIF =1/(1-x),xtitle(R2) xline(0.9,lp(dash)) yline(10,lp(dash)) xlabel(0.1(0.1)1) ylabel(10 100 200 300) //画出各解释变量的VIF函数
//选择项"xtitle(R2)"指示横轴的标题为"R2"
//选择项"xline(0.9,lp(dash))"与"yline(10,lp(dash))"分别表示在横轴0.9与纵轴10的位置画一条虚线
//选择项"xlabel(0.1(0.1)1)"表示在横轴上,从0.1至1,每隔0.1的位置给出标签
//选择项"ylabel(10 100 200 300)"则表示在纵轴上10、100、200与300的位置给出标签
use grilic.dta,clear
qui reg lnw s expr tenure smsa rns
estat vif //计算各变量的方差膨胀因子
gen s2 = s^2
reg lnw s s2 e expr tenure smsa rns //结果显示有两个VIF大，存在多重共线性
reg s2 s //将s2对s进行回归,R^2高
//多重共线性的一个可能的解决方法是将变量标准化,即减去均值,除以标准差
sum s //计算变量s的均值与标准差
gen sd=(s-r(mean))/r(sd) //将其标准化,并记标准化变量为sd
gen sd2 =sd^2
reg lnw sd sd2 expr tenure smsa rns
estat vif //计算各变量的方差膨胀因子：此时不存在多重共线性
reg sd2 sd ////将sd2对sd进行回归,R^2低
reg lnw sd expr tenure smsa rns //sd2在上面的回归中不显著,去掉sd2但保留sd,再次进行回归
//sd为标准化的变量,故sd变化一个单位,等价于s变化一个标准差
reg lnw s expr tenure smsa rns //未将变量s标准化的回归
//是否将变量、标准化,对于回归结果(回归系数、标准误)没有任何实质性影响
9.7 极端数据
"极端观测值"(outliers)、"高影响力数据"(influential data)
use nerlove.dta,clear
reg lntc lnq lnpl lnpk lnpf
replace lnq=lnq*100 if _n==1 //将第一个观测值的产量对数(1nq)乘以100: _n"表示第n个观测值,"_n==1"表示第1个观测值.
reg lntc lnq lnpl lnpk lnpf //再次进行回归,结果显示：人为制造一个极端值后,回归系数的估计值变化很大,而且所有系数都变得不显著
reg lntc lnq lnpl lnpk lnpf if _n > 1 //去人造极端值掉,再对比回归结果，恢复“正常”
定义第 i 个观测数据对回归系数的"影响力"或"杠杆作用"(leverage)为\operatorname{lev}_{i} \equiv \boldsymbol{x}_{i}^{\prime}\left(\boldsymbol{X}^{\prime} \boldsymbol{X}\right)^{-1} \boldsymbol{x}_{i}，\operatorname{lev}_{i}越大,则\hat{\boldsymbol{\beta}}-\hat{\boldsymbol{\beta}}^{(i)}的变化越大。\operatorname{lev}_{i}的均值为 K/n，如果某些数据的\operatorname{lev}_{i}比平均值 K/n 高很多，则可能对回归系数有很大影响。
use nerlove.dta,clear
qui reg lntc lnq lnpl lnpk lnpf
predict lev,leverage //作完回归后,计算所有观测数据的影响力lev,并记为变量lev(可自行命名)
sum lev
dis r(max)/r(mean) //看lev的最大值是其平均值的多少倍
gsort -lev //将观测值按变量lev的降序排列.如果使用命令"sort lev",则为按升序排列
list lev in 1/3 //列出lev取值最大的三个数据
replace lnq=lnq*100 if _n==1 ////人为制造极端值，将第一个观测值的产量对数(1nq)乘以100
qui reg lntc lnq lnpl lnpk lnpf
predict lev1,lev //计算所有观测数据的影响力lev,并记为变量lev1
sum lev1
dis r(max)/r(mean) //结果显示倍数高，存在极端值
9.8 虚拟变量
生成虚拟变量：
gen d=(year>=1978) //"()"表示对括弧内的表达式"yeax>=1978"进行逻辑判断.如果此表达式为真,则取值为1;反之,取值为0.
假设有 30 个省的名字储存于变量 province,希望为每个省设立一个虚拟变量,分别记为 prov1, prov2, …,prov30",则可使用如下 Stata 命令:
tabulate province,generate(prov)
//"tabulate"表示将变量按其取值列表;
//选择项"generate(prov)"表示根据此变量的不同取值生成以"prov"开头的虚拟变量.
//由此生成的这些虚拟变量,将按照变量province的字母顺序而排序.
reg y x1 x2 prov2-prov30 //在进行回归时可使用变量的简略写法
//此回归中，为了避免虚拟变量陷阱而略去了第一个省的虚拟变量pxov1.
9.9 经济结构变动的检验
对于时间序列而言,模型系数的稳定性(stability)是很重要的问题.
如果存在"结构变动"(structural break),但未加考虑,也是一种模型设定误差.
use consumption.dta,clear
twoway connect c y year, msymbol(circle) msymbol(triangle) //考察中国1978一2013年"居民人均消费"(c)与"人均国内总产值"(y)的年度(year)时间趋势图,纵轴为c和y,横轴为year
twoway connect c y year, msymbol(circle) msymbol(triangle) xlabel(1980(10)2010) xline(1992) //将上述命令略加改进,去掉右边空白，并在1992年处画一条垂直线
传统"邹检验"(Chow,1960)
reg c y
scalar ssr=e(rss)
//"sca1ar"表示标量,将此回归的残差平方和(e(rss))记为标量ssr.
reg c y if year < 1992 //对1992年之前的子样本进行回归
scalar ssr1=e(rss)
reg c y if year >= 1992
scalar ssr2=e(rss)
di((ssr-ssr1-ssr2)/2)/((ssr1+ssr2)/32) //计算邹检验的F统计量
使用虚拟变量法进行结构变动的检验
gen d=(year>1991) //生成虚拟变量d(对于1992年及以后,d=1;反之,d=0)
gen yd=y*d //生成虚拟变量d与人均收入y的互动项yd
reg c y d yd //引入d与yd,进行全样本OLS回归
test d yd //检验d与yd的联合显著性
//使用虚拟变量法所得F统计量也为15.39,与传统邹检验完全相同
//上述结构变化检验仅在球形扰动项(同方差、无自相关)的情况下才成立.为此,下面进行异方差与自相关的检验
qui reg c y
estat imtest,white //怀特检验，"imtest"指information matrix test(信息矩阵检验)(7.5.3)
//怀特检验的结果可在5%的水平上拒绝"同方差"的原假设.
tsset year //为了进行自相关的BG检验,首先设定变量year为时间变量
estat bgodfrey //BG检验,考察是否存在一阶自相关
//可在1%水平上强烈拒绝"无自相关"的原假设.
//此模型的扰动项存在异方差与自相关,故应使用异方差自相关稳健的标准误,通过虚拟变量法检验结构变动.
dis 36^(1/4) //计算HAC标准误的截断参数,结果表示应设为3
newey c y d yd,lag(3) //进行Newey-West回归
test d yd //检验虚拟变量d及其互动项yd的联合显著性
9.10 缺失数据与线性插值(linear interpolation)
线性插值的基本假设是变量以线性的速度均匀地变化.
如果变量 y 有指数增长趋势(比如 GDP),则应先取对数,再用 lny 进行线性插值,以避免偏差.
如果需要以原变量 y 进行回归,可将线性插值的对数值 lny 再取反对数(antilog),即计算 exp(lny).
y 对 x 的线性插值 \hat{y}:\frac{\hat{y}-y_{0}}{x-x_{0}}=\frac{y_{1}-y_{0}}{x_{1}-x_{0}}
use consumption.dta,clear
//为了演示目的,假设1980年、1990年、2000年及2010年的人均GDP数据缺失.
gen y1=y //生成缺失这些年份数据的人均GDP变量,并记为y1
replace y1=. if year==1980|year ==1990|year==2000|year==2010
ipolate y1 year,gen(y2) //直接用yl对year进行线性插值,并将结果记为y2
//"ipolate"表示interpolate,将变量y1对变量year进行线性插值,并将插值的结果记为新变量y2
//由于人均GDP有指数增长趋势,故更好的做法是,先对y1取对数,进行线性插值,再取反对数,并将结果记为y3.
gen lny1=log(y1) //先对y1取对数
ipolate lny1 year,gen(lny3) //进行线性插值
gen y3=exp(lny3) //再取反对数,并将结果记为y3
list year y y2 y3 if year==1980|year==1990|year==2000|year==2010 //对比这两种方法的效果
//直接插值的结果y2倾向于高估真实值y,而且整体估计效果不如先取对数再插值的结果y3(1980年的结果是个例外).
10.工具变量法
内生性：解释变量与扰动项相关
内生性的主要来源包括：
遗漏变量偏差
联立方程偏差(双向因果关系)(simultaneity bias)
测量误差偏差(measurement error bias)
解决内生性的主要方法之一为工具变量法
10.4 二阶段最小二乘法(Two Stage Least Square, 2SLS/TSLS)
y=\beta_{0}+\beta_{1} x_{1}+\beta_{2} x_{2}+\beta_{3} w+\varepsilon
x1 与 x2 为内生变量
w 为外生解释变量(与扰动项 ε 不相关)
假设有三个有效工具变 z1,z2,z3
2SLS 的第一阶段回归：
将两个内生解释变量(x1,x2)对所有外生变量(包括工具变量 z1,z2,z3 及外生解释变量 w)进行回归:
\begin{array}{c}x_{1}=\alpha_{0}+\alpha_{1} z_{1}+\alpha_{2} z_{2}+\alpha_{3} z_{3}+\alpha_{4} w+u \\ x_{2}=\gamma_{0}+\gamma_{1} z_{1}+\gamma_{2} z_{2}+\gamma_{3} z_{3}+\gamma_{4} w+v\end{array}
2SLS 的第二阶段回归：
拟合值分别记为\hat{x}_{1} 与 \hat{x}_{2}
y=\beta_{0}+\beta_{1} \hat{x}_{1}+\beta_{2} \hat{x}_{2}+\beta_{3} w+\xi
//2SLS的Stata命令格式为
ivregress 2sls y xl x2 (x3=z1 z2),robust first
//"y"为被解释变量
//"x1x2"为外生解释变量
//"x3"为内生解释变量
//"z1,z2"为方程外的工具变量
//ivregress 工具变量法的总命令，包含多种估计方法，如2SLS,最大似然法，GMM等
//2sls 为ivregress的子命令
//()等号左边为内生变量，右边为工具变量(加上外生变量，相当于2SLS的第一阶段回归)
//选择项"robust"表示使用异方差稳健的标准误(默认为普通标准误)
//选择项"first"表示显示第一阶段的回归结果
10.5 弱工具变量
estat firststage //作完2SLS回归后,检验弱工具变量
//经验规则(rule of thumb)是,如果此检验的F统计量大于10,则可拒绝"存在弱工具变量"的原假设,不必担心弱工具变量问题.
ivregress liml y x1 x2 (x3 = z1 z2) //对弱工具变量更不敏感的"有限信息最大似然估计法"(Limited Information Maximum Likelihood Estimation,LIML)
10.6 对工具变量外生性的过度识别检验
过度识别检验的原假设为“H0：所有工具变量都是外生的”。
如果拒绝该原假设，则认为至少某个工具变量不是外生的，与扰动项相关。
estat overid //作完2SLS估计后,进行过度识别检验
即使接受了过度识别的原假设，也并不能证明这些工具变量的外生性。
这是因为，过度识别检验成立的大前提是，该模型至少是恰好识别的。
此大前提无法检验，只能假定其成立。
比如，如果只有一个内生变量，则在进行过度识别检验时，我们隐含地假定至少有一个工具变量是外生的，然后检验所有其他工具变量的外生性。
10.7 对解释变量内生性的豪斯曼检验:究竟该用 OLS 还是 IV
豪斯曼检验(Hausman specification test)
"H0：所有解释变量均为外生变量"。如果 H0 成立，则 OLS 与工具变量法都一致
假设在 H0 成立的情况下，OLS 最有效率，故不适用异方差的情形
reg y x1 x2
estimates store ols //(存储OLS的结果,记为o1s)
ivregress 2sls y x1 (x2 = z1 z2) //(假设x2为内生变量,z1,z2为IV)
estimates store iv //(存储2SLS的结果,记为iv)
hausman iv ols,constant sigmamore //(进行传统豪斯曼检验)
//选择项"sigmamore"表示统一使用更有效率的估计量(即OLS)所对应的残差来计算扰动项方差
//选择项"constant"表示β_IV与β_OLS中都包括常数项(默认不含常数项)
在异方差的情况，应使用改进的“杜宾-吴-豪斯曼检验”（Durbin-Wu-Hausman Test，简记 DWH）
estat endogenous //2SLS估计后，进行异方差稳健的DWH检验
10.9 工具变量法的 Stata 实例
//作为参照系，首先进行OLS回归，并使用稳健标准误。
reg lnw s expr tenure rns smsa,r //作为参照系，首先进行OLS回归，并使用稳健标准误。
//关键解释变量为s（教育年限），而expr,tenure,rns,smsa为控制变量
//由于遗漏变量“能力”与教育年限正相关，故“能力”对工资的贡献也被纳入教育的贡献，因此高估了教育的回报率。
//引入智商（iq）作为能力的代理变量（proxy），再进行OLS回归。
reg lnw s iq expr tenure rns smsa,r
//由于用iq来度量能力存在测量误差，故ig是内生变量。
//使用变量（med，kww）作为iq的工具变量。
//母亲的教育年限（med）与KWW测试成绩（kww）都与iq正相关；并假设med与kww为外生。
ivregress 2sls lnw s expr tenure rns smsa (iq = med kww),r first //进行2SLS回归，使用稳健标准误，并显示第一阶段的回归结
estat overid //过度识别检验
//结果显示(p = 0.6972)，无法拒绝原假设：所有工具变量外生（与扰动项不相关）
//弱工具变量检验需计算第一阶段回归的普通（非稳健）F统计量，故首先使用普通标准误重新进行2SLS估计。
quietly ivregress 2sls lnw s expr tenure rns smsa (iq=med kww)
estat firststage
//由于检验第一阶段回归的两个工具变量系数联合显著性的F统计量为14.91，超过10，故认为不存在弱工具变量。
ivregress liml lnw s expr tenure rns smsa (iq = med kww),r //示范：使用对弱工具变量更不敏感的有限信息最大似然法(LIML)
//LIML的系数估计值与2SLS非常接近,从侧面印证了"不存在弱工具变量"
//豪斯曼检验，原假设"所有解释变量均为外生"，即不存在内生变量。
//传统的豪斯曼检验建立在同方差的前提下，回归中不使用稳健标准误（没有用选择项"r"）。
quietly reg lnw iq s expr tenure rns smsa
estimates store ols
quietly ivregress 2sls lnw s expr tenure rns smsa (iq = med kww)
estimates store iv
hausman iv ols,constant sigmamore
//p值Prob > chi2 = 0.0499，可在5%的显著性水平上拒绝"所有解释变量均为外生"的原假设，认为iq为内生变量。
estat endogenous //异方差稳健的DWH检验
//F统计量与一个chi2统计量的值都小于0.05，故认为iq为内生解释变量
//汇报结果：将以上各种估计法的系数及标准误列在同一表格中
qui reg lnw s expr tenure rns smsa,r
est sto ols_no_iq
qui reg lnw iq s expr tenure rns smsa,r
est sto ols_with_iq
qui ivregress 2sls lnw s expr tenure rns smsa (iq=med kww),r
est sto tsls //不能用数字开头
qui ivregress liml lnw s expr tenure rns smsa (iq=med kww),r
est sto liml
estimates table ols_no_iq ols_with_iq tsls liml,b se
//table 列表形式
//选择项“b”表示显示回归系数，选择项“se”表示显示标准误。
//用一颗星表示10%的显著性水平,两颗星表示5%的显著性水平,三颗星表示1%的显著性水平
estimates table ols_no_iq ols_with_iq tsls liml,star(0.1 0.05 0.01)
绝绝子非官方命令"estout"
同时显示回归系数、标准误与表示显著性的星号
ssc install estout //（ 盗版stata不可用，容易卡死 ）下载非官方命令"estout"
esttab ols_no_iq ols_with_iq tsls liml,se r2 mtitle star(* 0.1 ** 0.05 *** 0.01) //123颗星分别表示10%,5%,1%的显著性水平
esttab ols_no_iq ols_with_iq tsls liml using iv.rtf ,se r2 mtitle star(* 0.1 ** 0.05 *** 0.01) //输出到以iv命名Word文档（扩展名"rtf"表示rich text format），文档中倒数第二行根据实际情况看是否修改为Robust Standard errors in parentheses
11.二值选择模型
11.1 二值选择模型
twoway function Probit=normal(x),range(-5 5)||function Logit=exp(x)/(1+exp(x)),range(-5 5) lpattern(dash) ytitle(累积分布函数)
//标准正态分布与逻辑分布的累积分布函数
11.9 二值选择模型的 Stata 命令与实例
//二值选择模型的Stata命令为
probit y xl x2 x3,r //(Probit模型)
logit y xl x2 x3,r or //(Logit模型)
//选择项"x"表示使用稳健标准误(默认为普通标准误)
//选择项"or"表示显示几率比(odds ratio),而不显示回归系数
//完成Probit或Logt估计后,可进行预测,计算准确预测的百分比,或计算边际效应:
predict y1 //(计算发生概率的预测值,记为y1)
estat clas //(计算准确预测的百分比,clas表示classification)
margins,dydx(*) //(计算所有解释变量的平均边际效应;"*"代表所有解释变量)
margins,dydx(*) atmeans //(计算所有解释变量在样本均值处的边际效应)
margins,dydx(*) at(x1=0) //(计算所有解释变量在x1=0处的平均边际效应)
margins,dydx(x1) //(计算解释变量x1的平均边际效应)
margins,eyex(*) //(计算平均弹性,其中的两个"e"均指elasticity)
margins,eydx(*) //(计算平均半弹性,x变化一单位引起y变化百分之几)
margins,dyex(*) //(计算平均半弹性,x变化1%引起y变化几个单位)
use titanic.dta,clear
list
//假设观测值的重复次数记录于变量freg,则可通过在命令的最后加上"[fweight=freq]"来实现此加权计算或估计;其中"fweight"指"frequency weight'"(频数权重).
sum [fweight=freq] //各变量的统计特征,泰坦尼克号平均存活率为0.32
//分别计算小孩、女士以及各等舱旅客的存活率
sum survive if child [fweight=freq]
sum survive if female [fweight=freq]
sum survive if class1 [fweight=freq]
sum survive if class2 [fweight=freq]
sum survive if class3 [fweight=freq]
sum survive if class4 [fweight=freq]
//作为参照系,首先使用OLS估计线性概率模型
reg survive child female class1 class2 class3 [fweight=freq],r
//由于所有乘客分为四类（class1-class4），故只能放入三个虚拟变量（class1-class3）；而将虚拟变量class4（船员）作为参照类别，不放入回归方程。
//使用Logit进行估计:
logit survive child female class1 class2 class3 [fweight=freq],nolog //选择项"nolog"表示不显示MLE数值计算的迭代过程
logit survive child female class1 class2 class3 [fweight=freq],nolog r //使用稳健标准误进行Logt估计
//汇报几率比而非系数(logit+选择性or 等于 直接用logistic)
logit survive child female class1 class2 class3 [fweight=freq],or nolog //logit+选择性or（or：odds ratio）
logistic survive child female class1 class2 class3 [fweight=freq],nolog //直接用logistic
margins,dydx(*) //计算Logt模型的平均边际效应
margins,dydx(*) atmeans //计算在样本均值处的边际效应
estat clas //计算Logit模型准确预测的比率
predict prob //根据Logt模型的回归果,预测每位乘客的存活概率,并记为变量prob
//由此,可以计算给定某种特征旅客的生存概率:
list prob survive freq if class1==1 & child==0 & female==1 //计算Ms.Rose(头等舱、成年、女性)的存活概率
list prob survive freq if class3 ==1 & child==0 & female==0 //计算Mr.Jack(三等舱、成年、男性)的存活概率
probit survive child female class1 class2 class3 [fweight=freq],nolog //进行Probit估计
margins,dydx(*) //Probit模型的平均边际效应
estat clas //Probit模型的预测准确度
predict prob1 //使用Probit模型预测每位个体的存活概率,记为变量prob1
corr prob prob1 [fweight =freq] //考察prob1与prob(Logit模型预测结果)的相关性(本例中Probit与Logit模型对个体存活概率的预测结果相关系数高达0.9997,可以视为无差异)
附加例子
sysuse auto //系统自带数据集
logit foreign price mpg //foreign为01变量
\operatorname{Logit}\left(\frac{p}{1-p}\right)=\alpha+\beta_{1} \times x_{1}+\beta_{2} \times x_{2}+\varepsilon
Logit: foreign =\alpha+\beta_{1} \times price +\beta_{2} \times m p g+\varepsilon
\hat{\mathrm{y}}=\log \left(\frac{\hat{p}}{1-\hat{p}}\right)=-7.64811+0.000266 \times price +0.233835 \times m p g
\hat{p}=\frac{e^{\hat{y}}}{1+e^{\hat{y}}}
Likelihood function： \widehat{p_{1}} \times \widehat{p_{2}} \times \cdots \times\left(1-\widehat{p_{k+1}}\right) \times \cdots \times\left(1-\widehat{p_{n}}\right)
Log likelihood function：\ln \left(\widehat{p_{1}}\right)+\ln \left(\widehat{p_{2}}\right)+\cdots+\ln \left(\widehat{p_{k}}\right)+\ln \left(1-\widehat{p_{k+1}}\right)+\cdots+\ln \left(1-\widehat{p_{n}}\right)
Log likelihood = -36.462189 最优拟合直线所对应的对数似然函数取值
McFadden’s Pseudo R squared R^{2}=\frac{L L(\text { null })-L L(\text { model })}{L L(\text { null })}
LL（mode）：模型的对数似然函数
LL（null）：零模型的对数似然函数
零模型：没有自变量，只有截距项的 Logit 回归模型
logit foreign //零模型:没有自变量,只有截距项的Logit回归模型
R^{2}=\frac{L L(\text { null })-L L(\text { model })}{L L(\text { null })}=\frac{-45.03321-(-36.462189)}{-45.03321}=0.1903
Likelihood Ratio test 似然比检验(检验整个方程显著性)
$2 \times[L L( model )-L L( null )] \sim \chi^{2}(k) \quad k 自变量个数 $
检验模型的解释力度是否优于零模型
H_{0}: L L( model )=L L( null )
Likelihood Ratio 越大，P 值越小，越能拒绝原假设
本例中，logit foreign price mpg，两个解释变量，k=2
\begin{aligned} \operatorname{LR} \text { chi2(2) } & =2 \times[L L(\text { model })-L L(\text { null })] \\ & =2 \times[-36.462-(-45.033)] \\ & =17.142\end{aligned}
Prob > chi2 = 0.0002 ，所以应该拒绝原假设。即认为模型的解释力度在零模型的基础上有显著的提高
比较两个模型的拟合优度
模型1: foreign =\alpha+\beta_{1} \times price +\beta_{2} \times m p g+\varepsilon
模型2: foreign =\alpha+\beta_{1} \times price +\beta_{2} \times m p g+\beta_{3} \times trunk +\beta_{4} \times weight +\varepsilon
logit foreign price mpg
estimates store model1 //保存回归结果
logit foreign price mpg trunk weight
estimates store model2
lrtest model1 model2 //比较两个模型的拟合优度
LR chi2(2) = 38.64
Prob > chi2 = 0.0000
说明模型 2 它的拟合程度要显著的好于模型 1 的拟合程度
系数的解读
\hat{\mathrm{y}}=\log \left(\frac{\hat{p}}{1-\hat{p}}\right)=-7.64811+0.000266 \times price +0.233835 \times m p g
\hat{\mathrm{y}}=\log \left(\frac{\hat{p}}{1-\hat{p}}\right)=\beta_{0}+\beta_{1} \times price +\beta_{2} \times m p g
当自变量 X1 增加一个单位时，因变量\log \left(\frac{\hat{p}}{1-\hat{p}}\right)会增加 β1 个单位
\frac{\hat{p}}{1-\hat{p}}=e^{\beta_{0}+\beta_{1} \times \text { price }+\beta_{2} \times m p g}
当 price 取值为 𝑎, \frac{\hat{p}}{1-\hat{p}}=e^{\beta_{0}+\beta_{1} \times a+\beta_{2} \times m p g} \quad odds 1
当 price 取值为 𝑎+1, \frac{\hat{p}}{1-\hat{p}}=e^{\beta_{0}+\beta_{1} \times(a+1)+\beta_{2} \times m p g} \quad odds 2
通过比值比（odds ratio）进行解释：
odds2 除以 odds1: odds ratio =\frac{\left(\frac{\hat{p}}{1-\hat{p}}\right)_{o d d s 2}}{\left(\frac{\hat{p}}{1-\hat{p}}\right)_{o d d s 1}}=\frac{e^{\beta_{0}+\beta_{1} \times(a+1)+\beta_{2} \times m p g}}{e^{\beta_{0}+\beta_{1} \times a+\beta_{2} \times m p g}}=e^{\beta_{1}}
x 增加 1 个单位，预测概率的比值比（odds ratio）增加e^{\beta_{1}}个单位
通过发生比(odds)进行解释：
\left(\frac{\hat{p}}{1-\hat{p}}\right)_{\text {odds } 2}=e^{\beta_{1}} \times\left(\frac{\hat{p}}{1-\hat{p}}\right)_{\text {odds } 1}
x 增加 1 个单位，预测概率的发生比(odds)变成原发生比的e^{\beta_{1}}倍
\frac{o d d s 2-o d d s 1}{o d d s 1}=e^{\beta_{1}}-1
x 增加 1 个单位，预测概率的发生比(odds)变化的百分比是e^{\beta_{1}}-1，或发生比变化了 $100 \times\left(e^{\beta_{1}}-1\right) %$
回归系数总结“公式”
当系数大于 0 时：
自变量每增加一个单位，因变量对应的发生比增加e^{\beta}-1
当系数小于 0 时：
自变量每增加一个单位，因变量对应的发生比减少 $1-e^{\beta}$
logistic foreign price mpg //logistic 直接显示odds tario，即e^β
回归系数的标准误
衡量回归系数估计值与回归系数真实值之间的离散程度
Logit 模型，Z 统计量：\mathrm{H}_{0}: \beta=0
最大似然估计方法
Z=\frac{\hat{\beta}-0}{S E(\hat{\beta})} ~ 正态分布
而线性模型 ~ t 分布
Z 大于 1.65，在 10% 水平下显著
Z 大于 1.96，在 5% 水平下显著
Z 大于 2.58，在 1% 水平下显著
回归系数的置信区间：
回归系数的置信区间有 95% 的概率覆盖了系数真实值
12.面板数据
| 时间序列数据 | 同一个体，同一变量，不同时间 |
|---|---|
| 横截面数据 | 不同个体，同一变量，同一时间 |
| 面板数据 | 不同个体，同一变量，不同时间 |
面板数据分类
n 个个体，共有 T 期
短面板(short panel)：T 较小,而 n 较大,较为常见
长面板(long panel)：T 较大,而 n 较小
动态面板(dynamic panel)：解释变量包含被解释变量的滞后值
静态面板(static panel)：解释变量不包含被解释变量的滞后值
平衡面板(balanced panel)：每个时期在样本中的个体完全一样，每一个个体在每个时间内都没有缺失数据
非平衡面板(unbalanced panel)：每个时期在样本中的个体不完全一样，某些个体在某些时间内有缺失数据
出现缺失数据,怎么办?
1 可以补全的数据
可以用 0 补全:专利数据、研发费用
不随时间变化的数据
2 不可以补全的数据
企业资产、每年销售量、GDP…
伪面板数据：Pseudo-panels observe cohorts, i.e. stable groups of individuals, rather than individuals over time. Individual variables are replaced by their intra-cohort means.
一般为某个群体的特征值构成的面板数据（实际研究对象为群组中的每个个体）
但并不是所有某个群体的特征值构成的面板数据都是伪面板数据（比如这些群体本身就是研究个体）
轮换面板数据：每个时间所统计的个体中，都有一定量的个体被新的个体替换
一般是为了保持样本量不变
缺点：替换的个体数据，可能会存在个体异质性偏差
多样本多维时间数据：不同变量各自的时间序列数据堆叠在一起，不属于面板数据
面板数据模型的三种不可观测因素
Y_{i t}=\beta_{0}+\beta_{1} \times X_{1 i t}+\beta_{2} \times X_{2 i t}+\cdots+\beta_{n} \times X_{n i t}+\varepsilon_{i t}
不可观测因素\varepsilon_{i t}：因变量 Y 没有被自变量 X 所解释的部分
普通线性回归模型中不拆分，面板数据中拆分为：
\varepsilon_{i t}=\mu_{i}+v_{t}+u_{i t}
| \mu_{i} | 个体不可观测因素 | 不随时间变化而变化的个体因素 |
|---|---|---|
| v_{t} | 时间不可观测因素 | 不随个体变化而变化的时间因素 |
| u_{i t} | 其他随时间个体变化的不可观测因素 | 既随时间变化，也随个体变化 |
混合模型： 不存在个体和时间不可观测因素 \mu_{i}=v_{t}=0
当我们个体效应或者时间效应不全等于 0 的时候
不能使用混合模型!
可使用：
固定效应和随机效应都是针对因素回归模型的，每种因素模型同时具有固定和随机效应模型
混合模型中不存在固定效应和随机效应
固定效应模型：
Y_{i t}=\beta_{0}+\beta_{1} \times X_{1 i t}+\beta_{2} \times X_{2 i t}+\cdots+\beta_{n} \times X_{n i t}+\mu_{i}+v_{t}+u_{i t}
u_{i} 与 X 存在相关关系\quad 或者\quad v_{t} 与 X 存在相关关系
随机效应模型：
\mathrm{Y}_{i t}=\beta_{0}+\beta_{1} \times X_{1 i t}+\beta_{2} \times X_{2 i t}+\cdots+\beta_{n} \times X_{n i t}+\mu_{i}+v_{t}+u_{i t}
u_{i} 与 X 不存在相关关系\quad 并且\quad v_{t} 与 X 不存在相关关系
\operatorname{cov}\left(\mu_{i}, X\right)=0 \quad \& \quad \operatorname{cov}\left(v_{t}, X\right)=0
但扰动项会自相关，OLS 不是最有效率的
固定效应模型和随机效应模型 都假设：
自变量 X 对因变量 Y 的影响,对每个个体都是一样的（回归系数 β 对每一个个体都是恒定不变的）
变系数回归模型（上述假设不成立）：
Y_{i t}=\beta_{0}+\beta_{1 i} \times X_{1 i t}+\beta_{2 i} \times X_{2 i t}+\cdots+\beta_{n i} \times X_{n i t}+\mu_{i}+v_{t}+u_{i t}
β 有下角标 i
但大部分研究还是假设回归系数 β 是不随个体或者时间变化而变化的
总结：
混合模型
混合模型和普通 OLS 几乎一模一样，结果解读、形式也和 OLS 模型一模一样
假设 1:误差项与自变量不能存在相关关系
\mathrm{Y}_{i}=\beta_{0}+\beta_{1} \times X_{1 i}+\beta_{2} \times X_{2 i}+\cdots+\beta_{n} \times X_{n i}+u_{i}
\operatorname{cov}\left(u_{i}, X\right)=0
假设 2:回归系数 β 对不同个体和时间是固定不变的
假设 3:误差项相互独立\operatorname{cov}\left(u_{i}, u_{j}\right)=0
存在个体效应或时间效应，不能使用混合模型
固定效应
1.最小二乘虚拟变量法 Least squares dummy variable （LSDV）
个体效应（“单向固定效应”One-way FE）：
\mathrm{Y}_{i t}=\beta_{0}+\beta_{1} \times X_{1 i t}+\cdots+\beta_{n} \times X_{n i t}+\mu_{i}+u_{i t}
\mathrm{Y}_{i t}=\beta_{0}+\beta_{1} \times X_{1 i t}+\cdots+\beta_{n} \times X_{n i t}+\beta_{\mu 1} D_{\mu 1}+\beta_{\mu 2} D_{\mu 2}+\cdots+\beta_{\mu i-1} D_{\mu i-1}+u_{i t}
时间效应（“单向固定效应”One-way FE）：
\mathrm{Y}_{i t}=\beta_{0}+\beta_{1} \times X_{1 i t}+\cdots+\beta_{n} \times X_{n i t}+v_{t}+u_{i t}
Y_{i t}=\beta_{0}+\beta_{1} \times X_{1 i t}+\cdots+\beta_{n} \times X_{n i t}+\beta_{v 1} D_{v 1}+\beta_{v 2} D_{v 2}+\cdots+\beta_{v t-1} D_{v t-1}+u_{i t}
个体效应和时间效应（“双向固定效应”Two-way FE）：
\mathrm{Y}_{i t}=\beta_{0}+\beta_{1} \times X_{1 i t}+\cdots+\beta_{n} \times X_{n i t}+\mu_{i}+v_{t}+u_{i t}
\mathrm{Y}_{i t}=\beta_{0}+\beta_{1} \times X_{1 i t}+\cdots+\beta_{n} \times X_{n i t}+\beta_{\mu 1} D_{\mu 1}+\beta_{\mu 2} D_{\mu 2}+\cdots+\beta_{\mu i-1} D_{\mu i-1}+\beta_{p 1} D_{v 1}+\beta_{v 2} D_{v 2}+\cdots+\beta_{v t-1} D_{v t-1}+u_{i t}
reg y x1 x2 ... i.个体效应变量 i.时间效应变量 //i.变量：已此变量为分组依据构造01变量，可以只添加其中一个，或两个都添加
2.组内估计量法
个体效应：
\mathrm{Y}_{i t}=\beta_{0}+\beta_{1} \times X_{1 i t}+\cdots+\beta_{n} \times X_{n i t}+\mu_{i}+u_{i t}
第一步：以每个个体为一组，求平均值
所有的时间变量进行平均，所以没有时间下标 t
个体异质性\mu_{i}不随时间改变，平均值仍为\mu_{i}本身
\overline{Y_{i}}=\beta_{0}+\beta_{1} \times \overline{X_{1 i}}+\cdots+\beta_{n} \times \overline{X_{n i}} \quad+\mu_{i}+\overline{u_{i}}
第二步：两个式子相减
\left(\mathrm{Y}_{i t}-\overline{Y_{i}}\right)=\beta_{1} \times\left(X_{1 i t}-\overline{X_{1 i}}\right)+\cdots+\beta_{n} \times\left(X_{n i t}-\overline{X_{n i}}\right)+\left(u_{i t}-\overline{u_{i}}\right)
式子中已不存在个体效应
只要新扰动项u_{i t}-\overline{u_{i}}与新解释变量\mathbf{x}_{i t}-\overline{\mathbf{x}}_{i}不相关，则 OLS 一致，称为“固定效应估计量”（Fixed Effects Estimator），记为\hat{\boldsymbol{\beta}}_{\mathrm{FE}}。
\hat{\boldsymbol{\beta}}_{\mathrm{FE}}主要使用每位个体的组内离差信息，也称“组内估计量”（within estimator）。
时间效应：
\mathrm{Y}_{i t}=\beta_{0}+\beta_{1} \times X_{1 i t}+\cdots+\beta_{n} \times X_{n i t}+v_{t}+u_{i t}
以每个时间为一组，求平均值
\overline{Y_{t}}=\beta_{0}+\beta_{1} \times \overline{X_{1 t}}+\cdots+\beta_{n} \times \overline{X_{n t}} \quad+v_{t}+\overline{u_{t}}
\left(Y_{i t}-\bar{Y}_{t}\right)=\beta_{1} \times\left(X_{1 i t}-\overline{X_{1 t}}\right)+\cdots+\beta_{n} \times\left(X_{n i t}-\overline{X_{n t}}\right)+\left(u_{i t}-\overline{u_{t}}\right)
式子中已不存在时间效应
个体效应和时间效应：
①Y_{i t}=\beta_{0}+\beta_{1} \times X_{1 i t}+\cdots+\beta_{n} \times X_{n i t}+\mu_{i}+v_{t}+u_{i t}
以每个个体为一组，求平均值
②\bar{Y}_{i}=\beta_{0}+\beta_{1} \times \overline{X_{1 i}}+\cdots+\beta_{n} \times \overline{X_{n i}}+\mu_{i}+\bar{v}+\overline{u_{i}}
以每个时间为一组，求平均值
③\bar{Y}_{t}=\beta_{0}+\beta_{1} \times \overline{X_{1 t}}+\cdots+\beta_{n} \times \overline{X_{n t}}+\bar{\mu}+v_{t}+\overline{u_{t}}
将所有数据求平均值，y x1 x2 ... xn u v 都求平均值
④\bar{Y}=\beta_{0}+\beta_{1} \times \overline{X_{1}}+\cdots+\beta_{n} \times \overline{X_{n}}+\bar{\mu}+\bar{v}+\bar{u}
①-②-③+④：
\left(Y_{i t}-Y_{i}-\overline{Y_{t}}+\bar{Y}\right)=\beta_{1} \times\left(X_{1 i t}-\overline{X_{1 i}}-\overline{X_{1 t}}+\overline{X_{1}}\right)+\cdots+\beta_{n} \times\left(X_{n i t}-\overline{X_{n i}}-\overline{X_{n t}}+\overline{X_{n}}\right)+\left(u_{i t}-\overline{u_{i}}-\overline{u_{t}}+\bar{u}\right)
式子中已不包含个体效应和时间效应
3.一阶差分法
①y_{i t}=\mathbf{x}_{i t}^{\prime} \boldsymbol{\beta}+\mathbf{z}_{i}^{\prime} \boldsymbol{\delta}+u_{i}+\varepsilon_{i t}
②y_{i, t-1}=\mathbf{x}_{i, t-1}^{\prime} \boldsymbol{\beta}+\mathbf{z}_{i}^{\prime} \boldsymbol{\delta}+u_{i}+\varepsilon_{i, t-1} 一阶滞后
①-②：
y_{i t}-y_{i, t-1}=\left(\mathbf{x}_{i t}-\mathbf{x}_{i, t-1}\right)^{\prime} \boldsymbol{\beta}+\left(\varepsilon_{i t}-\varepsilon_{i, t-1}\right)
使用 OLS 即得到“一阶差分估计量”（First Differencing Estimator），记为\hat{\boldsymbol{\beta}}_{\mathrm{FD}}
一致性条件比组内估计量要弱。（动态面板中使用的多）
静态面板一般使用组内估计量法。
xtset 个体变量 时间变量 //告诉stata哪个变量是个体变量（第一个）和时间变量（第二个）
xtreg y x,fe //单因素个体效应模型,fe:fixed effect 固定效应回归
xtreg y x,fe i(时间变量) //单因素时间效应模型，没有考虑个体效应
xtreg y x i.时间变量,fe //（有缺陷的）双因素回归：xtreg控制个体效应，LSDV控制时间效应
reghdfe y x,absord(个体变量 时间变量) //双因素固定效应模型，reghdfe需要手动安装
//xtreg、LSDV、reghdf得到的系数和t值都一样
检验固定效应模型 VS 混合模型
Y_{i t}=\beta_{0}+\beta_{1} \times X_{1 i t}+\cdots+\beta_{n} \times X_{n i t}+\mu_{i}+u_{i t}
Y_{i t}=\beta_{0}+\beta_{1} \times X_{1 i t}+\cdots+\beta_{n} \times X_{n i t}+\beta_{\mu 1} D_{\mu 1}+\beta_{\mu 2} D_{\mu 2}+\cdots+\beta_{\mu i-1} D_{\mu i-1}+u_{i t}
原假设: \beta_{\mu 1}=\beta_{\mu 2}=\cdots=\beta_{\mu i-1}=0
备择假设: /\beta_{\mu 1}, \beta_{\mu 2}, \ldots, \beta_{\mu i-1}不全为 0
接受原假设（F 显著大于 0）：可以使用混合模型
拒绝原假设（P 小于 0.05）：应该使用固定效应模型
xtreg y x,fe //结果的最后一行的F test，检验个体效应是否存在
xtreg y x,fe i(时间变量) //结果的最后一行的F test，检验时间效应是否存在
xtreg y x i.时间变量,fe //结果的最后一行的F test，检验个体效应是否存在
testparm i.(时间变量) //检验所有时间的01变量是否都为0
随机效应模型
y_{i t}=\mathbf{x}_{i t}^{\prime} \boldsymbol{\beta}+\mathbf{z}_{i}^{\prime} \boldsymbol{\delta}+u_{i}+\varepsilon_{i t}
个体效应u_{i}与解释变量均不相关，故 OLS 一致。
扰动项由\left(u_{i}+\varepsilon_{i t}\right)组成，不是球型扰动项，故 OLS 不是最有效率的。
法 1：FGLS 估计
法 2：MLE 最大似然估计（扰动项服从正态分布）
法 3：组间估计量（损失较多信息，不常用）
xtreg y x,re //单因素个体效应，re：random effect
xtreg y x,re i(时间变量) //单因素时间效应
使用固定效应还是随机效应？
豪斯曼检验
传统的豪斯曼检验不适用于异方差的情形，须使用异方差稳健的豪斯曼检验
非平衡面板数据
不影响计算离差形式的组内估计量（within estimator）
固定效应模型的估计可照样进行
对于随机效应模型，非平衡面板数据也无实质影响，可照常进行 FGLS 估计
但缺失数据的原因不可以是内生的，不可以人为非随机的丢弃部分个体
大总结
第一步：检验个体效应、时间效应是否存在
检验个体因素是否存在：
xtset idcode year
xtreg Y X1 X2 X3, fe
结果最后一行F检验中的P值：P<0.05存在个体因素，P>0.05不存在个体因素
检验时间因素是否存在：
xtset idcode year
xtreg Y X1 X2 X3, fe i(year)
结果最后一行F检验中的P值：P<0.05存在时间因素，P>0.05不存在时间因素
检验个体和时间因素是否存在：
xtset idcode year
xtreg Y X1 X2 X3 i.year, fe
结果最后一行F检验中的P值：P<0.05存在个体因素，P>0.05不存在个体因素
testparm i.year
结果最后一行F检验中的P值：P<0.05存在时间因素，P>0.05不存在时间因素
只要有一个检验结果表示存在个体/时间因素，就保留
第二步：选择固定效应还是随机效应
个体单因素：
xtset idcode year
xtreg Y X1 X2 X3, fe
estimates store fixed
xtreg Y X1 X2 X3, re
estimates store random
hausman fixed random
结果最后一行P值：P<0.05 固定效应，P>0.05随机效应
时间单因素：
xtset idcode year
xtreg Y X1 X2 X3, fe i(year)
estimates store fixed
xtreg Y X1 X2 X3, re i(year)
estimates store random
hausman fixed random
结果最后一行P值：P<0.05 固定效应，P>0.05随机效应
双因素：推荐使用固定效应模型(stata中不容易实现，可以用R)
第三部：各模型 stata 指令
模型1：混合模型
reg Y X1 X2 X3
模型2：个体单因素固定效应模型
xtset idcode year
xtreg Y X1 X2 X3, fe
模型3：个体单因素随机效应模型
xtset idcode year
xtreg Y X1 X2 X3, re
模型4：时间单因素固定效应模型
xtset idcode year
xtreg Y X1 X2 X3, fe i(year)
模型5：时间单因素随机效应模型
xtset idcode year
xtreg Y X1 X2 X3, re i(year)
模型6：双因素固定效应模型
xtset idcode year
xtreg Y X1 X2 X3 i.year, fe
或者
reg Y X1 X2 X3 i.idcode i.year
模型7：双因素随机效应模型
使用R
面板模型的 Stata 实例
设定面板数据
xtset panelvar timevar
//“xtset”告诉 Stata 数据为面板数据
//“panelvar”个体变量，取值须为整数且不重复，相当于将样本中每位个体进行编号。
//“timevar”为时间变量。
//假如“panelvar”本来是字符串(比如，国家名字 country)，可用以下命令转换为数字型变量：
encode country, gen(cntry) //选择项“gen(cntry)”表示将新生成的数字型变量记为 cntry，变量 cntry 就以“1, 2, 3, …”来指代不同的国家。
//显示面板数据统计特性
xtdes //(显示面板数据的结构，是否为平衡面板)
xtsum //(显示组内、组间与整体的统计指标)
xtline varname //(对每位个体分别显示该变量的时间序列图；如希望将所有个体的时间序列图叠放在一起，可加上选择项overlay)
家庭联产承包责任制(household responsibility system)与中国农业增长
use lin_1992.dta,clear
xtset province year //设定 province 与 year 为面板(个体)变量及时间变量，结果显示为平衡的面板数据(strongly balanced)
xtdes //显示数据集的结构
xtsum ltvfo ltlan ltwlab ltpow ltfer hrs mipric1 giprice mci ngca //显示数据集各变量的统计特征
//overall 整体
//between 组间
//within 组内
xtline ltvfo //显示被解释变量 ltvfo 在 28 个省的时间趋势图
混合回归
reg y x1 x2 x3,vce(cluster id)
//“id”指用来确定每位个体的变量。
//选择项“vce(cluster id)”表示以变量 id 作为聚类变量来计算聚类稳健的标准误。
reg ltvfo ltlan ltwlab ltpow ltfer hrs mipric1 giprice mci ngca,vce(cluster province) //选择项“vce(cluster province)”表示使用以“province”为聚类变量的聚类稳健标准误。
estimates store OLS //将此结果储存，并记为“OLS”
reg ltvfo ltlan ltwlab ltpow ltfer hrs mipric1 giprice mci ngca //使用普通标准误，无“vce”选择项。系数不变，但普通标准误会低估std.err.
固定效应
**固定效应模型(组内估计量)
xtreg y x1 x2 x3,fe r
//选择项“fe”表示“fixed effects”(固定效应估计量)，默认为“re”表示“random effects”(随机效应估计量)。
//选择项“r”表示使用聚类稳健标准误；在使用命令xtreg时，Sata已经知道这是面板数据，故使用选择项“r"或“vce(cluster id)”，都是使用聚类稳健标准误。
**LSDV 法
reg y x1 x2 x3 i.id,vce(cluster id)
//“id”表示用来确定个体的变量，“i.id”表示根据变量“id”而生成的虚拟变量。
//选择项“vce(cluster id)”表示使用聚类稳健的标准误
//此处若使用选择项“r”，表示异方差稳健的标准误，并不是使用聚类稳健标准误。因为使用的是reg命令，stata认为处理的是横截面数据。
**使用组内估计量，并记其估计结果为“FE_robust”：
xtreg ltvfo ltlan ltwlab ltpow ltfer hrs mipric1 giprice mci ngca,fe r //输出结果包括常数项(_cons)，是所有个体效应ui的平均值。
estimates store FE_robust
**使用混合回归还是个体固定效应模型？
xtreg ltvfo ltlan ltwlab ltpow ltfer hrs mipric1 giprice mci ngca,fe //不加选择项“r”，普通标准误。输出结果含一个F检验,其原假设为"可以接受混合回归"。但由于未使用聚类稳健标准误，此F检验并不有效
estimates store FE //(将估计结果记为“FE”)
**LSDV 法：
reg ltvfo ltlan ltwlab ltpow ltfer hrs mipric1 giprice mci ngca i.province,vce(cluster province) //不少个体虚拟变量在5%水平上显著，可拒绝“所有个体虚拟变量都为0”的原假设，认为存在个体固定效应，不应使用混合回归。
estimates store LSDV //将估计结果记为“LSDV”
**一阶差分法(FD)
//Stata没有专门执行一阶差分法的命令，"xtserial"检验组内自相关，使用选择项“output”可附带提供一阶差分法的估计结果
//xtserial命令安装方法： https://blog.csdn.net/weixin_61790598/article/details/126108280
xtserial ltvfo ltlan ltwlab ltpow ltfer hrs mipric1 giprice mci ngca,output
estimates store FD //(将此结果记为“FD”)
**双向固定效应(Two-way FE):在固定效应模型中考虑时间效应
xtreg ltvfo ltlan ltwlab ltpow ltfer hrs mipric1 giprice mci ngca t,fe r //为节省待估参数，加入时间趋势项t(结果显示该时间趋势项并不显著，可以考虑去掉)
estimates store FE_trend //(将估计结果记为“FE_trend”)
** 展示：加入年度虚拟变量 （本方法复杂，后有简易方法）
tab year,gen(year_) //定义年度虚拟变量
//tabulate 列表，根据year的不同取值，生成一系列以“year_”开头的虚拟变量
xtreg ltvfo ltlan ltwlab ltpow ltfer hrs mci ngca year_2-year_18,fe r //进行含时间虚拟变量的双向固定效应估计
//加入年度虚拟变量后， 由于价格变量 mipric1 与 giprice 在各省都一样，无法包括在回归方程中，以避免严格多重共线性 。
estimates store FE_TW //（将结果记为“FETW”）：
//检验所有年度虚拟变量的联合显著性
test year_2 year_3 year_4 year_5 year_6 year_7 year_8 year_9 year_10 year_12 year_13 year_14 year_15 year_16 year_17 year_18 //(结果显示，强烈拒绝“无时间固定效应”的原假设，应在模型中包括时间固定效应。)
** 简易方法 ：直接用以下命令来估计双向固定效应模型(不必先生成时间虚拟变量)
xtreg ltvfo ltlan ltwlab ltpow ltfer hrs mci ngca i.year,fe r
//“i.year”表示根据变量 year 的不同取值来生成年度虚拟变量
随机效应
//随机效应估计的 Stata 命令
xtreg y x1 x2 x3,re r theta
//选择项“re”为默认选项(可省略)
//选择项“r”表示使用聚类稳健标准误。xtreg时等同于使用选择项“vce(cluster id)”
//选择项“theta”表示显示用于进行广义离差变换的θ值
//随机效应模型进行 MLE 估计
xtreg y x1 x2 x3,mle
//进行随机效应(RE)的估计
xtreg ltvfo ltlan ltwlab ltpow ltfer hrs mci ngca,re r theta
estimates store RE_robust //(将结果记为“RE_robust”)
//检验使用混合回归，还是个体随机效应模型？
//检验个体随机效应的LM 检验，原假设为不存在个体的随机效应。如拒绝原假设，则模型中应有反映个体特性的随机扰动项，不应使用混合回归。
//LM检验在"进行随机效应(RE)的估计"命令之后才能进行
xttest0 //LM 检验 (拒绝“不存在个体随机效应”的原假设( p值 0.0000)，认为应选择“随机效应”)
//展示：使用普通标准误的随机效应估计(将结果记为“RE”)，去掉选择项“r”，回归系数不变，但会低估真实标准误
xtreg ltvfo ltlan ltwlab ltpow ltfer hrs mci ngca,re
estimates store RE
//对照：对随机效应模型进行最大似然 MLE 估计
xtreg ltvfo ltlan ltwlab ltpow ltfer hrs mci ngca,mle nolog
//结果中最后一行的LR检验：个体随机效应是否存在。结果显示拒绝原假设 Prob >= chibar2 = 0.000，认为存在个体随机效应，不应进行混合回归
//随机效应 MLE 的系数估计值与 FGLS 有所不同，但在性质上依然类似
estimates store MLE //(将结果记为“MLE”)
固定效应还是随机效应：豪斯曼检验
**豪斯曼检验的 Stata 命令为
xtreg y x1 x2 x3,fe //(固定效应估计)
estimates store FE //(存储结果为FE)
xtreg y x1 x2 x3,re //(随机效应估计)
estimates store RE //(存储结果为RE)
hausman FE RE,constant sigmamore //(豪斯曼检验)
//选择项“constant”表示在比较系数估计值时包括常数项(默认不含常数项)
//选择项“sigmamore”表示统一使用更有效率的那个估计量(即随机效应估计量)的方差估计
hausman FE RE,constant sigmamore //此结果与书上不一样，因为生成FE和RE的解释变量不一致，应修改FE如下
//FE：xtreg ltvfo ltlan ltwlab ltpow ltfer hrs mipric1 giprice mci ngca,fe
xtreg ltvfo ltlan ltwlab ltpow ltfer hrs mci ngca,fe
estimates store FE_hausman
hausman FE_hausman RE,constant sigmamore //此结果与书上相同
//Prob > chi2 = 0.0000拒绝原假设，应使用固定效应模型，而非随机效应模型
**稳健的豪斯曼检验
//原假设：随机效应模型正确
//非官方命令 xtoverid,“overid”指“overidentification test”(过度识别检验)
//命令安装方式： https://www.bilibili.com/read/cv30012101/?spm_id_from=333.999.0.0
quietly xtreg ltvfo ltlan ltwlab ltpow ltfer hrs mci ngca,r
xtoverid //在使用命令 xtoverid 之前，须先以稳健标准误执行xtreg命令
//p值为 0.0000，拒绝随机效应，应使用固定效应模型
组间估计量
xtreg ltvfo ltlan ltwlab ltpow ltfer hrs mci ngca,be //"be": use between-effects estimator
//根据豪斯曼检验，应选择固定效应，而组间估计量仅在随机效应情况下才一致，故组间估计结果不可信。
将主要方法的回归系数及标准误列表
esttab OLS FE_robust FE_trend FE RE,b se mtitle
//esttab外部命令见第8章，选择项“b”表示显示回归系数，选择项“se”表示显示标准误
//mtitles[(list)]: specify model titles to appear in table header
//结果中第四列由于相同原因，与书上不一致，修改如下：
esttab OLS FE_robust FE_trend FE_hausman RE,b se mtitle //此结果与书上相同，但陈强老师勘误中指出应使用上一条命令，书上的第四列有误
//变量的度量精确度可能会影响显著性
//制度变量一般都具有内生性，使用面板的工具变量法解决内生性
13.平稳时间序列
"严格平稳过程"：要求随机过程的有限维分布不随时间推移而改变
13.1 时间序列的自相关
“平稳序列”(stationary)与“非平稳序列”(non-stationary)
k 阶自协方差(autocovariance of order k)
\gamma_{k} \equiv \operatorname{Cov}\left(y_{t}, y_{t+k}\right)=\mathrm{E}\left[\left(y_{t}-\mu\right)\left(y_{t+k}-\mu\right)\right]\quad \mu \equiv \mathrm{E}(y)
反映同一变量( y)相隔 k 期之间的自相关程度
k 阶样本自协方差
\hat{\gamma}_{k} \equiv \frac{1}{T-k} \sum_{t=1}^{T-k}\left(y_{t}-\bar{y}\right)\left(y_{t+k}-\bar{y}\right) \quad \bar{y} \equiv \frac{1}{T} \sum_{t=1}^{T} y_{t}
k 阶自相关系数(autocorrelation of order k)
\rho_{k} \equiv \operatorname{Corr}\left(y_{t}, y_{t+k}\right) \equiv \frac{\operatorname{Cov}\left(y_{t}, y_{t+k}\right)}{\operatorname{Var}\left(y_{t}\right)}
k 阶样本自相关系数
\hat{\rho}_{k} \equiv \hat{\gamma}_{k} / \hat{\gamma}_{0} \quad \hat{\gamma}_{0} \equiv \frac{1}{T-1} \sum_{t=1}^{T}\left(y_{t}-\bar{y}\right)^{2}
对于严格平稳过程，\rho_{k}不依赖于具体时间，仅是滞后阶数 k 的函数，称为“自相关函数”(Autocorrelation Function，简记 ACF)
将\left(k, \rho_{k}\right)画成图，为“自相关图”(correlogram)
解决 gdp_china.dta 汉字乱码
//先cd切换到dta所在目录
unicode analyze gdp_china.dta
unicode encoding set gb18030
unicode translate gdp_china.dta
use gdp_china.dta,clear
tsset year ////”tsset"表示"time series set".它告诉Stata,该数据集为时间序列,且时间变量为year.
tsline y,xlabel(1980(10)2010)
//“tsline”表示画时间趋势图，在此等价于命令“line y year”(year为时间变量)
//“xlabel(1980(10)2010)”表示在横轴1980-2010期间，每隔10年做个标注(label)
//GDP存在指数增长（exponential growth）的趋势
gen lny=log(y) //将GDP取对数，把指数趋势变为线性趋势
tsline lny,xlabel(1980(10)2010)
//GDP 对数存在线性趋势，但依然不平稳(期望值不断增长)
gen dlny=d.lny //将 GDP 对数进行一阶差分
tsline dlny,xlabel(1980(10)2010)
GDP 对数差分约等于 GDP 的增长率：
\begin{aligned} \Delta \ln y_{t} & \equiv \ln y_{t}-\ln y_{t-1}=\ln \left(\frac{y_{t}}{y_{t-1}}\right) \\ & =\ln \left(\frac{y_{t-1}+\Delta y_{t}}{y_{t-1}}\right)=\ln \left(1+\frac{\Delta y_{t}}{y_{t-1}}\right) \approx \frac{\Delta y_{t}}{y_{t-1}}\end{aligned}
gen g=(y-l.y)/l.y //计算 GDP 的增长率(记为g) l.y滞后算子，表示上一期的y
tsline dlny g,xlabel(1980(10)2010) lpattern(dash)
//lpattern(dash),第一条曲线为虚线，第二条曲线默认为实线
//lpattern(dash dot)，第一条曲线为虚线，第二条曲线为点线，中间空格隔开
linepatternstyle Description
----------------------------------------
solid solid line
dash dashed line
dot dotted line
dash_dot
shortdash
shortdash_dot
longdash
longdash_dot
blank invisible line
"formula" e.g., "-." or "--.." etc.
-----------------------------------------
A formula is composed of any combination of
l solid line
_ (underscore) a long dash
- (hyphen) a medium dash
. short dash (almost a dot)
# small amount of blank space
------------------------------------------
corrgram dlny //“corrgram”表示 correlogram，即画自相关图
//第1列滞后阶数(LAG)
//第2列表示自相关系数(AC)，倒数第2列即为自相关图(以线条长度表示各阶自相关系数)
//第3列PAC为偏自相关系数(Partial Autocorrelation)，倒数第1列为PAC的图示
//Q统计量(Q)及相应的p值(Prob>Q),其原假设为"从一阶开始到各阶自相关系数均为0"
ac dlny,lags(20) //画自相关图的另一命令
//“ac”表示 autocorrelation
//选择项“lags(20)”表示画 1-20 阶的自相关图；默认所画的最高阶数为min{floor(n/2)-2,40}，其中floor(n/2)为不超过n/2的最大整数。
//实心点表示各阶自相关系数
//阴影部分表示置信区域(confidence band)
//如果某阶自相关系数落在置信区域之外,则说明该阶自相关系数显著地不为0(拒绝原假设):反之,则不显著(接受该阶自相关系数为0的原假设)
“单变量时间序列”(univariate time series)：仅关心某变量的未来值，不理会因果关系，只考虑相关关系
一阶自回归模型(AR(1))：使用过去值预测当前值
y_{t}=\beta_{0}+\beta_{1} y_{t-1}+\varepsilon_{t} \quad(t=2, \cdots, T)
扰动项 \varepsilon_{t} 为白噪声, 故 \operatorname{Cov}\left(\varepsilon_{t}, \varepsilon_{s}\right)=0, \forall t \neq s
① \widehat{\Delta \ln y_{t}}=0.0437698+0.5362727 \Delta \ln y_{t-1}
② \widehat{\Delta \ln y_{2013}}=0.083309
③ \widehat{\ln y_{2013}}=\ln y_{2012}+\widehat{\Delta \ln y_{2013}}
④ 2013 年 GDP 的预测值 \widehat{y_{2013}}=\exp \left(\ln y_{2012}+\widehat{\Delta \ln y_{2013}}\right)
⑤ 预测误差=y_{2013}-\widehat{y_{2013}}
**以OLS估计lny_t的一阶自回归模型
reg dlny l.dlny if year<2013,r //使用2013年前的数据回归，1式
//假设扰动项ε_t无自相关，故使用异方差稳健的标准误即可，不必使用异方差自相关稳健的HAC标准误
predict dlny1 //计算回归方程的拟合值,并记为dlny1
list dlny1 if year==2013 //2式
dis exp(lny[35]+dlny1[36]) //4式；“x[n]”表示变量x的第n个观测值
//因为样本容量为36,"lny[35]"表示变量 lny 的第 35 个观测值(即 2012 年),而"dlny1[36]"表示变量 dlny1 的第 36 个观测值(即 2013 年)
dis y[36] //查看2013年的实际GDP
dis y[36]-exp(lny[35]+dlny1[36]) //计算预测误差，5式
13.3 高阶自回归
p 阶自回归模型，记为 AR(p)：y_{t}=\beta_{0}+\beta_{1} y_{t-1}+\cdots+\beta_{p} y_{t-p}+\varepsilon_{t}
需要估计滞后阶数 p
方法 1：“由大到小的序贯 t 规则” (general-to-specific sequential t rule)
方法 2：使用信息准则
实践中，可结合以上两种方法来确定 $ p
$，如二者结果不一致，为了保守起见(尽量避免遗漏变量偏差)，可取二者滞后阶数的大者。
可检验模型的残差是否存在自相关（比如，使用 Q 检验）：如果残差存在自相关，则须扩大滞后阶数。
**确定滞后阶数
quietly reg dlny l.dlny if year<2013,r //一阶自回归AR(1)
estat ic //计算information criterion
reg dlny l(1/2).dlny if year<2013,r //二阶自回归AR(2)，l(1/2).dlny表示1阶和2阶的滞后
estat ic
//回归结果显示dlny的二阶滞后项L2.d1ny依然在1%水平上显著,根据序贯t规则,滞后阶数P应至少大于或等于2.
//信息准则结果显示AR(2)的AIC与BIC均比AR(1)略有上升.故根据信息准则,应选择p=1,即AR(1)模型.
//两者结果不一致，保守起见，可取二者滞后阶数的大者。
reg dlny l(1/3).dlny if year<2013,r ////三阶自回归AR(3)，l(1/3).dlny表示1阶、2阶和3阶的滞后
estat ic
//回归结果显示dlny的三阶滞后项不显著很不显著，根据序贯t规则，应选择p=2
//信息准则结果显示AR(3)的AIC与BIC均比AR(2)都有上升.故根据信息准则,应选择p=2,即AR(2)模型.
//综上，应选择p=2
没有做：使用命令 corrgram 对残差进行 Q 检验也表明，AR(1)的残差存在自相关，而 AR(2)的残差无自相关(参见习题)。
**使用 AR(2)模型预测 GDP，并与 AR(1)的预测效果对比
quietly reg dlny l(1/2).dlny if year<2013,r
predict dlny2
dis exp(lny[35]+dlny2[36])
dis y[36]-exp(lny[35]+dlny2[36])
//对于 2013 年的 GDP，AR(2)模型的预测误差为-680.78688 亿元，即高估了 680.78688 亿元；
//AR(1)模型则高估了 895.90347 亿元。
//AR(2)的预测效果优于AR(1)，因为二阶滞后仍包含有用信息,AR2回归中二阶滞后项高度显著
13.4 自回归分布滞后模型(Autoregressive Distributed Lag Model)ADL
在 AR(p)模型中，为了提高预测力，也可引入(不止一个)其他解释变量
y_{t}=\beta_{0}+\underbrace{\beta_{1} y_{t-1}+\cdots+\beta_{p} y_{t-p}}_{A R(p)}+\underbrace{\gamma_{1} x_{t-1}+\cdots+\gamma_{q} x_{t-q}}_{D L(q)}+\varepsilon_{t}（只引入了一个其他解释变量）
p 为 y 的自回归阶数，而 q 为 x 的滞后阶数
假定扰动项\varepsilon_{t}为白噪声，则 OLS 一致
滞后阶数( p ,q )的选择，可使用信息准则(AIC 或 BIC)，或进行序贯检验，即使用 t 或 F 检验来检验最后一阶系数的显著性
ADL 为一个动态模型，y 有自回归，x 对 y 也有分布滞后的作用
解释变量 x_{t-1} 对于 y_{t} 的边际效应为 \gamma_{1}, 但这并非长期效应
\begin{array}{l}\text { 由于 }\left\{y_{t}\right\} \text { 与 }\left\{x_{t}\right\} \text { 为平稳序列, 故均值不随时间而变, 分别记为 } y^{*} \text { 与 } x^{*} \text { 。 }\end{array}
x^{*} 增加一单位对 y^{*} 的边际效应为\frac{d y^{*}}{d x^{*}}=\frac{\gamma_{1}+\cdots+\gamma_{q}}{1-\beta_{1}-\cdots-\beta_{p}}
这就是 x 永久性增加一单位对 y 的长期效应，也称“长期乘数”(long-run multiplier)
use border.dta, clear
tsset decade
reg border l(1/2).border l.drought diff age rival wall unified,r
**计算气候冲击对游牧边界的长期效应
dis -.6333046/(1-1.518284+.5586965)
13.5 误差修正模型(Error Correction Model，ECM)
13.6 移动平均与 ARMA 模型
另一类时间序列模型为“移动平均过程”(Moving Average Process，简记 MA)
一阶移动平均过程为 MA(1)：
y_{t}=\mu+\varepsilon_{t}+\theta \varepsilon_{t-1} （μ 为常数项，\left\{\varepsilon_{t}\right\}为白噪声，而\varepsilon_{t}的系数被标准化为 1，y 可看成是白噪声的移动平均）
q 阶移动平均过程 MA(q)：
y_{t}=\mu+\varepsilon_{t}+\theta_{1} \varepsilon_{t-1}+\theta_{2} \varepsilon_{t-2}+\cdots+\theta_{q} \varepsilon_{t-q} （假设{\varepsilon_{t}}为 iid 且服从正态分布，用 MLE 估计）
无穷阶移动平均过程 MA(∞)：
y_{t}=\mu+\varepsilon_{t}+\theta_{1} \varepsilon_{t-1}+\theta_{2} \varepsilon_{t-2}+\cdots=\mu+\sum_{j=0}^{\infty} \theta_{j} \varepsilon_{t-j} （将 y 的决定因素追溯到无穷远的过去）
MA(∞)收敛的充分条件：\left\{\theta_{j}\right\}_{j=0}^{\infty}为“绝对值可加总”(Absolutely Summable，简记 AS)，\sum_{j=0}^{\infty}\left|\theta_{j}\right|<\infty（有限）
将 AR(p)与 MA(q)结合起来，可得 ARMA(p, q)模型：
y_{t}=\beta_{0}+\underbrace{\beta_{1} y_{t-1}+\cdots+\beta_{p} y_{t-p}}_{A R(p)}+\underbrace{\varepsilon_{t}+\theta_{1} \varepsilon_{t-1}+\cdots+\theta_{q} \varepsilon_{t-q}}_{M A(q)}
13.7 脉冲响应函数
对于 y_{t}=\beta_{0}+\beta_{1} y_{t-1}+\varepsilon_{t}, 假设 \left|\beta_{1}\right|<1, 则此 \mathrm{AR}(1) 是 \mathrm{MA}(\infty)
\begin{aligned} y_{t} & =\beta_{0}+\beta_{1} y_{t-1}+\varepsilon_{t} \\ & =\beta_{0}+\beta_{1}\left(\beta_{0}+\beta_{1} y_{t-2}+\varepsilon_{t-1}\right)+\varepsilon_{t} \\ & =\left(\beta_{0}+\beta_{0} \beta_{1}\right)+\beta_{1}^{2} y_{t-2}+\beta_{1} \varepsilon_{t-1}+\varepsilon_{t} \\ & =\left(\beta_{0}+\beta_{0} \beta_{1}\right)+\beta_{1}^{2}\left(\beta_{0}+\beta_{1} y_{t-3}+\varepsilon_{t-2}\right)+\beta_{1} \varepsilon_{t-1}+\varepsilon_{t} \\ & =\beta_{0}\left(1+\beta_{1}+\beta_{1}^{2}\right)+\beta_{1}^{3} y_{t-3}+\beta_{1}^{2} \varepsilon_{t-2}+\beta_{1} \varepsilon_{t-1}+\varepsilon_{t} \\ & =\cdots \\ & =\beta_{0}\left(1+\beta_{1}+\beta_{1}^{2}+\cdots\right)+\varepsilon_{t}+\beta_{1} \varepsilon_{t-1}+\beta_{1}^{2} \varepsilon_{t-2}+\beta_{1}^{3} \varepsilon_{t-3}+\cdots \\ & =\frac{\beta_{0}}{1-\beta_{1}}+\varepsilon_{t}+\beta_{1} \varepsilon_{t-1}+\beta_{1}^{2} \varepsilon_{t-2}+\beta_{1}^{3} \varepsilon_{t-3}+\cdots\end{aligned}
平稳的 AR(1)看成是过去所有扰动项的总效应之和，离现在越远的扰动项其影响力呈几何级数递减。
“动态乘子”(dynamic multiplier)：
\operatorname{IRF}(j) \equiv \frac{\partial y_{t+j}}{\partial \varepsilon_{t}}=\beta_{1}^{j} 当第 t 期的扰动项\varepsilon_{t}变化 1 单位时(而其他期的扰动项均不变)，对相隔 j 期的y_{t+j}的影响
动态乘子与绝对时间 t 无关，是时间间隔 j 的函数
将\left(\partial y_{t+j} / \partial \varepsilon_{t}\right)视为 j 的函数，称为“脉冲响应函数”(Impulse Response Function，简记 IRF)
刻画了y_{t+j}对\varepsilon_{t}的 1 单位脉冲(impulse)的响应(response)
脉冲响应图：\left(j, \partial y_{t+j} / \partial \varepsilon_{t}\right)
**将 AR(p)视为一维的向量自回归（Vector Autoregression，简记 VAR）
varbasic x y z,lags(numlist) irf
//“varbasic”为估计 VAR 模型的便捷命令，而“x y z”为 VAR 模型所包含的变量(此例中只有一个变量y)
//选择项“lags(numlist)”表示滞后阶数，默认为“lags(1 2)”或“lags(1/2)”，即滞后二阶
// 选择项“irf”表示画脉冲响应图
varbasic dlny if year<2013,lags(1) irf //估计 dlny 的 AR(1)模型,画 IRF 图
//使用命令varbasic的 估计系数 与命令“reg dlny l.dlny”完全相同，只是命令 varbasic 不提供异方差稳健标准误的选择项(时间序列一般不存在异方差问题)
//某阶的置信区间包含0，则表示其不显著
varbasic dlny if year<2013,irf //估计 dlny 的 AR(2)模型，画 IRF 图
//没有指定滞后阶数，默认使用2阶
//AR(2)模型的脉冲响应函数不再单调递减，更具动态特征，先下降，变为负数后再反弹上升，又下降并趋于0
一维的向量自回归（Vector Autoregression，简记 VAR）
13.8 向量自回归过程(Vector Autoregression, VAR)
“多变量时间序列”(multivariate time series)：同时关心几个变量的预测，将这些变量放在一起，作为一个系统来预测，使得预测相互自洽(mutually consistent)
假设有两个时间序列\left\{y_{1 t}, y_{2 t}\right\}分别作为两个回归方程的被解释变量，解释变量为这两个变量的 p 阶滞后值，构成二元的 VAR(p)系统：：
\left\{\begin{array}{c}y_{1 t}=\beta_{10}+\beta_{11} y_{1, t-1}+\cdots+\beta_{1 p} y_{1, t-p}+\gamma_{11} y_{2, t-1}+\cdots+\gamma_{1 p} y_{2, t-p}+\varepsilon_{1 t} \\ y_{2 t}=\beta_{20}+\beta_{21} y_{1, t-1}+\cdots+\beta_{2 p} y_{1, t-p}+\gamma_{21} y_{2, t-1}+\cdots+\gamma_{2 p} y_{2, t-p}+\varepsilon_{2 t}\end{array}\right.
\left\{\varepsilon_{1 t}\right\} 与 \left\{\varepsilon_{2 t}\right\}均为白噪声(无自相关)
两个方程的扰动项之间可以存在“同期相关性”(contemporaneous correlation)：
\operatorname{Cov}\left(\varepsilon_{1 t}, \varepsilon_{2 s}\right)=\left\{\begin{array}{ll}\sigma_{12} & \text { 若 } t=s \\ 0 & \text { 其他 }\end{array}\right.
\begin{aligned}\left(\begin{array}{l}y_{1 t} \\ y_{2 t}\end{array}\right)= & \left(\begin{array}{l}\beta_{10} \\ \beta_{20}\end{array}\right)+\left(\begin{array}{l}\beta_{11} \\ \beta_{21}\end{array}\right) y_{1, t-1}+\cdots+\left(\begin{array}{c}\beta_{1 p} \\ \beta_{2 p}\end{array}\right) y_{1, t-p} \\ & +\left(\begin{array}{l}\gamma_{11} \\ \gamma_{21}\end{array}\right) y_{2, t-1}+\cdots+\left(\begin{array}{l}\gamma_{1 p} \\ \gamma_{2 p}\end{array}\right) y_{2, t-p}+\left(\begin{array}{c}\varepsilon_{1 t} \\ \varepsilon_{2 t}\end{array}\right)\end{aligned}
\left(\begin{array}{l}y_{1 t} \\ y_{2 t}\end{array}\right)=\left(\begin{array}{l}\beta_{10} \\ \beta_{20}\end{array}\right)+\left(\begin{array}{ll}\beta_{11} & \gamma_{11} \\ \beta_{21} & \gamma_{21}\end{array}\right)\left(\begin{array}{l}y_{1, t-1} \\ y_{2, t-1}\end{array}\right)+\cdots+\left(\begin{array}{ll}\beta_{1 p} & \gamma_{1 p} \\ \beta_{2 p} & \gamma_{2 p}\end{array}\right)\left(\begin{array}{l}y_{1, t-p} \\ y_{2, t-p}\end{array}\right)+\left(\begin{array}{l}\varepsilon_{1 t} \\ \varepsilon_{2 t}\end{array}\right)
定义各向量，得：\boldsymbol{y}_{t}=\Gamma_{0}+\Gamma_{1} \boldsymbol{y}_{t-1}+\cdots+\Gamma_{p} \boldsymbol{y}_{t-p}+\varepsilon_{t}
形式与 AR(p)相似，故名“VAR(p)”
{\varepsilon_{t}}是白噪声的推广，称为“向量白噪声过程”(vector white noise process)，或“新息过程”(innovation process)
解释变量\left\{\boldsymbol{y}_{t-1}, \cdots, \boldsymbol{y}_{t-p}\right\}依赖于\left\{\boldsymbol{\varepsilon}_{t-1}, \boldsymbol{\varepsilon}_{t-2}, \cdots\right\}，而{\varepsilon_{t}}与\left\{\boldsymbol{\varepsilon}_{t-1}, \boldsymbol{\varepsilon}_{t-2}, \cdots\right\}不相关，可视所有解释变量为前定变量，与当期扰动项\varepsilon_{t}不相关，因此可用 OLS 对每个方程分别进行一致估计
在 VAR 建模时，需确定变量的滞后阶数，及包含几个变量
滞后阶数的选择
方法一、使用信息准则
方法二、检验最后一阶系数的显著性(由大到小的序贯规则)
方法三、检验 VAR 模型的残差是否为白噪声（是否有自相关），如果存在自相关，应加入更高阶的滞后变量。
13.9 VAR 的脉冲响应函数
n 元 VAR(p)系统：\boldsymbol{y}_{t}=\Gamma_{0}+\Gamma_{1} \boldsymbol{y}_{t-1}+\cdots+\Gamma_{p} \boldsymbol{y}_{t-p}+\varepsilon_{t}
VAR(p)系统也可写成“向量移动平均过程”(Vector Moving Average Process) VMA(∞)的形式：
y_{t}=\alpha+\varepsilon_{t}+\psi_{1} \varepsilon_{t-1}+\psi_{2} \varepsilon_{t-2}+\cdots=\alpha+\sum_{i=0}^{\infty} \psi_{i} \varepsilon_{t-i}
\boldsymbol{\psi}_{0} \equiv \boldsymbol{I}_{n}, 而 \boldsymbol{\psi}_{j} 为 n 维方阵
边际效应 \frac{\partial \boldsymbol{y}_{t+s}}{\partial \boldsymbol{\varepsilon}_{t}^{\prime}}=\boldsymbol{\psi}_{s}
13.10 格兰杰因果检验
考虑 ADL(p, p)模型：
y_{t}=\gamma+\sum_{m=1}^{p} \alpha_{m} y_{t-m}+\sum_{m=1}^{p} \beta_{m} x_{t-m}+\varepsilon_{t}
原假设H_{0}: \beta_{1}=\cdots=\beta_{p}=0 x 的过去值对预测 y 的未来值有无帮助
如拒绝 H0 ，称 x 是 y 的“格兰杰因”(Granger cause)
格兰杰因果关系并非真正意义上的因果关系，充其量只是动态相关关系，表明一个变量是否对另一变量有“预
测能力”(predictability)。
格兰杰因果检验仅适用于平稳序列，或者有协整关系的单位根过程
对于不存在协整关系的单位根变量，则只能先差分，得到平稳序列后再进行格兰杰因果检验
**假设变量为x y z
varsoc x y z, maxlag(#) //计算不同滞后期的信息准则
//“soc”表示selection-order criteria
//“maxlag(#)”表示最大滞后期，默认值为4
//结果中主要看AIC和SBIC，最小值已由*标注
varbasic x y z,lags(numlist) irf //估计 VAR 模型的便捷命令
//选择项“lags(numlist)”表示滞后阶数，默认为“lags(1 2)”或“lags(1/2)”，即滞后二阶
//“irf”表示画(未正交化)脉冲响应图，默认为“oirf”(画正交化脉冲响应图)
var x y z, lags(numlist) exog(w1 w2) //估计 VAR 的正式命令
//选择项“lags(numlist)”表示滞后阶数，默认为“lags(1/2)”，即滞后二阶
//如果要滞后三阶，可使用选择项“lags(1/3)”
//“exog(w1 w2)”表示在 VAR 模型中引入外生变量w1,w2
varlmar //对残差是否存在自相关进行LM 检验,ar表示Autocorrelation
varstable,graph //估计 VAR 后，通过特征值检验该 VAR 系统是否为平稳过程
//如果所有特征值都在单位圆内部，则为平稳过程
//选择项“graph”表示画出特征值的几何分布图
varwle
//估计 VAR 后，对每个方程以及所有方程的各阶系数的联合显著性进行沃尔德检验
//“wle”表示 Wald lag-exclusion statistics
vargranger //估计 VAR 后，进行格兰杰因果检验
**有关VAR脉冲响应函数的命令
irf create irfname, set(filename) step(#) replace order(varlist) //计算脉冲响应函数，命名为irfname
//选择项“set(filename)”表示建立脉冲文件“filename”，使之成为当前的脉冲文件(make filename active)，并将脉冲结果“irfname”存入此脉冲文件“filename”(若未使用选择项“set(filename)”指定脉冲文件，则将脉冲响应结果存入当前的脉冲文件)
//一个脉冲文件"filename"可存储多个脉冲响应结果
//“step(#)”表示考察截止#期的脉冲响应函数，默认为“step(8)”
//“replace”表示替代已有的同名脉冲响应结果irfname (如果有)
//“order(varlist)”指定变量排序，默认使用估计VAR时的变量排序计算正交化IRF
irf graph irf,impulse(varname) response(varname) //画脉冲响应图(未正交化)
//第一个irf表示总命令，第二个irf表示画脉冲响应图
//“impulse(varname)”用于指定脉冲变量
//“response(varname)”用来指定反应变量,默认画出所有变量的脉冲响应图
irf graph oirf,impulse(varname) //画正交化的脉冲响应图
//“正交化的脉冲响应函数”(Orthogonalized Impulse ResponseFunction，简记 OIRF)
//如将以上两个命令中的“irf graph”改为“irf table”，则将相应信息列表（可查看取值）而非画图
fcast compute prefix,step(#)
//估计 VAR 后，计算被解释变量未来#期的预测值，并把预测值赋予被解释变量加上前缀“prefix”(自行确定)的变量名
fcast graph varlist,observed
//运行命令“fcast compute”后，将变量“varlist”的预测值画图
//选择项“observed”表示与实际观测值相比较
macro_swatson.dta
use macro_swatson.dta,clear
//由于通胀率 inf 可能不平稳，考虑其一阶差分dinf与失业率 unem 构成的二元 VAR 系统
tsline dinf unem,lpattern(solid dash)
varsoc dinf unem //根据信息准则，估计此 VAR 系统的阶数
//当p=2时(上表打星号“*”者)，AIC 与 BIC 信息准则最小化
var dinf unem,lags(1/2) //估计二阶向量自回归模型，也可不加lags(1/2)，因为不加就默认为2阶
varwle //检验各阶系数的联合显著性
//结果显示无论单一方程，还是两个方程作为整体，各阶系数均高度显著
varlmar //检验残差是否为白噪声，即残差是否存在自相关
//P值大于0.1，可接受残差“无自相关”的原假设
varstable,graph //检验 VAR 系统是否为平稳过程，并画图
//结果显示所有特征值均在单位圆之内，故 VAR 模型满足平稳性条件
vargranger //考察变量 dinf 与 unem 之间的格兰杰因果关系
//第一行为：解释变量为dinf，假设unem前的各阶系数均为0（excluded）
//P值约为0，拒绝原假设：unem前的各阶系数均为0，可认为unem为dinf的格兰杰原因
//第二行all表示：解释变量为dinf，假设所有其他解释变量前的各阶系数均为0（本例中除了dinf，只有unem一个其他解释变量，所以结果与第一行相同）
//结果显示二者互为格兰杰原因，格兰杰因果检验无法提供变量排序的信息，因此计算正交化脉冲响应函数时，需要改变变量的排序，进行稳健性检验
irf create iu, set(macro) //计算脉冲响应函数，将脉冲文件命名为“macro”，并将脉冲结果命名为“iu”(表示变量排序为 dinf，unem)
irf graph oirf,yline(0) //画正交化的脉冲响应图
//选择项“yline(0)”表示在纵轴 y=0 处画一条水平线
//Graphs by irfname,impulse variable,and response variable:小图的标题命名顺序为“脉冲名称、冲击变量、响应变量”
//左下小图标题为“iu,unem,dinf”，表明此图为根据脉冲结果 iu，冲击变量 unem，响应变量 dinf 所画的脉冲响应图,失业率 unem 的一个标准差的正向冲击，将使未来一期的 dinf 下降,但未来二期的 dinf 即反弹，然后此影响逐渐消失归于 0
irf create ui, order(unem dinf) //变换变量次序，考察正交化脉冲响应函数的稳健性,此命令将脉冲结果记为ui
//选择项“order(unem dinf)”表示变量 unem 排在 dinf 之前
irf graph oirf,i(unem) r(dinf) yline(0) //比较两种变量排序下的脉冲响应图(unem→dinf)
//在不同变量排序下，变量 dinf 对于 unem 冲击的脉冲响应差别不大，只是在反应幅度上略有不同
irf graph oirf,i(dinf) r(unem) yline(0) //比较两种变量排序下的脉冲响应图(dinf→unem)（另一个方向）
//变量排序对于从 dinf 至 unem 的脉冲响应幅度有较大影响，但变动方向上依然类似
//左图所有脉冲响应函数均不显著（0都在置信区间内），右图第一期的脉冲响应函数显著
**仅用1999年以前的数据估计VAR模型，然后预测1999年1季度-2002年1季度的10 个季度，并与实际观测值比较
varbasic dinf unem if quarter<tq(1999q1),lags(1/2) nograph
//“tq(1999q1)”表示季度数据格式
//选择项“nograph”表示不画脉冲响应图
fcast compute f_,step(10) //预测未来 10 个季度的变量取值，分别记为“f_dinf”与“f_unem”,分别为对 dinf 与 unem 的预测值
fcast graph f_dinf f_unem,observed lpattern(dash) //将 dinf 与 unem 的预测值画图
//选择项“observed”表示显示变量的实际观测值
//选择项“lpattern(dash)”表示以虚线来画变量的预测值(以便区别于实线的实际观测值)
//结果显示对通胀率变动的预测准确度优于对失业率的预测
//预测的时期越长，则预测的精确度越低
13.12 时间趋势项
set obs 100 //假设样本容量为 100
gen t=_n //设定第n个观测值的t为n
gen t2=t^2
corr t t2 //计算时间趋势t及其平方项t^2的相关系数
//存在严重的多重共线性
//常见的做法是仅包含线性趋势项，以避免多重共线性
宏观经济变量通常都有时间趋势，比如 y 与 x 都有时间趋势，故简单地将 y 对 x 进行回归将发现二者存在显著关系，而事实上只是因为共同的时间趋势所驱动。这是一种伪回归”(spurious regression)
将遗漏的时间趋势 t 加入回归方程，可消除此伪回归线现象
y_{t}=\alpha+\beta x_{t}+\gamma t+u_{t}
扰动项 ut 不再包含时间趋势，不会与 xt 相关，故 OLS 一致
**对 GDP 对数(lny)建模，估计 lny 的 AR(2)模型，并加上时间趋势项
use gdp_china.dta,clear
gen t=_n //定义时间趋势项t
reg lny l(1/2).lny t if year<2013,r
//lny 的两阶滞后以及时间趋势项都高度显著
predict lny3
dis exp(lny3[36])
dis y[36]-exp(lny3[36])
13.13 季节调整(seasonal adjustment)
包含季节效应的时间序列不能直接计算环比增长率
季节调整通过估计“季节因子”(seasonal factor)来进行
“加法季节因子”(additive seasonal factor)与“乘法季节因子”(multiplicative seasonal factor)
“加法模型”(additive model)
Y_{t}=T C_{t}+S_{t}+I_{t}
Y 为原序列，TC 为趋势循环要素， S 为季节要素，而 I 为不规则要素
“乘法模型”(multiplicative model)
Y_{t}=T C_{t} \times S_{t} \times I_{t}
“对数加法模型”
\ln Y_{t}=\ln T C_{t}+\ln S_{t}+\ln I_{t}
回归法
首先生成月度(或季度)虚拟变量，然后把时间序列对这些虚拟变
量进行 OLS 回归，所得残差就是经季节调整后的序列。
use airpassengers.dta,clear
tsset time
tsline airpassengers
**生成月度虚拟变量
gen month=month((dofm(time))) //从时间变量 time 提取月度信息，记为变量 month
//第一个month为自定义变量名，第二个month为stata的函数，从时间变量time中提取月份
tab month,gen(m) //使用命令 tab 生成月度虚拟变量
//选择项“gen(m)”表示，根据变量 month 的不同取值，生成相应的虚拟变量，记为 m1, m2 , ... , m12，分别对应于 12 个月
reg airpassengers m2-m12 //以 1 月份为参照值，把变量 airpassengers 对第 2-12 月的月度虚拟变量进行回归
//七、八两月的虚拟变量(m7 与 m8)均在 5%水平上显著为正；其他月份的虚拟变量不显著
predict air_sa,r //计算上述回归的残差项(记为 air_sa)
sum airpassengers
gen airpassengers_sa = air_sa+r(mean) //季节调整序列: OLS 残差项的平均值一定为 0，故需把原序列的均值加回去
tsline airpassengers_sa airpassengers,lpattern(dash) //将回归法的季节调整序列与原序列画图
13.14 日期数据的导入
年度数据，“tsset year”将变量 year 设为时间变量即可(假设时间变量为 year)
“1949-10-01”或“1949/10/01”的时间变量，在导入 Stata 后，可能被视为“字符串”(string)，而非“数
字型”(numeric)，无法直接对其进行运算(比如，滞后一期)
**日度数据(daily data)，可使用如下命令转换为“整数日期变量”(integer date variable)
gen newvar=date(varname,"YMD")
//函数“date”表示转换为日期变量
//“varname”为原来的时间变量
//“newvar”为新定义的时间变量,将以“整数日期”的形式显示
//"YMD"告诉 Stata，原始数据的格式为“年-月-日”。如数据格式为“月-日-年”，则应为"MDY"，以此类推
//在 Stata 内部，所有日期变量的存储格式均为“elapsed dates”(Stata 称为 Stata Internal Format，简记 SIF)，即从 1960 年 1 月 1日以来过了多少天
format newvar %td //以通常的日期格式(Human Readable Format，简记 HRF)显示
//“%td”中的“d”即表示“date”
**月度数据(monthly data)
gen newvar=monthly(varname,"YM")
//函数“monthly”表示转换为月度变量
//"YM"告诉 Stata，原始数据的格式为“年-月”
format newvar %tm //以日期格式在 Stata 中显示
**季度数据(quarterly data)
gen newvar=quarterly(varname,"YQ")
//函数“quarterly”表示转换为季度变量
//"YQ"告诉Stata，原始数据的格式为“年-季”
format newvar %tq //以日期格式在 Stata 中显示
//“%tq”中的“q”即表示“quarter”
**在原始数据中，年、月、日分别以数字型(numeric)变量“Y, M, D”来表示，可使用以下命令将其合成单一的日期变量
gen newvar=mdy(M,D,Y)
help date //有关日期数据的更多更多说明
14.单位根与协整
“非平稳序列”(non-stationary time series)
一、确定性趋势(deterministic trend)
y_{t}=\beta_{0}+\beta_{1} t+\varepsilon_{t}
二、结构变动(structural break)
y_{t}=\left\{\begin{array}{ll}\alpha_{1}+\beta_{1} x_{t}+\varepsilon_{t}, & \text { 若 } t<\bar{t} \\ \alpha_{2}+\beta_{2} x_{t}+\varepsilon_{t}, & \text { 若 } t \geq \bar{t}\end{array}\right.
三、随机趋势(stochastic trend)
y_{t}=y_{t-1}+\varepsilon_{t}
\begin{array}{l}y_{1}=y_{0}+\varepsilon_{1} \\ y_{2}=y_{1}+\varepsilon_{2}=y_{0}+\varepsilon_{1}+\varepsilon_{2} \\ y_{3}=y_{2}+\varepsilon_{3}=y_{0}+\varepsilon_{1}+\varepsilon_{2}+\varepsilon_{3} \\ \quad \vdots \\ y_{t}=y_{t-1}+\varepsilon_{t}=y_{0}+\varepsilon_{1}+\cdots+\varepsilon_{t}=y_{0}+\sum_{s=1}^{t} \varepsilon_{s}\end{array}
如果\varepsilon_{1}增加一单位,所有\left\{y_{1}, y_{2}, \cdots, y_{t}, \cdots\right\}都将增加一个单位
\left\{\varepsilon_{t}\right\}为此模型的"随机趋势"：来自\left\{\varepsilon_{t}\right\}的任何扰动对\left\{y_{t}\right\}都有永久效应(permanent effect),影响力不随时间而衰减
\operatorname{Var}\left(y_{t}\right)=\operatorname{Var}\left(\sum_{s=1}^{t} \varepsilon_{s}\right)=\sum_{s=1}^{t} \operatorname{Var}\left(\varepsilon_{s}\right)=t \sigma_{\varepsilon}^{2}
当 t \rightarrow \infty 时, \operatorname{Var}\left(y_{t}\right) \rightarrow \infty (方差发散), 故 \left\{y_{t}\right\} 非平稳
“带漂移的随机游走”(random walk with drift)：包含常数项
y_{t}=\beta_{0}+y_{t-1}+\varepsilon_{t}
\mathrm{E}\left(y_{t}\right)=\beta_{0}+y_{t-1}
\Delta y_{t}=\beta_{0}+\varepsilon_{t}
随机游走的差分为平稳序列，随机游走也称为“差分平稳”(difference stationary)序列
AR(1)模型：y_{t}=\beta_{0}+\beta_{1} y_{t-1}+\varepsilon_{t}，\beta_{1}=1则为随机游走
平稳的时间序列为“零阶单整”(Integrated of order zero)，记为 I(0)
如果时间序列的一阶差分为平稳过程，称为“一阶单整”(Integrated of order one)，记为 I(1)，也称为“单位根过程”(unit root process)
如果时间序列的 d 阶差分为平稳过程，称为“d 阶单整”(Integrated of order d)，记为 I(d)
对于 I(0)序列，由于它是平稳的，故长期而言有回到其期望值的趋势。这种性质称为“均值回复”(mean-reverting)。
非平稳的 I(1)序列会“到处乱跑”(wander widely)，没有上述性质。比如，随机游走的方差越来越大，趋向无穷。
I(0)序列对过去行为只有有限记忆，即发生在过去的扰动项对未来的影响随时间而衰减。
I(1)序列则对过去行为有无限长的记忆，即任何过去的冲击都将永久地改变未来的整个序列。
如果时间序列 \left\{y_{t}\right\} 的 d 阶差分为平稳的 \operatorname{ARMA}(p, q) 过程，则称 \left\{y_{t}\right\} 为 \operatorname{ARIMA}(p, d, q) 过程。(1) 自回归系数的估计量不服从渐近正态分布，t 检验失效
$
y_{t}=y_{t-1}+\varepsilon_{t} $
y_{1}=\varepsilon_{1}, \quad y_{2}=\varepsilon_{1}+\varepsilon_{2},\quad \cdots \cdots,\quad y_{t}=\varepsilon_{1}+\cdots+\varepsilon_{t}=\sum_{s=1}^{t} \varepsilon_{s}
program randwalk,rclass //(定义程序“randwalk”，并以 r()形式储存结果)
drop _all //(删去内存中已有数据)
set obs 1000 //(确定样本容量为 1000)
gen eps=rnormal() //(产生服从标准正态分布的扰动项)
gen y=sum(eps) //(假设y0=0,定义随机游走)
gen t=_n //(定义时间变量，第t期即第i个观测值)
tsset t //(将数据设为时间序列，以便使用滞后算子)
reg y L.y //(回归)
return scalar b1=_b[L.y] //(记 OLS 系数为 b1)
end //(程序结束)
simulate beta=r(b1),seed(10101) reps(1000):randwalk
kdensity beta //把系数的核密度图画出来，结果显示左偏，并非正态分布
(2) 两个相互独立的单位根变量可能出现伪相关或伪回归
两个单位根过程\quad y_{t}=y_{t-1}+u_{t} ; \quad x_{t}=x_{t-1}+v_{t}
假设 y_{0}=0, x_{0}=0, 则 y_{t}=\sum_{s=1}^{t} u_{s}, x_{t}=\sum_{s=1}^{t} v_{s}
drop _all //(删去内存中已有数据)
set obs 10000 //(确定样本容量为 10,000)
set seed 1234 //(确定随机数的种子为 1234)
gen u=rnormal() //(产生服从标准正态分布的扰动项ut)
gen y=sum(u) //(定义随机游走yt)
set seed 12345 //(确定随机数的种子为12345,保证xt与yt相互独立)
gen v=rnormal() //(产生服从标准正态分布的扰动项vt)
gen x=sum(v) //(定义随机游走xt)
reg y x
//回归系数在1%水平上显著,存在"伪回归"
//存在随机趋势，所以高度相关，如果为平稳过程，不会如此
reg u v
//扰动项的回归系数高度不显著
**画出yt与xt的时间趋势图
gen t=_n //(定义时间变量t)
line y x t,lpattern(dash)
//尽管yt与xt相互独立，但二者都是单位根过程，存在相似的时间趋势，从而出现伪相关与伪回归，从而误导统计推断
如何避免伪相关或伪回归？
方法之一：先对 I(1)变量作差分，得到平稳 I(0)序列，再作回归。
方法之二：“协整”(cointegration)。须先检验是否存在单位根。
14.5 单位根检验
Augmented Dickey-Fuller 单位根检验
\Delta y_{t}=\beta_{0}+\delta y_{t-1}+\gamma_{1} \Delta y_{t-1}+\gamma_{2} \Delta y_{t-2}+\cdots+\gamma_{p-1} \Delta y_{t-p+1}+\gamma t+\varepsilon_{t}
H_{0}: \delta=0 \quad vs \quad H_{1}: \delta<0
dfuller y,lags(p) regress noconstant drift trend
//选择项“lags(p)”表示包含 p阶滞后差分项，默认为“lags(0)”，对应于 DF 检验
//选择项“regress”表示显示回归结果
//选择项“noconstant drift trend”(三者最多选一项，不能并用),含义见上表
Schwert(I989)建议的最大滞后阶数p_{\max }=\left[12 \cdot(T / 100)^{1 / 4}\right]
use nelson_plosser.dta,clear
tsline lrgnp lun if year>=1890,lp(dash) xlabel(1890(10)1970)
//变量 lrgnp 的取值始于 1890 年，变量 lun 的取值始于 1909 年,“if year>=1890”来限制画图的范围(使之更加美观)
//图中显示实际 GNP 对数(虚线)有明显的上升趋势，可能有时间趋势项
dfuller lrgnp,trend //对于实际 GNP 对数，首先考虑带常数项与时间趋势项的DF检验
//DF 统计量 Z(t)为–2.026 > –3.489(左边单侧检验)，故可在5%的水平上接受“存在单位根”的原假设
**扰动项可能存在自相关，考虑更高阶的 ADF 检验
di 12*(62/100)^(1/4)
//10.65,p=10,lag中的参数选择9
dfuller lrgnp,lags(9) trend reg
//reg表示汇报ADF检验回归的结果
//时间趋势项(_trend)很显著( p值为 0.017)，但最后一阶滞后项(L9D.)在5%的水平上并不显著
//依次令p=9,8,7,...,3，进行 ADF 检验，最后一阶滞后项仍不显著
dfuller lrgnp,lags(1) trend reg //令p=2再进行 ADF 检验
//最后一阶滞后项(LD.)在 1%的水平上显著地不等于0
//ADF 统计量 Z(t)为–2.994 > –3.490，无法在 5%的水平上拒绝单位根的原假设
**进一步检验 lrgnp 的一阶差分是否为平稳过程
dfuller d.lrgnp
//ADF 统计量 Z(t)为–5.322 < –3.566，可在 1%的水平上拒绝单位根的原假设，认为lrgnp 为平稳过程
**考察失业率对数 lun 是否含有单位根
dfuller lun,lags(3) reg //根据序贯t规则，选择p=4,进行不带时间趋势的 ADF 检验
//ADF 统计量 Z(t)为–3.588 < –3.542，可在 1%的水平上拒绝单位根的原假设，认为失业率对数 lun 为平稳过程
14.7 协整的思想与初步检验
如果多个单位根序列拥有“共同的随机趋势”(common stochastic trend)，则可对这些变量作适当的线性组合而消去此随机趋势，从而得到平稳序列
use macro_3e.dta,clear
tsset time
tsline fygm3 fygt1,lpattern(dash)
//考察 fygm3 (3 个月国库券利率)与 fygt1 (1 年期国库券利率)的时间趋势
两个 I(1)过程\left\{y_{t}\right\},\left\{x_{t}\right\}：
\left\{\begin{array}{l}y_{t}=\alpha+\beta w_{t}+\varepsilon_{t} \\ x_{t}=\gamma+\delta w_{t}+u_{t}\end{array}\right.w_{t} 为随机游走, w_{t}=w_{t-1}+v_{t} ; \varepsilon_{t}, u_{t}, v_{t} 均为白噪声
\delta y_{t}-\beta x_{t}=(\alpha \delta-\beta \gamma)+\left(\delta \varepsilon_{t}-\beta u_{t}\right)
EG-ADF 检验：判断变量间是否存在协整关系
Stata 的默认方法：“迹检验”(trace test)是似然比检验，为单边右侧检验，即 trace 越大，越倾向于拒绝原假设
协整分析的起点是单位根检验。如确定所有变量都为单位根过程，可进一步检验这些 I(1)变量是否存在协整关系
VAR 模型：\boldsymbol{y}_{t}=\alpha+\Phi_{1} \boldsymbol{y}_{t-1}+\Phi_{2} \boldsymbol{y}_{t-2}+\cdots+\Phi_{p} \boldsymbol{y}_{t-p}+\boldsymbol{\varepsilon}_{t}
VAR 模型对应的向量误差修正模型(Vector Error Correction Model，VECM):
\Delta y_{t}=\alpha+\Gamma_{0} y_{t-1}+\Gamma_{1} \Delta y_{t-1}+\cdots+\Gamma_{p-1} \Delta y_{t-p+1}+\varepsilon_{t}
假设变量为x,y,z
(1) 检验协整秩
vecrank x y z,lags(#) max trend(none) trend(trend)
//命令vecrank的输出结果将列出"h=0,1,…,n-1"的一系列检验,并以星号(*)标出所接受h值,即协整秩.
//选择项“lags(#)”表示对应的 VAR 模型中滞后的阶数，默认为“lags(2)”
//选择项“max”表示也进行最大特征值检验，默认仅进行迹检验
//选择项“trend(none)”表示不包括常数项或时间趋势
//选择项“trend(trend)”表示包括常数项与时间趋势；默认包括常数项，但不包括时间趋势
(2) 估计协整关系
** MLE 估计 VECM 模型
vec x y z,lags(#) rank(#) trend(none) trend(trend) sindicators(varlist)
//选择项“lags(#)”表示对应的 VAR 模型中滞后的阶数，默认为“lags(2)”
//选择项“rank(#)”表示协整秩的阶数，默认为“rank(1)”
//选择项“trend(none)”表示不包括常数项或时间趋势
//选择项“trend(trend)”表示包括常数项与时间趋势；默认包括常数项，但不包括时间趋势
//选择项“sindicators(varlist)”表示加入季节虚拟变量
(3) 诊断性检验
veclmar //对残差是否存在自相关进行LM 检验
vecstable,graph //检验 VECM 系统是否为平稳过程
//选择项“graph”表示画出特征值的几何分布图
(4) 脉冲响应函数与预测
//估计 VECM 模型后，可计算脉冲响应函数或进行预测，命令与VAR 模型相同(见第 13 章)
//一般主要关注长期均衡关系(协整关系)，不太关心短期调整过程。
use mpyr.dta,clear
tsline logmr logy r,lpattern(solid dash shortdash) xlabel(1900(10)1990)
//从图形上大致考察(logmr, logy, r)是否存在协整关系
//选择项“lpattern(solid dash shortdash)”表示分别用实线、虚线与短虚线来画图
//“label(1900(10)1990)”表示，在横轴从 1900-1990 年，每隔 10 年做个标注(label)
**确定协整秩，即有多少个线性无关的协整关系
varsoc logmr logy r //计算信息准则，检验该系统所对应的 VAR 模型的滞后阶数
//打星号者为根据不同准则所选的滞后阶数
**进行协整秩检验
vecrank logmr logy r,lags(2) trend(trend) max
//选择项“lags(2)”表示对应的 VAR 模型滞后二阶(也是默认值，可省略)；
//选择项“trend(trend)”表示既包括常数项，也包括时间趋势项；
//选择项“max”表示显示最大特征值统计量
//结果中上半部分为迹检验，下半为特征值检验
**MLE 方法估计向量误差修正模型(VECM)
vec logmr logy r,lags(2) rank(1)
//选择项“lags(2)”表示对应的 VAR 模型滞后二阶(也是默认值，可省略)
//选择项“rank(1)”表示协整秩为 1(也是默认值，可省略)
//下半部分Cointegrating equations为协整关系，对协整方程的估计结果
//根据此协整方程，协整向量为(1, -0.98, 0.11)
//logmr 的系数被标准化为 1，故标准误缺失
//其他两个变量(logy 与 r)的协整系数均在 1%水平上显著
veclmar //检验 VECM 模型的残差是否存在自相关
//可接受“无自相关”的原假设
vecstable,graph //检验 VECM 系统是否稳定
//除了 VECM 模型本身所假设的单位根之外，伴随矩阵的所有特征值均落在单位圆之内，故是稳定系统
reg logmr logy r //直接用 OLS 估计长期均衡关系(EG-ADF 两步法)
//OLS 估计值与 MLE 估计结果较接近
//理论上，MLE 估计更有效率
估计的货币需求函数为：
\widehat{\log m r_{t}}=-0.73+0.98 \log y_{t}-0.11 r_{t}
货币需求的收入弹性为 0.98，货币需求的利率半弹性为-0.11，符合经济理论的预期
15.如何做实证研究
实证研究的步骤：准备阶段、选题、探索性研究、收集数据、建立计量模型、选择计量方法、解释回归结果、论文写作、与同行交流、提交论文或投稿
准备阶段
首先，必须掌握一定的经济理论，以获得观察经济现象的必要视角(perspective)、参照系(reference 或 benchmark)与分析工具(analytical tools)。
其次，为进行实证研究，须掌握一定的计量方法与统计软件(比如 Stata)。
选题—具体的“研究问题”(research question)
(1)具体：在“X 对 Y 有何作用”的句型中，应能明确 X 对 Y 具体是什么。
(2)有趣：你的研究问题为什么重要？别人会感兴趣吗？为什么我们要在乎你的问题(Why should we care)？知道问题的答案后，能影响人们对世界某方面的看法吗？
(3)新颖：论文的核心价值在于其创新性，即做出了文献中所没有的边际贡献。这种边际贡献可以是研究了新的现象、使用了新的(更好的)计量方法，或者使用了新的数据集。
(4)可行：即使研究问题很具体、很有趣、很新颖，如找不到相应的数据，则不可行
经济学中文期刊“四大金刚”：《经济研究》、《经济学(季刊)》、《世界经济》、《管理世界》；
经济学英文期刊的“top 5”：American Economic Review，Econometrica，Journal of Political Economy，Quarterly Journal of Economics，Review of Economic Studies。
探索性研究—新颖性与可行性
(1) 通过文献回顾评估选题的新颖性
(2) 确定所需数据是否可得
收集与整理数据
可关注文献中同类研究的数据来源，然后溯本及源
建立计量模型
回归分析一般只能说明变量之间的相关性，要对变量之间的因果关系做出判断，常需要经济理论。
理想情形：从理论模型推导出计量模型(econometric model)，即待估计的回归方程。
设定模型时，应尽量使用常识(common sense)与经济理论(economic theory)
模型既不能过于简单(解释变量过少)，也不宜过于复杂，而应当保持适当的简洁(keep it sensibly simple)。
选择计量方法
对于一般数据，通常先作 OLS，看结果，作为参照系。
被解释变量为虚拟变量，可使用 Probit 或 Logit。
如果是面板数据，应考虑固定效应、随机效应、时间效应等。
如果是时间序列，须先判断是否含单位根，再决定使用相应的计量方法。
从回归分析的相关关系升华到因果关系，是很大的飞跃，需要使用适当的计量方法来识别这种因果关系。
解释回归结果
计量结果可能很复杂，但真正重要的信息通常不多，比如回归系数(含符号)、 p 值，以及样本容量、拟合优度等统计量。
“统计显著性”(statistical significance)主要通过 p 值来考察，如果 p 值小于或等于 0.05，则该系数在统计上显著不等于零；反之，则在统计上不显著，可将此系数视为零(不存在)。
“经济显著性”(economic significance)主要通过系数的绝对值来考察，须特别注意变量的取值单位。
统计上显著而经济上不显著，意味着解释变量对被解释变量的影响很小(经济上不显著)，尽管这种影响的幅度被估计得很精确(统计上显著)。
假设检验，最需关注的只是原假设以及 p 值
诊断性检验
任何计量方法都有适用前提；如前提不成立，则无法使用此计量方法(可能导致不一致估计)。
稳健性检验
“稳健性检验”(robustness check)或“敏感度分析”(sensitivityanalysis)
通过改变样本区间(或去掉极端值)、函数形式、计量方法、控制变量、变量定义、数据来源等，来考察计量结果的稳定性。
只有稳健的结果才有说服力，稳健性检验已成为高质量实证论文不可或缺的一部分。
论文写作
引言(Introduction)
传统套路：提出研究问题后，先回顾已有文献的相关研究及不足之处，顺势引出本文的研究方法与主要贡献(比如，填补了文献的空白)。
传统套路的优点是，较有逻辑性，自然呈现学术发展的脉络；其缺点在于读者需有一定耐心，先回顾文献，才知道本文的主要工作。
现代套路：提出问题后，直奔主题，介绍本文的研究方法与主要结论，再回头介绍本研究与现有文献的关系。
这两种套路各有优缺点，适合不同的论文，但现代套路似乎日益流行。
**文献回顾(Literature Review)
**文献回顾的根本目的是为了厘清本文的研究与已有文献的关系，以凸显本文的边际贡献及在文献中的地位。
**计量模型与估计方法(Econometric model and estimation)
**通常有一个基准(baseline 或 benchmark)的计量模型，然后在此基础上对模型设定有所变化，比如增减或替换变量。
此部分需着重说明论文的估计策略(estimation strategy)，即究竟应使用什么计量方法来识别主要变量之间的因果关系。
**写作风格
**写作过程也是使思路更加清晰的过程。
仔细观察经典论文的文章结构与风格，并注意模仿。
写作伦理
“引用”(citation)与“抄袭”(plagiarism)有何区别？
二者最本质的区别在于，引用给出了信息的出处；而抄袭未提供出处，让读者误以为是作者的原创。
抄袭可定义为“将已经存在的思想或产品‘偷来’作为自己的思想或产品”。
抄袭是一种严重违背学术规范与职业道德的行为，可能导致无可挽回的严重后果。
即使是间接引用，比如用简洁的语言概述前人的思想，或将已有模型作了小的改动，也应及时注明其出处。
欢迎来到这里！
我们正在构建一个小众社区，大家在这里相互信任，以平等 • 自由 • 奔放的价值观进行分享交流。最终，希望大家能够找到与自己志同道合的伙伴，共同成长。
注册 关于