时间过的真快，转眼间，又是一年。
    翻了翻论坛，截止到2017年的最后一天，好像一共只发了5篇帖子。本来想在原来帖子的基础上加以完善补充，但是发现时间超限，无法再进行编辑了，只好再发一贴。
    这一贴主要针对以下两篇帖子进行Python版程序的部分补充。
        1、 面板数据非线性回归模型简介（门限回归，马尔可夫体制转换，平滑转换）
        https://bbs.pinggu.org/thread-1144995-1-1.html        2、面板数据非线性回归模型建模方法及其应用
        https://bbs.pinggu.org/thread-2574681-1-1.html


    最近开始学习Python，于是想着用Python把以前的程序都改写一下。鉴于工作后的这六年时间里，已经换了6台电脑，资料都不知道具体被转移到哪里了，只能重头再来。最关键的是，专业课丢下太久，再看文献觉得甚是生疏。
    为了给即将1岁9个月大的儿子做好榜样，只好硬着头皮，慢慢拾起来。于是，有了下面的第一篇改编。毕竟刚接触Python不久，程序里肯定还存在一些尙待优化的地方，不足之处还望各位读者多多包涵。
    以下为基于Python3.6改写的Hansen（1999）静态面板数据门限回归模型的估计及检验程序，结构略有调整，经验证结果完全一致。程序分为两大部分，第一部分为各种自定义函数；第二部分为程序的数据导入及参数设置。下面对第二部分进行简要的介绍，以方便读者了解程序结构，并根据需要选择下载。

#=======================================================================
#============PART2: 数据导入及模型估计检验区域，请根据实际需求修改  ======================
#=======================================================================

# step1 ：导入数据
# 注意事项：1、数据需与程序在同一文件夹内；
#                 2、数据的组织形式为：每一列为一个变量，每一列先"N(个体)"后"T(时期)"，即
    #         Y           X1  ...         Xk
    #    ---------------------------------------------------------
    #    |Y(N1T1)  X1(N1T1) ... Xk(N1T1)|
    #    |Y(N1T2)  X1(N1T2) ... Xk(N1T2)|
    #    |  ...     ...     ...  ...                    |     个体1
    #    |Y(N1Tt)  X1(N1Tt) ...  Xk(N1Tt) |
    #    ----------------------------------------------------------
    #    |Y(N2T1)  X1(N2T1) ... Xk(N2T1)|
    #    |Y(N2T2)  X1(N2T2) ... Xk(N2T2)|
    #    |...  ...  ...    ...                   ...    |     个体2
    #    |Y(N2Tt)  X1(N2T2) ... Xk(N2Tt) |   
    #    -----------------------------------------------------------

data = np.asmatrix(np.loadtxt('invest.txt'))  
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

# step2 ：设置变量，只需根据txt中的数据排列顺序，更改所需的列数。需注意第一列从【 0 】开始
i = data[:,0]  # 设定被解释变量Y
q = data[:,1]  # 设定不受门限效应影响的解释变量 X_no
c = data[:,2]  # 设定受门限效应影响的解释变量 X_yes
d = data[:,3]  # 设定门限变量 q
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


# step3 ：设置参数
t=15                 #  时期数（需设置）
nt = len(data[:,0])  # （自动计算，勿改）
n = int(nt/t)        # （自动计算，勿改）
qn = 400             # 门限变量的分位数数量（需设置）
conf_lev = .95       # 置信区间（需设置）
boot_1_ = 10          # 单门限模型自举次数（需设置）
boot_2_ = 10          # 双门限模型自举次数（需设置）
boot_3_ = 10          # 三门限模型自举次数（需设置）
# 如需4门限，或更多门限，只需新增变量即可，如boot_4_ = 300，也可在函数中直接输入
trim_1_ = .01        # 单门限模型门限变量剔除比例（需设置）
trim_2_ = .01        # 双门限模型门限变量剔除比例（需设置）
trim_3_ = .05        # 三门限模型门限变量剔除比例（需设置）
max_lag = 1          # 解释变量中，最大的滞后阶数（需设置）
tt = t-max_lag       # （自动计算，勿改）
ty = n*(t-max_lag-1) # （自动计算，勿改）
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


