*chp16.2 面板校正标准误
* xtpcse y x1 x2 x3, hetonly
* 其中，hetonly 表示存在组间异方差，但不存在组间同期相关；more既存在组间异方差，也存在组间同期相关
* e.g.
* lnp代表价格效应
* lnpmin 代表替代效应（可以购买更低价的香烟）
* lny 代表收入效应
clear all
cd "…data"
use mus08cigar.dta, clear
* 为考虑时间效应，生成时间趋势变量t
* 然后用LSDV法估计双向固定效应模型（作为对比，先不考虑组间异方差或同期相关）
gen t = year - 62
reg lnc lnp lnpmin lny i.state t, vce(cluster state)
est store OLS
// 回归结果显示，有些州虚拟变量很显著，即存在固定效应，应允许每个州拥有自己的截距项
// 时间效应则不太显著，p值高达0.19
// 然而，上表中的稳健标准误未考虑可能存在的组间异方差与组间同期相关
// 为此，使用面板校正标准误
xtpcse lnc lnp lnpmin lny i.state
est store PCSE
* 从上表可知，命令 xtpcse 的估计系数与命令 reg 的估计系数完全一致，只是标准误不同。
* 在本例中，面板校正标准误似乎反而更小。
*- 16.3 仅解决组内自相关的FGLS
* xtpcse y x1 x2 x3, corr(ar1) corr(psar1)
* 选择项 corr(ar1) 对应于 ρi = ρ 的情形；
* 而选择项 corr(psar1) 则允许每个面板有自己的 ρi, psar1 指 panel-specific arl.
* 如果 T 并不比 n 大很多，则建议使用选择项 corr(ar1), 即约束每个面板（个体）的自回归系数 ρ 均相等。
* 因为时间维度 T 可能无法提供足够的信息来分别估计每个面板自己的 ρi.
* 如果既不指定 corr(ar1)，也不指定 corr(psar1),则进行 OLS 回归。
* 无论使用 OLS 还是 Prais-Winsten 估计法
* 命令 xtpcse 均提供面板校正标准误（在组间异方差与同期相关的情况下仍成立）
* 继续上面的实例，考虑存在组内自相关的情形，但要求各组的自回归指数相同：
xtpcse lnc lnp lnpmin lny i.state t, corr(ar1)
est store AR1
* 从上表可知，是否考虑组内自相关，对估计系数有较大影响。
* 下面考虑允许各组自回归系数的不同的组内自相关情形：
xtpcse lnc lnp lnpmin lny i.state t, corr(psar1)
est store PSAR1
* 列出以上各方法的系数估计值和标准误
local m "OLS PCSE AR1 PSAR1"
esttab `m' , nogap compress replace ///
b(%6.3f) s(N r2_a, fmt(%5.0f %6.3f)) ///
se(%6.3f) ///
keep(`keep') drop(`drop') ///
star(* 0.1 ** 0.05 *** 0.01) ///
mtitle("OLS" "PCSE" "AR1" "PSAR1") ///
addnotes("*** 1% ** 5% * 10%")
* 从上表可知，AR1 与 PSAR1 的系数估计值比较接近，但与 OLS 的结果差别很大。
* 选择 AR1 还是 OLS，需要对组内自相关进行检验
*-16.4 全面FGLS
* 虽然命令 xtpcse 提供了组间异方差与同期相关稳健的面板校正标准误
* 但在进行 FGLS 估计时仅针对组内自相关，并未考虑组间异方差或同期相关
* 更为全面的 FGLS 估计则同时考虑了这三个因素
* 为此，可先进行 OLS 估计，然后使用残差 eit 来估计 εit 的协方差矩阵，依次进行FGLS估计。还可进行迭代 FGLS 估计。
* xtgls y x1 x2 x3, panels(option) corr(option) igls
* panels(iid) 假定不同个体的扰动项为独立同分布的
* panels(het) 假定不同个体的扰动项相互独立但可以有不同的方差
* panels(cor) 假定不同个体的扰动项同期相关且有不同的方差
* corr(ar1) 对应于 ρi = ρ 的组内自相关情形
* corr(psar1) 则允许每个面板有自己的回归系数
* igls 表示使用迭代式 FGLS
* e.g.
* 考虑同时存在组间异方差、同期相关以及组内自相关情形
xtgls lnc lnp lnpmin lny i.state t, panels(cor) cor(ar1)
* 其次，考虑各组自回归系数不同的情形
xtgls lnc lnp lnpmin lny i.state t, panels(cor) cor(psar1)
* 由以上两表可知，各种 FGLS 估计结果不尽相同。
* 究竟使用何种 FGLS 估计，则取决于对组间异方差、组内自相关与组间同期相关的检验。
* 总的来说，"OLS + 面板校正标准误" 最为稳健，全面的 FGLS 估计最有效率，而仅解决组内自相关的 FGLS 则介于二者之间。
*-16.5 组间异方差的检验 (xttest3)
* Greene 提供了对组间异方差的wald 检验
* 原假设 不同个体的扰动项方差均相等
* 此 wald 检验可通过 xttest3 实现
* 该命令只能在运行命令 xtreg, fe 或 xtgls 之后才能使用
* ssc install xttest3
qui xtreg lnc lnp lnpmin lny i.state t, r fe
xttest3
qui xtgls lnc lnp lnpmin lny i.state t
xttest3
* 无论在xtreg, fe 或 xtgls 之后使用 xttest3，均强烈拒绝同方差假设，认为存在组间异方差
*-chp16.6 组内自相关的检验 (xtserial)
* Wooldridge 提供了对组内自相关的 wald 检验
* 原假设 不存在组内自相关
* net install st0039
* xtserial y x1 x2 x3, output
* 由于该命令不接受以 i.state 表示州虚拟变量，故首先生成州虚拟变量，以state1-state10 来表示
tab state, gen(state)
* xtserial lnc lnp lnpmin lny state2-state10 t
* 结果拒绝 “不存在一阶组内自相关” 的原假设
* 由于命令 xtserial 不允许使用 i.state
* 故只能使用 tab 生成的州虚拟变量 state2-state10
*- 16.7 组间同期相关的检验
* 考虑原假设“不存在组间同期相关”。
* 如果原假设成立，则根据残差计算的不同个体扰动项的相关系数接近于0
* Greene 提供了一个对组间同期相关的BP LM检验，可由非官方命令 xttest2 实现。
* xtreg, fe
* xtgls
* ivreg2
clear
use mus08cigar.dta, clear
gen t = year - 62
qui xtreg lnc lnp lnpmin lny t, fe
xttest2
* LM 检验的缺陷在于，仅适合长面板。
* 而检验组间同期相关的 xtscd 也适用于 n 大 T 小的短面板，其中 “csd” 代表 “cross-sectional dependence”。
* 其原假设也是 “不存在组间同期相关”。
* 它包括三种检验方法，分别由 Friedman, Frees, 以及 Pesaran提出，前两种为半参数检验，而后一种为通常的参数检验。
* 命令 xtcsd 须在运行命令 xtreg 之后才能使用。
* ssc install xtcsd, replace
* xtcsd, pesaran abs show
* xtcsd, friedman abs show
* xtcsd, frees abs show
* 其中，选择项 “show” 表示显示残差的相关系数矩阵，而选择项 “abs” 表示显示该矩阵非对角线元素的绝对值之平均。
* 回到数据集 mus08cigar.dta 的例子。
qui xtreg lnc lnp lnpmin lny t, fe
xtcsd, pes
xtcsd, fri
xtcsd, fre abs show
* 以上三种检验的 p 值均小于0.01，故强烈拒绝“无组间同期相关”的原假设，认为存在组间同期相关。
* 上表最后一行显示，非主对角线元素绝对值的平均值高达 0.466.
* 命令“xtcsd, fre abs show” 所显示的残差相关系数矩阵与命令 “xttest2” 所显示完全相同。
* 另外，还可以把面板数据排成“宽型”，然后使用似不相关回归来检验组间截面相关。
*-16.8 变系数模型
* “变系数模型”分为两大类，取决于将“可变系数”视为常数还是随机变量。
* 1. 视为常数
reg lnc lnp lnpmin lny i.state i.state#c.lny t, vce(cluster state)
* i.state 表示州虚拟变量
* i.state#c.lny 表示州虚拟变量与lny的互动项
* c 表示 continuous，用以指明 lny 为连续型变量
* 这些变量在stata中被称为 “因子变量”，是变量窗口所列变量的扩展
* 2. 随机系数模型
* xtrc y x1 x2 x3, beta
* beta 表示，显示对每组系数的估计。
* 如果拒绝原假设，则认为应使用变系数模型。
xtrc lnc lnp lnpmin lny, beta
*-16.9 面板工具变量法
* 1. 对固定效应模型先进行离差变换，再使用 IV
* xtivreg y [varlist1](varlist2 = varlist_iv), fe
* 2. 对固定效应模型先进行一阶差分，再使用 IV
* xtivreg y [varlist1](varlist2 = varlist_iv), fd
* 3. 对随机效应模型使用 IV
* xtivreg y [varlist1](varlist2 = varlist_iv), re
* stata 命令 xtivreg 并不提供 vce(robust) 来保证使用聚类稳健标准误。
* 但可使用 vce(robust) 以得到稳健的 自助标准误。
* 当工具变量个数多于内生解释变量个数时，对面板数据进行 GMM 估计更有效率。
* 可用 xtivreg2，该命令只能处理固定效应模型，先对模型进行 FE 或 FD变换后，再对变换后的模型使用 GMM。
ssc install xtivreg2
xtivreg2 y [varlist1](varlist2 = varlist_iv), fe gmm
xtivreg2 y [varlist1](varlist2 = varlist_iv), fd gmm
* 4. 面板IV的过度识别检验
* 对于面板工具变量的过度识别检验，可涌 xtoverid 实现。
* 该命令须在运行 xtreg, xtivreg, xtivreg2 或 xthtayalr 之后才能使用。
ssc install xtoverid, replace
* 16.10 豪斯曼-泰勒估计量（选读）
* xthtayalr depvar indepvars, endog(varlist) amacurdy
* 16.11 动态面板
* 1. 差分GMM
* 2. 水平GMM
* 3. 系统GMM
* 16.12 动态面板的 stata命令及实例
* xtabond depvar indepvars, lags(p) maxldep(q) twostep vce(robust) pre(varlist) endogenous(varlist) inst(varlist)
* depvar 被解释变量
* indepvars 为满足严格外生性的变量
* lags(p) 使用被解释变量的 p 阶滞后值作为解释变量
* maxldep(q) 表示最多使用q 阶被解释变量的滞后值作为工作变量，默认值为使用所有可能的滞后值
* twostep 表示使用GMM，默认方法为2sls
* pre(varlist) 前定变量
* endogenous(varlist) 内生变量
* inst(varlist) 额外的工具变量
* vce(robust) 稳健标准误
* 系统GMM的stata命令为
* xtdpdsys depvar [indepvars], lags(p) maxldep(q) twostep vce(robust) pre(varlist) endogenous(varlist) inst(varlist)
clear
cd …data
use mus08psidextract.dta, clear
xtabond lwage occ south smsa ind, lags(2) maxldep(3) pre(wks,lag(1,2)) endogenous(ms, lag(0,2)) endogenous(union, lag(0,2)) twostep vce(robust)
est store DiffGMM
* lags(2) 解释变量包含被解释变量的一阶与二阶滞后
* maxldep(3) 表示最多使用被解释变量的三个滞后值作为工具变量
* pre(wks, lag(1,2)) 表示变量wks 及其一阶滞后 L.wks 为前定解释变量，而使用其两个更高阶滞后值为工作变了
* endogenous(ms, lag(0,2)) 指定变量 ms 为内生解释变量
* 通过检验扰动项的差分是否存在一阶与二阶自相关，来检验原假设
estat abond
* 也可以检验扰动项的差分是否存在更高阶的自相关：
estat abond, artests(3)
* artests(3) 表示扰动项差分是否存在一阶至三阶自相关
* 其次，由于差分GMM使用了40个 iv，故需要进行过度识别检验
* 故此，运行 xtabond 命令，但略去选择项 vce(robust)
qui xtabond lwage occ south smsa ind, lags(2) maxldep(3) pre(wks, lag(1,2)) endogenous(ms, lag(0,2)) endogenous(union,lag(0,2)) twostep
estat sargan
// 上表显示，在0.05显著水平上，无法拒绝“所有工具变量均有效”的原假设，因为 p值=0.086 > 0.05
* 下面使用系统GMM来估计以上模型
xtdpdsys lwage occ south smsa ind, lags(2) maxldep(3) pre(wks, lag(1,2)) endogenous(ms, lag(0,2)) endogenous(union, lag(0,2)) twostep vce(robust)
estimates store SysGMM
* 对比差分GMM与系统GMM的系数估计值与标准误
estimates table DiffGMM SysGMM, b se
* 从上表可知，系统GMM的系数估计值与差分 GMM接近，但前者标准误更小
estat abond //对系统GMM扰动项的自相关性进行检验
* 上表显示，可在0.05水平上拒绝，扰动项差分的二阶自相关系数为0的假设。
* 系统GMM的适用性值得怀疑
* 继续进行过度识别检验
qui xtdpdsys lwage occ south smsa ind, lags(2) maxldep(3) pre(wks, lag(1,2)) endogenous(ms, lag(0,2)) endogenous(union, lag(0,2)) twostep
estat sargan
* 系统 GMM 虽然可能更有效率，但其成立的前提似乎不满足
* 为解决扰动项自相关问，在解释变量中加入被解释变量的三阶之后
* 对于内生性解释变量 ms, union仅使用一个更高阶滞后值为工具变量
xtdpdsys lwage occ south smsa ind, lags(3) maxldep(5) pre(wks, lag(1,2)) endogenous(ms, lag(0,1)) endogenous(union, lag(0,1)) twostep vce(robust)
estat abond
* 上表显示，可以接受扰动项无自相关的原假设
qui xtdpdsys lwage occ south smsa ind, lags(3) maxldep(5) pre(wks, lag(1,2)) endogenous(ms, lag(0,1)) endogenous(union, lag(0,1)) twostep
estat sargan
// 结果显示，可以在0.05水平上接受“所有工具变量都有效”的原假设。
* 对修正的模型设定，可以进行系统GMM估计。
* xtabond 与 xtdpdsys 不提供异方差稳健的 Hansen 统计量。
* 故使用非官方命令 xtabond2.
xtabond2 lwage L(1/2).lwage L(0/1).wks ms union occ south smsa ind, gmm(lwage, lag(2 4)) gmm(wks ms union, lag(2 3)) iv(occ south smsa ind) nolevel twostep robust
* gmm(lwage, lag(2 4) 表示使用变量lwage 的2-4阶滞后作为GMM式工具变量
* gmm(wks ms union, lag(2 3)) 使用变量wks, ms 与 union的2-3阶滞后作为GMM式工具变量
* iv(occ south smsa ind) 表示使用变量 occ, south, smsa 与ind 作为自身的iv工具变量
* xtabond2 得到的系数与标准误估计值与上文使用的 xtabond 的结果完全相同
*-16.13 偏差校正LSDV法
use mus08cigar.dta, clear
g t = year - 62
xtlsdvc lnc lnp lnpmin lny t, initial(ab) vcov(50) bias(3)
*-16.14 重复截面数据与组群分析
引自 第16章长面板与动态面板
2021-02-17 10:20:21
*-17 非线性面板
* 17.1 面板二值选择模型
* 对于面板数据，如果被解释变量为虚拟变量，则称为“面板二值选择模型”
* probit y x1 x2 x3, vce(cluster id) //混合probit回归
* logit y x1 x2 x3, vce(cluster id) //混合logit回归
* 17.2 面板二值选择模型的随机效应估计
* xtprobit y x1 x2 x3 //默认为随机效应probit
* xtlogit y x1 x2 x3 //默认为随机效应logit
* 17.3 面板二值选择模型的固定效应估计
* xtlogit y x1 x2 x3, fe
* 17.4 面板二值选择模型的stata实例
cd …data
use uprising.dta, clear
* 首先，估计固定效应面板logit估计
xtlogit uprising age pop temp sfamine sfamine1 sfamine2 sfamine3 relief relief1 relief2 relief3, fe nolog
est store FE
* 其次，估计混合回归logit模型
logit uprising age pop temp sfamine sfamine1 sfamine2 sfamine3 relief relief1 relief2 relief3 if dyn ~= 2, vce(cluster dyn) nolog
* 为了在固定和混合之间做选择，下面进行Hausman检验
* 但Hausman检验不允许在混合回归使用vce(cluster dyn)，故重新进行混合回归
qui logit uprising age pop temp sfamine sfamine1 sfamine2 sfamine3 relief relief1 relief2 relief3 if dyn ~= 2, nolog
est store POOLED
hausman FE POOLED
// 强烈拒绝混合回归原假设，应使用FE
xtlogit uprising age pop temp sfamine sfamine1 sfamine2 sfamine3 relief relief1 relief2 relief3 if dyn ~= 2, re nolog
est store RE
* 最后一行，LR检验强烈拒绝原假设 H0：ρ=0，认为应使用面板随机效应模型
* 应使用quadchk检验数值积分的稳健性
quadchk, noout //选择项 “noout” 表示不显示重估模型的结果
* 为在固定效应和随机效应间见选择，再次进行Hausman检验
hausman FE RE
* 可在0.05 水平上拒绝re的原假设，认为应使用固定效应模型
xtprobit uprising age pop temp sfamine sfamine1 sfamine2 sfamine3 relief relief1 relief2 relief3 if dyn ~= 2, re nolog
*-17.5 面板泊松回归
* 混合poisson回归：
* poisson y x1 x2 x3, vce(cluster id) irr //irr汇报发生率比
* 面板泊松回归：
* xtpoisson y x1 x2 x3, fe normal irr
// fe 表示固定效应，默认为re
// normal 表示在估计随机效应模型时，指定个体异质性服从正态分布，默认为gamma分布
// irr汇报发生率比
*-17.6 面板负二项回归
* 混合负二项回归
* nbreg y x1 x2 x3, vce(cluster id)
* 面板负二项回归
* xtnbreg y x1 x2 x3, fe irr
* 17.7 面板技术模型的stata示例
clear
cd "…data"
* 1) 数据描述
use mus18data.dta, clear
xtset id year
xtdes //面板数据描述
sum mdu lcoins ndisease female age lfam child
* 命令 xtsum 得到更详细的面板数据统计特征
* 2) 混合poisson回归，并使用聚类稳健标准误
poisson mdu lcoins ndisease female age lfam child, vce(cluster id) nolog
* 3) 随机效应的面板poisson回归
* 由于 xtpoisson, re 不提供聚类稳健标准误，故使用 vce(boot) 计算聚类自助标准误
xtpoisson mdu lcoins ndisease female age lfam child, re vce(boot, reps(50) seed(10101) nodots )
* reps(50) 表示重复抽样50次
* seed(10101) 表示以10101作为随机数的种子
* nodots 表示不显示重复抽样过程的点（dots）
// LR检验强烈拒绝 alpha=0的假设，拒绝混合，认为应使用随机效应
xtpoisson mdu lcoins ndisease female age lfam child, re normal nolog
* 前者对数似然函数高于后者，说明前者拟合的更好
quadchk, nooutput //对于某些系数，其相对差距超过了10-4
* 选择项 intpoints(#) 将数值积分的计算点数增加到16
xtpoisson mdu lcoins ndisease female age lfam child, re normal nolog intpoints(16)
* 然后再检验数值积分的精确度
quadchk, nooutput //所有系数的相对差距均小于10-4，故可认为稳健
* 4) 固定效应的面板poisson回归
* 1' 自助标准误
* vce(boot, reps(50) seed(10101) nodots)
* 2' xtpqml，提供聚类稳健标准误
xtpqml mdu lcoins ndisease female age lfam child, i(id)
// 必选项 i(id) 表示以变量id来确定面板单位
* 聚类稳健标准误大约是普通标准误的两倍，导致所有系数都不显著
* 5) 是否存在过离散
sum mdu, detail
dis r(Var)/r(mean)
* 6) 混合负二项回归，并使用聚类稳健标准误
nbreg mdu lcoins ndisease female age lfam child, vce(cluster id) nolog
// alpha 明显不含0，故使用nbreg可提高效率
* 7) 随机效应的面板负二项回归
xtnbreg mdu lcoins ndisease female age lfam child, nolog vce(boot, reps(50) seed(10101) nodots)
// 最后一行，LR检验强烈拒绝原假设，认为应使用随机效应的面板负二项回归
* 8) 固定效应的面板负二项回归
xtnbreg mdu lcoins ndisease female age lfam child, fe nolog vce(boot, reps(50) seed(10101) nodots)
// 固定效应的面板负二项回归可估计不随时间而变的变量系数
* 9) 使用 Hausman检验在fe 和 re中进行选择
qui xtnbreg mdu lcoins ndisease female age lfam child
est store re
qui xtnbreg mdu lcoins ndisease female age lfam child, fe
est store fe
hausman fe re
// 拒绝随机效应，支持固定效应
*-17.8 面板Tobit
* 混合Tobit回归、
* tobit y x1 x2 x3, ll(#) ul(#) vce(cluster id)
* 随机效应的面板Tobit回归
* xttobit y x1 x2 x3, ll(#) ul(#) tobit
* 进行混合Tobit 回归，并使用聚类稳健标准误
tobit med lcoins ndisease female age lfam child, ll(0) nolog vce(cluster id)
* 进行随机效应的面板Tobit回归
xttobit med lcoins ndisease female age lfam child, ll(0) nolog
* LR检验强烈拒绝原假设，认为存在个体效应，应使用re的面板Tobit回归
*-17.9 面板随机前沿模型
* 对于技术效率不随时间而变的固定效应模型，可直接用 reg 进行LSDV估计
* 对于随机效应模型，则
* xtfrontier lny lnx1 lnx2 lnx3, cost ti tvd
* lny 为对数形式的被解释变量
* lnx1 lnx2 lnx3 为对数形式的解释变量
* cost 为估计成本函数，默认估计生产函数
* ti 表示技术效率不随时间而变的 re模型
* tvd 表示估计随机效应的时变衰减模型
* 首先使用固定效应的LSDV法
use xtfrontier1.dta, clear
reg lnwidgets lnmachines lnworkers i.id, vce(cluster id)
* 上表显示，绝大多数的个体虚拟变量均很显著，证明存在个体效应
* 其次，估计技术效率不随时间而变的 re模型
xtfrontier lnwidgets lnmachines lnworkers, ti nolog
* 预测无效率项uit，将其记为uhat
predict uhat, u
list id uhat in 1/15
* 再次，估计随机效应的时变衰减模型
xtfrontier lnwidgets lnmachines lnworkers,tvd nolog
* 上表的eat表示η，估计值接近0，且p值为0.693，故可认为技术效率不随时间而变
* 由于该样本的平均时间跨度为10.4年，可能存在技术变迁，下面加入时间虚拟变量重新估计
xtfrontier lnwidgets lnmachines lnworkers i.t, tvd nolog
* 有三个时间虚拟变量在0.05水平上显著，时变系数也在0.1上显著引自 第17章非线性面板
2021-02-19 21:29:10
*-chp18 随机实验与自然实验
clear all
cd "…data"
* diff y, treat(varname) p(varname) cov(z1 z2) robust report test
use cardkrueger1994.dta, clear
gen gd = t * treated //定义交叉项 gd
reg fte gd treated t, r
* 使用 diff 估计
diff fte, t(treated) p(t) robust
* 引入快餐品牌虚拟变量
reg fte gd treated t bk kfc roys, r
* diff 引入其他解释变量
diff fte, t(treated) p(t) cov(bk kfc roys) robust
* 使用diff检验个变量在实验组与控制组的均值是否相等
set line 250
diff fte, t(treated) p(t) cov(bk kfc roys wendys) test引自 第18章随机实验与自然实验
2021-02-20 10:32:26
*-chp19 蒙特卡罗法与自助法
*-19.1 蒙特卡罗法的思想与用途
* 蒙特卡罗法实例：模拟中心极限定理
program onesample, rclass
drop _all
set obs 30
gen x = runiform()
sum x
return scalar mean_sample = r(mean)
simulate xbar = r(mean_sample), seed(101) reps(10000):onesample
hist xbar, normal xtitle("xbar from many sample")
*-19.3 蒙特卡罗法实例：服从卡方分布的扰动项
clear all
program chi2data, rclass
drop _all
set obs 30
gen double x = rchi2(1)
gen y = 1 + 2 *-x + rchi2(1) - 1
reg y x
return scalar t2 = (_b[x]-2) / _se[x]
return scalar r2 = abs(return(t2)) > invttail(28, 0.025)
simulate reject = r(r2), reps(1000) seed(101): chi2data
mean reject
*-19.4 蒙特卡罗积分
clear all
set obs 100
set seed 1234
gen double y = rnormal(0, 1)
gen double gy = exp(-exp(y))
sum gy
clear all
set obs 1000
set seed 1234
gen double y = rnormal(0, 1)
gen double gy = exp(-exp(y))
sum gy
*-19.5 最大模拟似然法与模拟矩估计
*-19.6 自助法的思想与用途
*-19.7 自助法的分类
*-19.8 使用自助法估计标准误
*-19.9 使用自助法进行区间估计
*-19.10 使用自助法进行假设检验
*-19.11 自助法的异质性（选读）
*-19.12 异方差情况下的自助法
*-19.13 面板数据与时间序列的自助法
*-19.14 自助法的Stata命令
* 1. 使用选择项“vce(bootstrap)”
* reg y x1 x2 x3, vce(boot, reps(400) seed(10101) nodots)
* boot 为 bootstrap 的简写
* reps(400) 表示 B=400
* seed(10101) 表示种子为10101，以保证可重复性
* nodots 表示执行中不显示抽样点
* 在执行 vce(bootstrap) 的Stata命令后，还可用以下命令获得自助法更多信息
* estat boot, all
* 2. 使用命令“bootstrap”
* 对于一般统计量使用自助法，可使用命令 bootstrap
* bootstrap explist [,options]:command
* "explist"提供了需要进行自助法估计量的统计量
* bootstrap, reps(400) seed(10101) nodots:reg y x1 x2 x3 //和vce(bootstrap)
* bootstrap _b _se, reps(400) seed(10101) nodots:reg y x1 x2 x3 //和vce(bootstrap) //求得回归系数和标准误
* 还可用"estat boot"获得自助法更多信息
* 3. 使用命令"bsample"
* 更复杂的情形可用bsample进行有放回抽样，然后用 simulate 进行重复抽样与计算
* bsample 50 得到样本容量为50的有放回抽样
*-19.15 使用自助法进行稳健的Hausman检验
cd "…data"
use grilic.dta, clear
program hausmantest, eclass
tempname b bols biv
reg lw iq s expr tenure rns smsa, r
matrix `bols' = e(b)
ivregress 2sls lw s expr tenure rns smsa (iq = med kww), r
matrix `biv' = e(b)
matrix `b' = `bols' - `biv'
ereturn post `b'
bootstrap _b, reps(400) seed(10101) nodots nowarn: hausmantest
test s iq expr tenure rns smsa _cons
dis chi2tail(1, 3.19)引自 第19章蒙特卡罗法与自助法
2021-02-21 21:47:53
*-chp20 平稳时间序列
*-20.1 时间序列的数字特征
*-20.2 自回归模型
*-20.3 移动平均模型
*-20.4 ARMA
*-20.5 自回归分布滞后模型
*-20.6 ARMA模型的Stata命令及实例
* 1. 自相关与偏相关
* corrgram y, lags(#)
* ac y, lags(#)
* pac y, lags(#)
* 2. ARMA
* arima y, ar(1/#) ma(1/#)
* ar(1/#) 表示第1至第#阶自回归
* ma(#) 表示第1至第#阶移动平均
* ARMA 的另一等价命令为：
* arima y, arima(#p, #d, #q)
* #p 表示自回归的阶数
* #q 表示移动平均的阶数
* #d 表示原序列{yi}需要经过几次差分才是平稳过程
* 检验残差项是否存在自相关：
* predict e1, res //计算残差，并命名为e1
* corrgram e1, lags(#) //检验残差是否存在第1至第#阶自相关的Q检验
* 3. ADL 与 ARMAX
* ARMAX 的 Stata 命令：
* arima y x1 x2 x3, ar(#) ma(#)
* 对非平稳序列 logpe 建立ARMA模型
* 首先定义其一阶差分为 d_logpe
cd "…data"
use pe.dta, clear
tsset year
g d_logpe = d.logpe
* 计算前10阶自相关与偏相关系数
corrgram d_logpe, lags(10)
ac d_logpe, lags(10)
pac d_logpe, lags(10)
* 第4阶自相关与偏相关系数显著不为0
* 考虑AR(4) 与 MA(4) 模型
arima d_logpe, ar(1/4) nolog
estat ic
* 检查其残差项是否存在自相关：
predict e1, res
corrgram e1, lags(10) //结果表明，可以接受无自相关原假设，检验至第10阶自相关都不显著
* 其次，估计MA(4)模型
arima d_logpe, ma(1/4) nolog
estat ic
predict e2, res
corrgram e2, lags(10) //也不存在自相关
* 根据信息准则，AR(4) 模型由于MA(4)模型由于MA
* 将AR(4) 与 MA(4) 模型中的第1阶与第3阶略去
arima d_logpe, ar(2 4) nolog
estat ic
arima d_logpe, ma(2 4) nolog
estat ic
* 20.7 误差修正模型
* 20.8 MA(∞)与滞后算子
* 20.9 向量自回归过程
* 20.10 VAR的脉冲响应函数
* 20.11 预测误差的方差分解
* 20.12 格兰杰因果检验
* 20.13 面板格兰杰因果检验
* 20.14 VAR的Stata命令及实例
* varsoc x y z, maxlag(#)
* 此命令用来计算不同滞后期的信息准则
* soc 表示 Selection-Order Criteria
* maxlag(#) 表示最大滞后期，默认值为4
* 在未使用 maxlag(#) 情况下，如果恰好选择最优滞后4期，则通过 maxlag(#) 增加最大滞后期数。
* varbasic x y z, lags(numlist) irf fevd
* 估计VAR模型的便捷名
* lags(numlist) 表示滞后阶数，默认为 lags(1 2)，即滞后二阶
* irf 表示画（未正交化）脉冲响应图
* fevd 话预测方差分解图，默认为 orif 画正交化脉冲响应图
* var x y z, lags(numlist) dfk small exog(w1 w2)
* 选择项 lags(numlist) 表示滞后阶数，默认为 lags(1 2)，即滞后二阶
* 如果样本容量较小，可使用选择项 dfk 进行自由度调整
* small 表示显示小样本的 t 或 F 统计量，而非大样本的标准正态或卡方统计量
* exog(w1 w2) 表示在VAR模型中引入外生变量w1, w2
* 在使用命令 var 完成 VAR模型的估计之后，Stata提供一系列“估计后命令”来检验VAR的性质与假设。
* varlmar //估计VAR后，对残差是否存在自相关进行LM检验
* varnorm //估计VAR后，检验残差是否服从正态分布
* varstable, graph //估计VAR后，通过特征值检验该VAR是否为平稳过程。如果所有特征值都在单位圆内，则为平稳过程。graph显示特征值的几何分布图。
* varwle //估计VAR后，对每个方程以及所有方程的各阶系数的联合显著性进行wald检验，其中 wle 表示 Wald lag-exclusion statistics.
* vargranger //估计VAR后，进行格兰杰因果检验
* irf create irfname, set(filename) step(#) replace
* 估计VAR后，将有关脉冲响应的结果为irfname（可自行命名）
* set(filename) 表示建立脉冲文件 filename，使之成为当前的脉冲文件(make filename active)，并将脉冲结果 irfname 存入此脉冲文件 filename （若未使用选择项 set(filename) ,则将脉冲响应结果存入当前的脉冲文件）
* step(#) 表示考察#期的脉冲响应函数，默认值为 step(8)
* replace 表示替代已有的同名脉冲响应结果 irfname（如果有）
* 一个脉冲文件 filename 可存储多个脉冲响应结果 irfname
* irf graph irf, impulse(varname) response(varname) noci
* 画脉冲响应图（未正交化）。
* impulse(varname) 用于指定脉冲变量
* response(varname) 用来指定反应变量，默认画出所有变量的脉冲响应图
* noci 表示不画置信区间，默认画置信区间
* irf graph cirf, impulse(varname) response(varname) //累积脉冲响应图（未正交化）
* irf graph orif, impulse(varname) response(varname) // 正交化的脉冲响应图
* irf graph coirf, impulse(varname) response(varname) //正交化的累积脉冲响应图
* irf graph fevd, impulse(varname) response(varname) //预测方差分解图
* 如果将以上的 irf graph 改为 irf table，则将相应信息列表而非画图
* fcast compute prefix, step(#)
* 估计VAR后，计算被解释变量未来#期的预测值，并把预测值赋予被解释变量假设前缀 prefix （自行确定）的变量名
* fcast graph varlist, observed
* 运行命令 fcast compute 后，将变量 varlist 的预测值画图，其中 observed 表示与实际观测值相比较。
use varexample.dta, clear
tsline inflation unrate fedfunds, lpattern("-" "-") xline(169)
* xline(169) 表示在x轴上 t = 169 的位置画一条线
* 预留近10年数据作为“样本外预测”
* 图20.3中的垂直线为样本内与样本外的分界线
* 样本内外，失业率的波动范围大体相同。
sum inflation unrate fedfunds if date <= tq(2002q1)
sum inflation unrate fedfunds if date > tq(2002q1)
* 为了估计VAR，先根据信息准则确定 VAR模型 的阶数
varsoc inflation unrate fedfunds if date < tq(2002q1), maxlag(13)
* maxlag(13) 表示最多滞后13阶
* LL - 对数似然函数
* LR - 似然比检验，即对最后一阶系数的联合显著性进行似然比检验
* df 与 q 表示似然比统计量的自由度与p值
* FPE 度量向前一期预测的均方误差
* 上表显示，不同信息准则所选择的滞后阶数并不一致(星号)
* 最简洁的SBIC - 2阶 （过于简洁）
* HQIC - 3阶
* FPE - 9阶
* AIC - 11阶 （损失较多样本容量）
* LR - 12阶
* SBIC 与 HQIC 提供了对真实滞后阶数的一致估计，而FPE与AIC可能高估滞后阶数
* 作为折中，选择滞后4阶。
* 滞后4阶的VAR模型依然不能保证扰动项为白噪声，故本例选择滞后5阶，刚好能保证扰动项为白噪声。
var inflation unrate fedfunds if date <= tq(2002q1), lags(1/5)
* 即使滞后5阶，此VAR模型依然包含了48个参数，系数如此之多，以致于无法解释经济含义
* 因此，在实证论文中，甚至不汇报 VAR的回归系数，而主要汇报脉冲响应函数，预测方差分解与格兰杰检验。
* 在上述命令 var 中，由于样本容量为164，故没有使用选择项 dfk 或 small 进行小样本自由度调整。
* 检验各阶系数的联合显著性
varwle
// 虽然单一方程的某阶系数不显著，但作为三个方面的整体，各阶系数均高度显著。
* 下面检验残差是否为白噪声，即残差是否存在自相关
varlmar
* 结果显示，可接受残差“无自相关”原假设，即认为扰动项为白噪声
* 进一步检验此 VAR系统是否稳定（为平稳过程），结果见下表。
varstable, graph
* 所有特征值均在单位圆内，故VAR系统文档，但有两个根十分接近单位圆，意味着有些冲击有较强的持续性
* 进一步检验VAR模型的残差是否服从正态分布
varnorm
* 绝大多数检验结果均在 5% 的显著性水平上拒绝这三个变量的扰动项服从正态分布的原假设
* 尽管扰动项不服从正态对VAR的影响不大，但残差项的非正态性暗示模型可能偏离了真实的数据生成过程（DGP），并使得变量未来值的预测区经变得不可信
* VAR模型的用途之一是预测。
* 下面，预测未来40个季度（10年）的变量取值，分别记为 f_fedfunds, f_inflation 与 f_unrate， 然后画图
fcast compute f_, step(40)
fcast graph f_inflation f_unrate f_fedfunds, observed lpattern("_")
* observed 表示显示变量的实际观测值
* lpattern("_") 表示以虚线来表示变量的预测值
* VAR无法预测经济危机，而且预测时期越长，精度越低
* 除了预测外，还想知道对某变量的冲击会对该变量及其他变量产生怎样的动态影响，需要用到正交化的脉冲响应函数，但正交化的脉冲响应函数依赖于变量的排序。
* 为此分别考察变量间的格兰杰因果检验与交叉相关图。
* 首先，考察三个变量间的格兰杰因果检验：
vargranger
// unrate 是 inflation 的格兰杰原因
// fedfunds 不是 inflation 的格兰杰原因
// 拒绝 unrate、fedfunds 都不是 inflation 的格兰杰原因
* granger因果关系并未给出唯一的变量作用次序，为此，进一步考察交叉相关图
xcorr inflation unrate if date <= tq(2002q1), name(iu)
xcorr inflation fedfunds if date <= tq(2002q1), name(ir)
xcorr unrate fedfunds if date <= tq(2002q1), name(ur)
graph combine iu ir ur
* 左上部分显示，inflation 与滞后10个季度的unrate最相关
* 右上部分显示，inflation与滞后1个季度的fedfunds最相关
* 左下部分显示，unrate与提前8个季度的fedfunds最相关
* 如果无法看出交叉相关系数的最大值，可通过 xcorr, table 将交叉相关系数列表
xcorr inflation unrate if date <= tq(2002q1), table
* 因此，在计算脉冲响应函数与预测方差分解时，仍需变换次序以考察结果的稳健性：
irf create iuf, set(macrovar) step(20)
* irf created 计算所有与脉冲响应函数有关的变量与统计量
* 将计算结果命名为 iuf，存入新建立的脉冲文件 macrovar.irf，将此脉冲文件激活为当前文件，并更新
* 选择项 step(20) 意味着计算20期的脉冲响应函数，默认为step(8)
* 此命令中，由于没有选择项 order(varlist) 来指定变量次序，故默认为上文使用命令 var 进行估计的变量次序，即 inflation, unrate, fedfunds
* 脉冲文件其实是普通的stata数据文件，不过扩展名为irf
des using macrovar.irf
* 显示单个脉冲响应图（inflation 对 unrate 正交脉冲的响应）
irf graph oirf, yline(0) i(unrate) r(inflation)
* 变量 inflation 的预测误差方差分解
irf table fevd, r(inflation) noci
// r(inflation) 表示以 inflation 为响应变量，noci 表示不显示置信区间
// 对 inflation 向前一季度的预测，预测方差完全来自 inflation 本身。
// inflation 主要受自身的影响，变量 unrate 与 fedfunds 的作用很小。
* 把上述结果用图表达出来
irf graph fevd, r(inflation)
* 考察 unrate 的方差分解
irf table fevd, r(unrate) noci
// 48% 来自 inflation
* 下面考察 fedfunds 的预测方差分解：
irf table fevd, r(fedfunds) noci
// 少部分来自自身
* 预测方差分解图
irf graph fevd
* 以上脉冲响应与方差分解结果依赖于次序变量
* 故将次序变为交叉相关图所推荐的次序 inflation, fedfunds, unrate，以考察结果的稳健性
* 为此，建立新的脉冲响应结果，命名为 ifu
irf create ifu, order(inflation fedfunds unrate) step(20)
* 比较在以上两种变量排序下， inflation 对 fedfunds 脉冲的响应
irf graph oirf, i(fedfunds) r(inflation) yline(0) noci
* 比较在两种变量排序下， inflation 的预测方差来自 fedfunds 的比重
irf graph fevd, i(fedfunds) r(inflation) noci
*-20.15 季节调整
* 1.回归法
use turksales.dta, clear
tsline sales
gen quarter = quarter((dofq(t)))
tab quarter, gen(q)
* 然后，以第1季度为参照值，回归
reg sales q2 - q4
* 为了获得经季节调整的序列，下面使用 predict 来获得上述回归的残差项
predict sales_sa, r
* ols 残差项平均值一定为0，故需要把原序列的均值加回去，并记季节调整序列为sales_sa_reg
sum sales
gen sales_sa_reg = sales_sa + r(mean)
tsline sales_sa_reg sales, lpattern("_")
* 2.移动平均比率法
tssmooth ma sales1 = sales, window(2 1 1)
// 上述命令表示对序列 sales 移动平均，并将结构记为 sales1
// window(2 1 1) 表示在移动平均时，包括过去2项，当期1项，以及未来1项
tssmooth ma sales2 = sales, window(1 1 2)
// window(1 1 2) 表示在移动平均时，包括过去1项，当期1项，以及未来2项
gen sales_tc = (sales1 + sales2) / 2
* 将原序列减去趋势循环序列 TC
gen sales_si = sales - sales_tc
* 通过循环语句 forvalues 来计算每一季度的平均值，即季节因子
forvalues i=1/4{
qui sum sales_si if q`i' == 1
scalar s`i' = r(mean)
* 将四个季度的平均值进行平均
scalar s_mean = (s1 + s2 + s3 + s4) / 4
* 以循环语句将每个季节因子减去季节因子的平均值，得到标准化的季节因子，然后将原序列减去标准化的季节因子，即得到季节调整序列。
gen sales_sa_ma = 0
forvalues i=1/4{
scalar s`i' = s`i' - s_mean
replace sales_sa_ma = sales - s`i' if q`i' == 1
* 最后，看一下移动平均比率法的季节调整序列与原序列的时间序列图
tsline sales_sa_reg sales_sa_ma sales, lpattern("_" "-")
gen sales_sim = sales / sales_tc
forvalues i=1/4{
qui sum sales_sim if q`i' == 1
scalar sm`i' = r(mean)
scalar s_meanm = (sm1 * sm2 * sm3 * sm4)^0.25
gen sales_sa_mam = 0 //初始化季节调整序列为0，以便后面更新
forvalues i=1/4{
scalar sm`i' = sm`i' / s_meanm
replace sales_sa_mam = sales / sm`i' if q`i' == 1
* 对比加法模型和乘法模型的季节调整结果
tsline sales_sa_ma sales_sa_mam, lpattern("_")
* 3.X12方法
db sax12im
tsline sales_sa_ma sales_d1l, lpattern("-")引自 第20章平稳时间序列
2021-02-22 21:26:56
*-chp 21 单位根所带来的问题
*-21.1 非平稳序列
*-21.2 ARMA的平稳性
*-21.3 VAR的平稳性
*-21.4 单位根所带来的问题
* 1. 自回归系数的估计值向左偏向于0
* 2. 传统的t检验失效
* 3. 两个相互独立的单位根变量可能出现伪回归或伪相关
clear all
drop _all //删去内存中已有的数据
set obs 10000
set seed 1234
gen u = rnormal() //产生服从标准正态分布的扰动项ui
gen y = sum(u)
set seed 12345
gen v = rnormal() //产生服从标准正态分布的扰动项vi
gen x = sum(v) //定义随机游走
reg y x
gen t = _n //定义时间变量t
line y x t, lpattern(dash)
* 如何避免伪相关或伪回归？
* 1. 先对变量做一阶差分，然后再回归。
* 2. 协整，但首先必须对是否存在单位根进行检验
*- 21.5 单位根检验与平稳性检验
* 1. DF检验
* 2. ADF检验
* 是否应该带常数项或时间趋势项，主要应从理论上考虑。
* 在进行ADF检验时，如果确定滞后阶数p是个问题。
* dfuller y, lags(p) regress noconstant drift trend
* lags(p) 表示包含p阶滞后差分项，默认为lags(0),对应于DF检验
* regress 表示显示回归结果
* 3. PP检验
* pperron y, noconstant trend regress lags(#)
* lags(#) 指定N-W滞后阶数
* 其他选择项含义与ADF检验(dfuller)相同
* 4. DF-GLS单位根检验
* ADF 与 PP检验的共同缺点是检验功效低
* 两步检验：
* 1. 用GLS估计原序列{yi}的常数项与事件趋势项，计算去势后的序列。
* 2. 对{yid}使用ADF检验，这个检验被称为 DF-GLS检验，是目前最有功效的单位根检验。
* dfgls y, notrend maxlag(#) ers
* 5. KPSS平稳性检验
* kpss y, notrend maxlag(#)
* 6. 单整阶数的确定
*-21.6 单位根检验的Stata实例
clear
cd "…data"
use macro_swatson.dta, clear
line inf quarter
dfuller inf
* DF统计量为-3.718 < -3.488， 可在0.01水平上拒绝“存在单位根”的原假设。
* DF检验中的扰动项可能存在自相关，古窑考虑更高阶的ADF检验。
* 首先，计算Schwert建议的最大滞后阶数。
di 12 * (167/100)^(1/4)
* 令 p_hat = 13，进行 ADF检验
dfuller inf, lags(12) reg
* 上表显示，最后一阶滞后项（L12D）在0.05 水平上并不显著。
* 依次令 p_hat = 12, 11, 10，进行 ADF检验，最后一阶滞后项仍不显著。
* 令 p_hat = 9，再进行ADF检验
dfuller inf, lags(8) reg
* 最后一阶滞后项在0.05 水平上显著不为0。
* ADF统计量 Z(t)显示，无法在0.05水平上拒绝“存在单位根”的原假设。
* 可认为inf存在单位根
* 解决 DF检验中扰动项自相关问题的另一方法为pp检验
pperron inf
* PP检验 也是左边单侧检验，而检验统计量均小于0.05临界值，故可在0.05水平上拒绝“存在单位根”的原假设。
* 进行最有功效的检验，DF-GLS检验
dfgls inf
// 上表显示，从1阶至13阶滞后， 均无法在 0.05 水平拒绝存在单位根假设
* 将原假设变为平稳序列，进行KPSS检验
kpss inf, notrend
dfgls dinf
* 最优滞后阶数介于 2 与 7 之间
* 在此区间，检验统计量均小于 0.05 的临界值
kpss dinf, notrend
* 21.7 面板单位根检验
* 1. LLC检验
* LLC检验假设不存在截面相关
* xtunitroot llc y, trend noconstant demean lags(#) lags(aic #) lags(bic #) lags(hqic #)
cd "…data"
use pennxrate.dta, clear
xtunitroot llc lnrxrate if g7, lags(aic 10)
* lags(aic 10) 表示将差分滞后项的最大滞后阶数设为10，并根据 AIC信息准则 选择最优滞后阶数
* 使用 demean 缓解此截面相关
xtunitroot llc lnrxrate if g7, lags(aic 10) demean
* 2. HT检验
* xtunitroot ht y, trend noconstant demean
use pennxrate.dta, clear
xtunitroot ht lnrxrate, demean
* 3. Breitung 检验
* xtunitroot breitung y, trend noconstant demean robust lags(#)
xtunitroot breitung lnrxrate if oecd, robust
* 4. IPS检验
* xtunitroot ips y, trend demean lags(#) lags(aic #) lags(bic #) lags(hqic #)
xtunitroot ips lnrxrate if oecd, demean
xtunitroot ips lnrxrate if oecd, lags(aic 8) demean
* 5. 费雪式检验
* xtunitroot fisher y, dfuller pperron demean lags(#)
xtunitroot fisher lnrxrate, dfuller drift lags(2) demean
* 6. Hadri LM 检验
* xtunitroot hadri y, trend demean robust kernel (bartlett #) kernel(perzen #) kernel(quadraticspectral #)
xtunitroot hadri lnrxrate if oecd, kernel(bartlett 5) demean
* 21.8 协整的思想与初步检验
* 21.9 Beveridge-Nelson 分解公式
* 21.10 协整的定义与最大似然估计
* 21.11 协整分析的stata实例
* vecrank y1 y2…yn, lags(#) //默认包括常数项，但不包括时间趋势
* vecrank y1 y2…yn, lags(#) trend(none) //不包括常数项或时间趋势
* vecrank y1 y2…yn, lags(#) trend(trend) //包括常数项或时间趋势
* 在做完协整秩检验，并确定 h≥1后，就可以对 VECM模型进行最大似然估计。
vec y1 y2 … yn, lags(#) rank(#) //默认设置包括常数项，但不包括时间趋势
vec y1 y2 … yn, lags(#) trend(none) //不包括常数项，也不包括时间趋势
vec y1 y2 … yn, lags(#) trend(trend) //包括常数项，也包括时间趋势
irf graph irf, impulse(varname) response(varname) noci
irf graph cirf, impulse(varname) response(varname)
irf graph oirf, impulse(varname) response(varname)
irf graph coirf, impulse(varname) response(varname)
irf graph fevd, impulse(varname) response(varname)
fcast graph varlist, observed
cd "…data"
use mpyr.dta, clear
line logmr logy year, lpattern("l" "_")
gen _logv = -logv
gen r10 = r / 10
tsline _logv r10, lpattern("_")
vecrank logmr logy r, trend(trend) max
* 其次，检验改系统所对应的VAR表示法的滞后阶数
varsoc logmr logy r
* 使用 Johansen 的 MLE方法估计该系统的向量误差修正模型(VECM):
vec logmr logy r, lags(2) rank(1)
* 用 OLS 估计长期均衡关系 (EG-ADF两步法)：
reg logmr logy r
* 检验 VECM模型 是否存在自相关，如果存在自相关，则预示着需要增加滞后阶数。
qui vec logmr logy r, lags(2) rank(1)
veclmar
// 可以接受无自相关假设
vecnorm
* 检验此 VECM系统是否稳定
vecstable, graph
* 考察 VECM模型的正交化脉冲响应函数（将脉冲文件命名为money）
irf create money, set(money) step(10) replace
irf graph oirf
qui vec logmr logy r if year < 1980, lags(2) rank(1)
fcast compute f_, step(10)
fcast graph f_logmr f_logy f_r, observed lpattern("_")引自 第21章单位根与协整
2021-02-23 08:23:32
*-chp22 自回归条件异方差模型
*-22.1 条件异方差模型的例子
*-22.2 ARCH模型的实质
*-22.3 ARCH模型的MLE估计
*-22.4 GARCH模型
*-22.5 何时使用ARCH或GARCH模型
*-22.6 ARCH 与 GARCH模型的扩展
*-22.7 ARCH 与 GARCH 的stata命令及实例
* arch y x1 x2, arch(1/3) // ARCH(3)
* arch y x1 x2, arch(1) garch(1) // GARCH(1 ,1)
* arch y x1 x2, ar(1) ma(1) arch(1) garch(1) //带ARMA(1, 1) 的 GARCH(1, 1)
* arch y x1 x2, arch(1) dist(1) //ARCH(1)，将扰动项服从 t 分布
* arch y x1 x2, arch(1) het(z1 z2) //ARCH(1)，将z1, z2加入条件方差方程
* arch y x1 x2, arch(1) garch(1) tarch(1) //GARCH(1, 1)加上TARCH(1)
* arch y x1 x2, earch(1) egarch(1) //EGARCH(1 ,1)
* arch y x1 x2, arch(1/3) archm //ARCH(3) 加上ARCH-M
* 首先，看日收益率的时间趋势
cd "…data"
use sp500.dta, clear
line r t
varsoc r, maxlag(8)
* 上表显示，大多数准则均选择AR(5)模型
reg r L(1/5).r
* 对OLS残差是否存在ARCH效应进行LM检验
estat archlm, lags(1/5)
* 通过画图更直观地考察 ols的残差平方是否存在自相关
predict e1, res
g e2 = e1^2
ac e2
pac e2
corrgram e2, lags(10)
varsoc e2
arch r L(1/5).r, arch(1/3) nolog
* 所有准则均显示，应考虑ARCH(3)模型
arch r L(1/5).r, arch(1/3) nolog
* ARCH项很显著，使用更简洁的 GARCH(1, 1)模型
arch r L(1/5).r, arch(1) garch(1) nolog
* 考虑到“好”“坏”消息可能不对称，因此，加上TARCH项：
arch r L(1/5).r, arch(1) garch(1) tarch(1) nolog
// tarch 项显著，存在不对称效应
arch r L(1/5).r, arch(1) archm nolog
arch r L(1/3).r, earch(1) egarch(1) nolog
kdensity r, normopt(lpattern("-"))
quietly var r, lags(1/5)
varnorm
arch r L(1/5).r, arch(1) garch(1) dist(t) nolog
* 最后，对GARCH(1, 1)模型的条件方差进行预测：
qui arch r L(1/5).r, arch(1) garch(1) nolog
predict, variance
line h t
*-22.8 多维 GARCH模型
* mgarch model eq…eq, arch(#) garch(#) constraints(#) het(varlist) dist(t) robust noconstant
* e.g.
cd "…data"
use irates4.dta, clear
* tbil 和 bond 两个变量的一阶差分，进行VAR(1)建模，并带ARCH(1)扰动项。
* 下面这几个 mgarch 命令不知道为啥都没运行出来
mgarch dvech(D.bond D.tbill = LD.bond LD.tbill), arch(1) nolog
mgarch dvech(D.bond D.tbill = LD.bond LD.tbill, noconstant) (D.tbill = LD.tbill, noconstant), arch(1) nolog
mgarch ccc(D.bond D.tbill = LD.bond LD.tbill), arch(1) nolog
mgarch dcc(D.bond D.tbill = LD.bond LD.tbill), arch(1) nolog
mgarch vcc(D.bond D.tbill = LD.bond LD.tbill), arch(1) nolog引自 第22章自回归条件异方差模型
2021-02-24 13:36:58
*-chp23 似不相关回归
*-23.1 单一方程估计与系统估计
*-23.2 似不相关回归的假定
*-23.3 SUR 的 FGLS 估计
*-23.4 SUR 的假设检验
*-23.5 似不相关回归的stata命令及实例
* sureg (depvar1 varlist1) (depvar2 varlist2) … (depvarn varlistn), isure corr
* isure 表示迭代至收敛
* corr 汇报对无同期相关的检验
* 如果要检验联合假设“第1个方程中变量x的系数与第2个方程中变量z，以及第3个方程中变量w的系数相等”，可使用：
* test ([depvar1]x = [depvar2]z) ([depvar1]x = [depvar3]w)
* 如果要在满足以上两个约束条件的情况下进行SUR估计，则：
* contraint1 [depvar1]x = [depvar2]z （定义第1个约束条件）
* contraint2 [depvar1]x = [depvar2]w （定义第1个约束条件）
* sureg (depvar1 varlist1) (depvar2 varlist2) … (depvarn varlistn),c(1)(2)
cd "…data"
use hsb2.dta, clear
reg3 (read female ses schtyp socst)(math female ses schtyp science), ols
est store OLS
* reg3 为 三阶段最小二乘法（3SLS），加上 "ols" 后则进行单一方程OLS估计。
* 由于被解释变量 read 与math为同一学生的阅读与数学成绩，故这两个方程的扰动项在理论上很可能存在相关性。
sureg (read female ses schtyp socst) (math female ses schtyp science), corr
est store SUR
* 可在1%水平上拒绝扰动项相互独立假设
sureg (read female ses schtyp socst) (math female ses schtyp science), i nolog
est store SUR_i
local m "OLS SUR SUR_i"
local mt "OLS SUR SUR_i"
esttab `m' `s', nogap compress replace ///
b(%6.3f) s(N r2, fmt(%5.0f %6.3f)) ///
se(%6.3f) ///
keep(`keep') drop(`drop') ///
star(* 0.1 ** 0.05 *** 0.01) ///
mtitle(`mt') ///
addnotes("*** 1% ** 5% * 10%")
* 检验系数是否相等
qui sureg (read female ses schtyp socst) (math female ses schtyp science)
test ([read]ses = [math]ses)
* p值为0.91，接受假设
* 没跑出来
cap contraint1 [read]ses = [math]ses
cap sureg (read female ses schtyp socst) (math female ses schtyp science), c(1)
*-23.6 变系数面板数据的SUR估计
use mus08cigar.dta, clear
reshape wide lnc lnp lnpmin lny, i(year) j(state)
set line 250
list in 1/2
sureg (lnc1 lnp1 lnpmin1 lny1) (lnc2 lnp2 lnpmin2 lny2) (lnc3 lnp3 lnpmin3 lny3) (lnc4 lnp4 lnpmin4 lny4) (lnc5 lnp5 lnpmin5 lny5) (lnc6 lnp6 lnpmin6 lny6) (lnc7 lnp7 lnpmin7 lny7) (lnc8 lnp8 lnpmin8 lny8) (lnc9 lnp9 lnpmin9 lny9) (lnc10 lnp10 lnpmin10 lny10), corr引自 第23章似不相关回归
2021-02-26 16:17:05
*-chp24 联立方程模型
*-24.1 联立方程模型的结构式与简化式
*-24.2 联立方程模型的识别
*-24.3 单一方程估计法
*-24.4 三阶段最小二乘法
*-24.5 三阶段最小二乘法的stata实例
cd "…data"
* reg3 (depvar1 varlist1) (depvar2 varlist2) … (depvarN varlistN), ols 2sls sure ireg3 exog(varlist) endog(varlist) inst(varlist)
* ireg3 迭代式 3sls估计
* exog 额外的外生变量
* endog 额外的内生变量
* inst 额外的工具变量
use klein.dta, clear
reg3 (consump wagepriv wagegovt) (wagepriv consump govt capital1), ols
est store OLS
* 加上 sls 后，表示单一方程估计
reg3 (consump wagepriv wagegovt) (wagepriv consump govt capital1), 2sls
est store Two_SLS
* 加上 2sls 后，表示单一2SLS方程估计
reg3 (consump wagepriv wagegovt) (wagepriv consump govt capital1), first
est store Three_SLS
* 加上 first 后，表示显示第一阶段回归的结果
reg3 (consump wagepriv wagegovt) (wagepriv consump govt capital1), ireg3
est store Three_SLS_iter
local m "OLS Two_SLS Three_SLS Three_SLS_iter"
local mt "OLS Two_SLS Three_SLS Three_SLS_iter"
esttab `m' `s', nogap compress replace ///
b(%6.3f) s(N r2, fmt(%5.0f %6.3f)) ///
se(%6.3f) ///
keep(`keep') drop(`drop') ///
star(* 0.1 ** 0.05 *** 0.01) ///
mtitle(`mt') ///
addnotes("*** 1% ** 5% * 10%")
*-24.6 结构VAR
*-24.7 SVAR的stata实例
/* savr y1 y2 y3, aconstraints(constrints_a) aeq(matrix_aeq) acns(matrix_acns) bconstraints(constrints_a) beq(matrix_aeq) bcns(matrix_acns) varconstraints(constrints_v) lags(numlist)
aconstraints(constrints_a) constrint所定义的约束
aeq(matrix_aeq) matrix定义的约束
acns(matrix_acns) matrix定义的跨参数约束
bconstraints(constrints_a) 以下表示三种对矩阵B施加的跨参数约束
beq(matrix_aeq)
bcns(matrix_acns)
varconstraints(constrints_v) 对简化VAR的参数进行约束
lags(numlist) 简化VAR的阶数 */
use lutkepohl2.dta, clear
matrix A = (1, 0, 0 \ ., 1, 0 \ ., ., 1)
matrix B = (., 0, 0 \ 0, ., 0 \ 0, 0, .)
matrix list A
matrix list B
varsoc dln_inv dln_inc dln_consump //确定VAR模型的滞后阶数
svar dln_inv dln_inc dln_consump, aeq(A) beq(B) nolog
* 估计 SVAR 的主要兴趣通常在于考察其脉冲响应函数。
irf create germany, set(germany)
irf graph sirf
* sirf 结构脉冲响应函数，等价于 正交化的脉冲响应函数
* 考察 SVAR模型 的预测误差方差分解。
irf graph sfevd
* sfevd 表示 结构预测误差方差分解。
* 以上为恰好识别的 SVAR模型。
* 注意到，矩阵A的a21元素不显著，
matrix A = (1, 0, 0 \ 0, 1, 0 \ ., ., 1)
matrix B = (., 0, 0 \ 0, ., 0 \ 0, 0, .)
matrix list A
svar dln_inv dln_inc dln_consump, aeq(A) beq(B) nolog
* 长期SVAR命令：
* svar y1 y2 y3, lrconstraints(constraints_lr) lreq(matrix_lreq) lrcns(matrix_lrcns) varconstraints(constraints_v) lags(numlist)
* 对长期效应矩阵C作如下约束
use m1gdp.dta, clear
* 首先，定义矩阵C
matrix C = (., 0\ ., .)
matrix list C
* 其次，根据信息准则，确定简化VAR的滞后阶数。
varsoc ln_gdp ln_m1
* 上表显示，大多数信息准则支持 VAR(3)。下面，估计长期SVAR模型
svar ln_gdp ln_m1, lreq(C) lags(1/3) nolog
* 上表显示，此长期SVAR模型为恰好识别，且C矩阵中的所有元素都很显著。引自 第24章联立方程模型
2021-02-27 15:32:08
*-chp25 非线性回归与门限回归
*-25.1 非线性最小二乘法
*-25.2 非线性回归的stata命令及实例
* NLS 的stata命令：
* nl (depvar=<sexp>) [if] [in] [weight] [, options]
* nl sexp_prog : depvar [varlist] [if] [in] [weight] [, options]
* nl func_prog @ depvar [varlist] [if] [in] [weight], {parameters(namelist)|nparameters(#)} [options]
* nl (y = {beta1} + {beta2} * exp({beta3 = 1} * x)), robust
* {} 表示待估参数
* {beta3 = 1} 表示β3 的初始值为1
cd "…data
use consumption_china.dta, clear
nl (c = {beta1} + {beta2} * y^{beta3 = 1}), nolog
* 使用稳健标准误重新估计
nl (c = {beta1} + {beta2} * y^{beta3 = 1}), r nolog
*-25.3 门限回归
*-25.4 面板数据的门限回归
*-25.5 门限回归的计算机操作引自 第25章非线性回归与门限回归
2021-03-01 14:40:13
