准差分法/Cochrane-Orcutt估计法的python实现问题

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

准差分法/Cochrane-Orcutt估计法的python实现问题

1.背景

最近在学习陈强教授的《计量经济学及Stata应用》(第二版)，因为某些原因需要用python来实现计量过程。在学到准差分法（CO估计法）时，没能在statsmodel中找到对应的函数，自己写了一段感觉不太对，请大神指点迷津。

2.内容

教材中涉及的内容直接po出来，免得又去翻书：

准差分法(quasi difference) / Cochrane-Orcutt估计法 （159页）

变换原模型使转换后的扰动项变成球形扰动项。

• 假设原模型：$y_t=\beta_1+\beta_2x_{t2}+\cdots+\beta_Kx_{tK}+\epsilon_t \quad (t=1,\cdots,n) \quad {(8.11)}$
  • 其中 $\epsilon_t$ 存在自相关，且一阶自相关：$\epsilon_t=\rho\epsilon_{t-1}+\mu_t$
    • 自回归系数$|\rho|<1$
    • $\mu_t$为白噪声
• 将原模型滞后一期，两边同乘 $\rho$ ：$\rho y_{t-1}=\rho(\beta_1+\beta_2x_{t-1,2}+\cdots+\beta_Kx_{t-1,K}+\epsilon_{t-1})$
• 方程组（n-1个) 的扰动项为球形扰动项，可消除异方差，且消除了自相关。$y_t-\rho y_{t-1}=(1-\rho)\beta_1+\beta_2(x_{t2}-\rho x_{t-1,2})+\cdots+\beta_K(x_{tK}-\rho x_{t-1,K})+(\underbrace{\epsilon_t-\rho\epsilon_{t-1}}_{\mu_t})\quad {(8.14)}$
  ……

常使用迭代法进行估计，具体步骤：

• 第1步，用OLS估计原模型，用残差 {$e$} 作辅助回归，得到 $\hat\rho^{(1)}$，再用 $\hat\rho^{(1)}$进行CO或PW估计
• 第2步，用CO或PW得到的新残差估计 $\hat\rho^{(2)}$，再用 $\hat\rho^{(2)}$ 进行CO或PW估计
• 依次类推，直至收敛（即相邻两轮的$\rho$与系数估计值之差足够小）。

Stata代码（168页）

.prais consumption temp price income,corc

python代码

# 使用CO方法 需要多次迭代
import statsmodels.api as sm
import pandas as pd
from statsmodels.stats.diagnostic import acorr_breusch_godfrey

# 读取原始数据
icecream = pd.read_stata('../2_Data/Data-2e/icecream.dta')
X = icecream[['temp','price','income']]
y = icecream['consumption']
X = sm.add_constant(X)

# 初始OLS回归
model = sm.OLS(y,X)
results = model.fit() 

# 计算误差项的各期滞后值：datafr ame格式
resid = pd.Datafr ame()
for i in  range(len(icecream)):
	lag_n =  f'lag_{i}'
	resid[lag_n] = results.resid.shift(i)

# Cochrane-Orcutt迭代过程
converged =  False
iterations =  0
max_iterations =  10  # 设置最大迭代次数
tolerance =  0.001  # 设置收敛容差
hat_rho = [0]
results_list = [ ]

while  not converged and iterations < max_iterations:

	# 拟合AR(1)模型来估计自相关系数
	# 残差数据处理
	lag_X =  'lag_{}'.format(iterations+1)
	lag_y =  'lag_{}'.format(iterations)
	resid_ols = resid.copy()[[lag_y, lag_X]]
	resid_ols.dropna(inplace=True)

	# 进行AR(1)模型拟合
	resid_X = resid_ols[lag_X]
	resid_y = resid_ols[lag_y]

	# X = sm.add_constant(X) # 不加入常数项效果更好
	rho = sm.OLS(resid_y, resid_X).fit()
	hat_rho.append(rho.params.iloc[0])

	# 使用Cochrane-Orcutt变换调整误差项
	X_adj = X.copy() - X.copy().shift(1)*hat_rho[-1]
	y_adj = y.copy() - y.copy().shift(1)*hat_rho[-1]

	# 重新进行OLS回归，去除t0项
	results_new = sm.OLS(y_adj[1:], sm.add_constant(X_adj[1:])).fit()

	# 检查是否收敛
	if  abs(hat_rho[-1]-hat_rho[-2]) < tolerance:
		converged =  True
	else:
		results_list.append(results_new)
		results = results_new # 更新模型为新迭代的结果
		iterations +=  1
		print(f"迭代 {iterations}: rho = {hat_rho[-1]:.6f}, Breusch-Godfrey统计量p值{acorr_breusch_godfrey(results,1)[1]:.6f}:")
print('【Cochrane-Orcutt】迭代结果------------------------')
print(f'迭代次数: {iterations}次，模型呈现{iterations-1}阶自相关。')
print('--------------------------------------------【END】')
print(results_list[-2].summary())

---

【错误】：CO估计法都是只滞后一期，上面代码实现的是不停后移滞后项，这肯定不对。

---

修改后的代码：

import statsmodels.api as sm
import pandas as pd

# 读取原始数据,进行初始OLS回归
icecream = pd.read_stata('../2_Data/Data-2e/icecream.dta')
X = icecream[['temp','price','income']]
y = icecream['consumption']
X = sm.add_constant(X)
model = sm.OLS(y,X)
results = model.fit()