# step4 ：变量的滞后及去均值处理
y  = lag_v(i,0)      # 被解释变量的滞后处理（滞后阶数为0时，实际上未滞后，只是为了对齐数据，剔除了一部分数据）
yt = tr(y)           # 被解释变量的去均值处理
cf = lag_v(c,1)      # 受门限效应影响的解释变量的滞后处理
ct = tr(cf)          # 受门限效应影响的解释变量的去均值处理
q1 = lag_v(q,1)      # 不受门限效应影响的解释变量的滞后处理
d1 = lag_v(d,1)      # 门限变量的滞后处理
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~



# step5 ：所有解释变量的横向组合
# 如模型中有多个不受门限效应影响的解释变量，则可以根据data的组织形式，加以填充。
# 如data中的第2、3、5列均为不受门限效应影响的解释变量，以下表达式可改为：
# x=np.hstack([np.array(data[1,:]),np.array(data[2,:]),np.array(data[4,:])])
# 还可根据需要进行平方、立方、或者乘积运算。

# 下行为设置不受门限效应影响的解释变量，可以根据需要改动，参考以上说明进行改动。
x=np.hstack([np.array(q1),np.array(q1)**2,np.array(q1)**3,d1,np.multiply(q1,d1)])   
# 以下命令行均为自动计算，无需修改。
k=x.shape[1]
xt=np.asmatrix(np.zeros([n*(tt-1),k]))
j=0
while j<=k-1:
    xt[:,j]=tr(x[:,j])
    j=j+1                 # 所有不受门限效应影响的解释变量的去均值处理

xxt = np.hstack([xt,ct])  # 将不受门限效应影响及受门限效应影响的解释变量，通通去均值处理后，进行横向的合并拼接
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


# step6 ：根据step3设置的分位数量qn及剔除比例trim_，进行门限变量的筛选。
# 以下均为自动计算，无需修改。
thresh=d1
dd=np.asmatrix(np.unique(np.asarray(thresh))).T # 门限变量去重并升序排列
qnt1=qn*trim_1_
sq=np.arange(trim_1_,(trim_1_+(1/qn)*(qn-2*qnt1+1)),(1/qn))
qq1=dd[np.floor(sq*dd.size-1).astype(int).tolist()]
qn1=qq1.size
cc=-2*np.log(1-np.sqrt(conf_lev))

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


# step7 ：模型估计及检验
# 说明：model(门限值【列矩阵】、剔除比例【数值】、自举次数【数值】)
# 0、无门限效应下，模型的残差平方和计算   
sse0 = sse_calc(yt,xxt,1)

# 1、单门限效应模型
rhat1,min_cfs1,max_cfs1,beta1,sehet1,sehomo1,lrt1,crits1,p_value1,rrr1,sse1= model(0,trim_1_,boot_1_)

# 2.1、双门限效应模型
rhat2,min_cfs2,max_cfs2,beta2,sehet2,sehomo2,lrt2,crits2,p_value2,rrr21,sse2= model(rhat1,trim_2_,boot_2_)

# 2.2、双门限效应模型的二次优化
rhat1,min_cfs1,max_cfs1,beta1,sehet1,sehomo1,lrt1,crits1,p_value1,rrr22,sse2= model(rhat2,trim_2_,boot_2_)

# 3、三门限效应模型
rhat3,min_cfs3,max_cfs3,beta3,sehet3,sehomo3,lrt3,crits3,p_value3,rrr3,sse3= model(np.vstack([rhat1,rhat2]),trim_3_,boot_3_)

# 4、四门限效应模型
rhat4,min_cfs4,max_cfs4,beta4,sehet4,sehomo4,lrt4,crits4,p_value4,rrr4,sse4= model(np.vstack([rhat1,rhat2,rhat3]),trim_3_,1)


    为了节省时间，自举次数设置的比较少，各位同学在使用过程中，可以根据实际情况自行设置。设置完参数以后，点击运行，结果及图形显示示例如下：









    代码只为学习Python使用，熟悉MATLAB或者R的可以去Hansen的个人网页直接下载数据和程序。
至于后续会不会有Python代码之二，就看时间允不允许吧~