stata---门限回归

参见视频：
http://www.tudou.com/programs/view/D2PLmztpI1k/
http://www.tudou.com/programs/view/CWmhnWUou7w/

Stata门限模型的操作和结果详细解读：
http://www.wenkuxiazai.com/doc/d6057d60c281e53a5802ffd7-2.html

*-模拟数据的生成：单一门槛
clear
set obs 100    //100个观察值
set seed 123456 
gen e = invnorm(uniform())   //干扰项，服从正态分布 e~N(0,1)
gen x = 3 * invnorm(uniform()) //x~N(0, 3^2)
gen t = _n    //[1,100]
gen y = .     //空
tsset t       //tsset指令是时间序列数据的估计命令
replace y = 1 + 2*x + e
replace y = 1 + -2*x + e if t>50
save xtthres_siml.data, replace

将数据生成一下

结果：

scatter y x   //绘制散点图

reg y x       //回归
est store full  //存储回归结果

x前面的系数是-0.5669306
p=0.005，p越大，错误的概率越大，通常p不能大于0.005
r^2=0.0769 越接近0，越不好；越接近1，越好

reg y x if t<=50
est store left

reg y x if t>50
est store right

*-虚拟变量回归：Chow Test的基本思想
*-如果知道断点 在t=50
gen d = (t>50)  //t=50
gen xd = x*d
reg y x xd
est store strue50

x前面的系数是1.949462，很接近2
xd前面的系数是-3.980677，加上x前面的系数很接近-2

*-如果不知道断点
gen d10 = (t>10)  //假设t=10
gen xd = x*d10
reg y x xd10
est store strue10
gen d20 = (t>20)  //假设t=30
gen xd = x*d20
reg y x xd20
est store strue20
gen d30 = (t>30)  //假设t=30
gen xd = x*d30
reg y x xd30
est store strue30
gen d40 = (t>40)  //假设t=40
gen xd = x*d40
reg y x xd40
est store strue40

local m "full strue10 strue20 strue30 strue40 strue50"
est_table 'm', mtitle('m')

*-网格搜索结构突变点
gen rss = .  //残差列表
forvalues i=5/95{ //i从5到95，不能从1到100，否则左侧或右侧没有数据
    gen di = (t>`i')
    gen xdi = x*di
    qui reg y x xdi
    qui replace rss = e(rss) in `i'  //替换rss中第i个值
    drop di xdi
}
sort rss
list rss t in 1/5  //输出最小的前5个结果

最佳门限值：t=50，残差平方和为74.68341

*-图示
line rss t, xline(50, lp(dash) lw(thick)) ///
                  lw(thick) xlabel(0(10)100) sort

以上代码总结如下：

*-模拟数据的生成：单一门槛
clear
set obs 100    //100个观察值
set seed 123456 
gen e = invnorm(uniform())   //干扰项，服从正态分布 e~N(0,1)
gen x = 3 * invnorm(uniform()) //x~N(0, 3^2)
gen t = _n    //[1,100]
gen y = .     //空
tsset t       //
replace y = 1 + 2*x + e
replace y = 1 + -2*x + e if t>50
save xtthres_siml.data, replace


*-基本统计分析
scatter y x   //绘制散点图
reg y x       //回归
est store full  //存储回归结果
reg y x if t<=50
est store left
reg y x if t>50
est store right


*-虚拟变量回归：Chow Test的基本思想
*-如果知道断点 在t=50
gen d = (t>50)  //t=50
gen xd = x*d
reg y x xd
est store strue50
*-如果不知道断点
gen d10 = (t>10)  //假设t=10
gen xd = x*d10
reg y x xd10
est store strue10
gen d20 = (t>20)  //假设t=30
gen xd = x*d20
reg y x xd20
est store strue20
gen d30 = (t>30)  //假设t=30
gen xd = x*d30
reg y x xd30
est store strue30
gen d40 = (t>40)  //假设t=40
gen xd = x*d40
reg y x xd40
est store strue40

local m "full strue10 strue20 strue30 strue40 strue50"
est_table 'm', mtitle('m')


*-网格搜索结构突变点
gen rss = .  //残差列表
forvalues i=5/95{ //i从5到95，不能从1到100，否则左侧或右侧没有数据
    gen di = (t>`i')
    gen xdi = x*di
    qui reg y x xdi
    qui replace rss = e(rss) in `i'  //替换rss中第i个值
    drop di xdi
}
sort rss
list rss t in 1/5  //输出最小的前5个结果

*-图示
line rss t, xline(50, lp(dash) lw(thick)) ///
                  lw(thick) xlabel(0(10)100) sort