# Cochrane-Orcutt迭代过程
converged =  False
iterations =  0
max_iterations =  10  # 设置最大迭代次数
tolerance =  0.001  # 设置收敛容差
hat_rho = [0]
results_list = [ ]

# 初始化第一轮的残差滞后数据
resid = pd.Datafr ame()
resid['lag_0'] = results.resid
resid['lag_1'] = results.resid.shift(1)
resid['mu'] =  0
# print(resid.head())

while not converged and iterations < max_iterations:
	# 拟合AR(1)模型来估计自相关系数
	resid_ols = resid.copy()
	resid_ols.dropna(inplace=True)

	# 进行AR(1)模型拟合
	resid_X = resid_ols['lag_1']
	resid_y = resid_ols['lag_0']
	rho = sm.OLS(resid_y, resid_X).fit()
	hat_rho.append(rho.params.iloc[0])
	print(f"迭代 {iterations+1}: rho = {hat_rho[-1]:.6f}")

	# 使用Cochrane-Orcutt变换后再进行回归，得到新残差项
	X_adj = X.copy() - X.copy().shift(1)*hat_rho[-1]
	y_adj = y.copy() - y.copy().shift(1)*hat_rho[-1]
	results_new = sm.OLS(y_adj[1:], sm.add_constant(X_adj[1:])).fit()

	# 检查是否收敛
	if  abs(hat_rho[-1]-hat_rho[-2]) < tolerance:
	    converged =  True
	else:
	    resid['mu'] = results_new.resid
	    resid.loc[1:,'lag_0']= hat_rho[-1] * resid.loc[1:,'lag_1'] + resid['mu'][1:]
	    resid['lag_1'] = resid['lag_0'].shift(1)
	    # print(resid.head(3))
	    # 更新模型为新迭代的结果
	    results = results_new
	    results_list.append(results_new)
	    iterations +=  1

print('【Cochrane-Orcutt】迭代结果------------------------')
print(results_list[-1].summary())

在else里面，对原模型8.11中的$\epsilon_t$做了还原，因为如果按上文步骤的第2步，做出来结果不收敛。所以用8.14的残差$\mu_t$和$\hat\rho$进行了如下计算：$\epsilon_t = \hat\rho\epsilon_{t-1}+\mu_t$
这里的$\epsilon_1$直接取的是初始OLS回归后的残差序列第1期的值。但还是算不出来教材里面的结果。
计算结果：

迭代 1: rho = 0.400633
迭代 2: rho = 0.396399
迭代 3: rho = 0.396524
【Cochrane-Orcutt】迭代结果------------------------
                            OLS Regression Results                            
==============================================================================
Dep. Variable:            consumption   R-squared:                       0.651
Model:                            OLS   Adj. R-squared:                  0.609
Method:                 Least Squares   F-statistic:                     15.54
Date:                Tue, 23 Apr 2024   Prob (F-statistic):           6.52e-06
Time:                        15:45:05   Log-Likelihood:                 60.905
No. Observations:                  29   AIC:                            -113.8
Df Residuals:                      25   BIC:                            -108.3
Df Model:                           3                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          0.1547      0.289      0.535      0.597      -0.441       0.750
temp           0.0036      0.001      6.445      0.000       0.002       0.005
price         -0.8905      0.811     -1.098      0.282      -2.560       0.779
income         0.0032      0.002      2.094      0.047    5.36e-05       0.006
==============================================================================
Omnibus:                        3.070   Durbin-Watson:                   1.546
Prob(Omnibus):                  0.216   Jarque-Bera (JB):                1.643
Skew:                           0.418   Prob(JB):                        0.440
Kurtosis:                       3.813   Cond. No.                     8.60e+03
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 8.6e+03. This might indicate that there are
strong multicollinearity or other numerical problems.

数据及代码

如果愿意复现上述过程，数据可以在[陈强教授的计量经济学及Stata主页]下载《计量经济学及Stata应用》(第二版)对应的数据集，然后修改上述代码中的对应数据的路径即可。

问题：

• 迭代收敛：“相邻两轮的$\rho$与系数估计值之差足够小”。
  • 应该是$\rho^{(i)}$和$\rho^{(i+1)}$之间的差，前后两轮系数估计值之间的差吧，不然$\rho$和系数估计值有可比性吗？
  • 收敛到多少算足够小？
• 代码问题：迭代过程到第2轮，发现$\rho$已经与第一轮足够收敛了，那就应该取第一轮的OLS结果，应该是这样理解吗？

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏1 回帖

关键词：Cochrane Cochran python CHR cut

xujingjun 发表于 2024-4-23 12:48

麻烦大佬帮忙看看这个哪里出问题了。可能是对这个迭代过程理解错误，我自己又修改了下，但是无论怎么调参数都算不出和书上一样的结果。

感谢分享。

newfei188 发表于 2024-4-24 10:57

意思是，我这个计算方法是正确的吗？

watalo 发表于 2024-4-22 12:32
准差分法/Cochrane-Orcutt估计法的python实现问题
1.背景
最近在学习陈强教授的《计量经济学及Stata应用》 ...

原本以为一切都无足轻重，认识你后，我终于有了手足的感觉；至交后，我终于知道了什么叫情同手足。都说朋友如手足，确实不假。所以记得天天洗手洗足哟！