Orcutt 包使用 OC 程序最后一步的残差计算 DW
The Orcutt Package Calculates DW using residuals from the Last Step of OC Procedure
This began as a question, but after reading the reference provided in the answer below, as well as the source code, the solution became clear. In case anyone else finds themselves in this position:
orcutt 包,版本 2.2 使用特殊程序计算其 CO 模型的 DW 统计数据。 MWE 使用 orcutt 包的例子来表明它的 Durbin-Watson 统计不是基于 OC 估计的残差。
library(orcutt)
data(icecream, package="orcutt")
dw_calc = function(x){sum((x[2:length(x)] - x[1:(length(x)-1)]) ^ 2) / sum(x ^ 2)}
lm = lm(cons ~ price + income + temp, data=icecream)
e = lm$residuals
dw_calc(e)
# 1.02117 <- Durbin-Watson Statistic
coch = cochrane.orcutt(lm)
e = coch$residuals
dw_calc(e)
# 1.006431 <- Durbin-Watson Statistic
coch
# Durbin-Watson statistic
# (original): 1.02117 , p-value: 3.024e-04
# (transformed): 1.54884 , p-value: 5.061e-02
orcutt 包报告 1.54884 但新残差的实际 DW 为 1.006431。报告的值 1.54884 来自收敛过程的最后一轮(参见 Hildreth-Lu)。详见下文:
lm = lm(cons ~ price + income + temp, data=icecream)
reg = lm(cons ~ price + income + temp, data=icecream)
convergence = 8
X <- model.matrix(reg)
Y <- model.response(model.frame(reg))
n<-length(Y)
e<-reg$residuals
e2<-e[-1]
e3<-e[-n]
regP<-lm(e2~e3-1)
rho<-summary(regP)$coeff[1]
rho2<-c(rho)
XB<-X[-1,]-rho*X[-n,]
YB<-Y[-1]-rho*Y[-n]
regCO<-lm(YB~XB-1)
ypCO<-regCO$coeff[1]+as.matrix(X[,-1])%*%regCO$coeff[-1]
e1<-ypCO-Y
e2<-e1[-1]
e3<-e1[-n]
regP<-lm(e2~e3-1)
rho<-summary(regP)$coeff[1]
rho2[2]<-rho
i<-2
while (round(rho2[i-1],convergence)!=round(rho2[i],convergence)){
XB<-X[-1,]-rho*X[-n,]
YB<-Y[-1]-rho*Y[-n]
regCO<-lm(YB~XB-1)
ypCO<-regCO$coeff[1]+as.matrix(X[,-1])%*%regCO$coeff[-1]
e1<-ypCO-Y
e2<-e1[-1]
e3<-e1[-n]
regP<-lm(e2~e3-1)
rho<-summary(regP)$coeff[1]
i<-i+1
rho2[i]<-rho
}
regCO$number.interaction<-i-1
regCO$rho <- rho2[i-1]
regCO$DW <- c(lmtest::dwtest(reg)$statistic, lmtest::dwtest(reg)$p.value,
lmtest::dwtest(regCO)$statistic, lmtest::dwtest(regCO)$p.value)
dw_calc(regCO$residuals)
regF<-lm(YB ~ 1)
tF <- anova(regCO,regF)
regCO$Fs <- c(tF$F[2],tF$`Pr(>F)`[2])
# fitted.value
regCO$fitted.values <- model.matrix(reg) %*% (as.matrix(regCO$coeff))
# coeff
names(regCO$coefficients) <- colnames(X)
# st.err
regCO$std.error <- summary(regCO)$coeff[,2]
# t value
regCO$t.value <- summary(regCO)$coeff[,3]
# p value
regCO$p.value <- summary(regCO)$coeff[,4]
class(regCO) <- "orcutt"
# formula
regCO$call <- reg$call
# F statistics and p value
df1 <- dim(model.frame(reg))[2] - 1
df2 <- length(regCO$residuals) - df1 - 1
RSS <- sum((regCO$residuals)^2)
TSS <- sum((regCO$model[1] - mean(regCO$model[,1]))^2)
regCO$rse <- sqrt(RSS/df2)
regCO$r.squared <- 1 - (RSS/TSS)
regCO$adj.r.squared <- 1 - ((RSS/df2)/(TSS/(df1 + df2)))
regCO$gdl <- c(df1, df2)
#
regCO$rank <- df1
regCO$df.residual <- df2
regCO$assign <- regCO$assign[-(df1+1)]
regCO$residuals <- Y - regCO$fitted.values
regCO
}
我不太了解这个领域,但似乎更有可能有多种 definitions/estimation 方法用于此统计信息。回到 ?orcutt
中引用的书(Verbeek M. (2004) 现代计量经济学指南,John Wiley & Sons Ltd,ISBN:978-88-08 -17054-5),并搜索 Google books 得到
上面示例中 orcutt
计算出的值与书中给出的值一致。在本书的前面(第 108 页)它说
In [the Cochrane-Orcutt procedure], $\rho$ and $\beta$ are recursively estimated until convergence, i.e. having estimated $\beta$ by EGLS (by $\beta^*$), the residuals are recomputed and $\rho$ is estimated again using the residuals from the EGLS step. With this new estimate of $\rho$, EGLS is applied again and one obtains a new estimate of $\beta$ ...
换句话说,您上面给出的 $\rho$ 估计值似乎只对应于 Orcutt-Cochrane 程序的第一步。
This began as a question, but after reading the reference provided in the answer below, as well as the source code, the solution became clear. In case anyone else finds themselves in this position:
orcutt 包,版本 2.2 使用特殊程序计算其 CO 模型的 DW 统计数据。 MWE 使用 orcutt 包的例子来表明它的 Durbin-Watson 统计不是基于 OC 估计的残差。
library(orcutt)
data(icecream, package="orcutt")
dw_calc = function(x){sum((x[2:length(x)] - x[1:(length(x)-1)]) ^ 2) / sum(x ^ 2)}
lm = lm(cons ~ price + income + temp, data=icecream)
e = lm$residuals
dw_calc(e)
# 1.02117 <- Durbin-Watson Statistic
coch = cochrane.orcutt(lm)
e = coch$residuals
dw_calc(e)
# 1.006431 <- Durbin-Watson Statistic
coch
# Durbin-Watson statistic
# (original): 1.02117 , p-value: 3.024e-04
# (transformed): 1.54884 , p-value: 5.061e-02
orcutt 包报告 1.54884 但新残差的实际 DW 为 1.006431。报告的值 1.54884 来自收敛过程的最后一轮(参见 Hildreth-Lu)。详见下文:
lm = lm(cons ~ price + income + temp, data=icecream)
reg = lm(cons ~ price + income + temp, data=icecream)
convergence = 8
X <- model.matrix(reg)
Y <- model.response(model.frame(reg))
n<-length(Y)
e<-reg$residuals
e2<-e[-1]
e3<-e[-n]
regP<-lm(e2~e3-1)
rho<-summary(regP)$coeff[1]
rho2<-c(rho)
XB<-X[-1,]-rho*X[-n,]
YB<-Y[-1]-rho*Y[-n]
regCO<-lm(YB~XB-1)
ypCO<-regCO$coeff[1]+as.matrix(X[,-1])%*%regCO$coeff[-1]
e1<-ypCO-Y
e2<-e1[-1]
e3<-e1[-n]
regP<-lm(e2~e3-1)
rho<-summary(regP)$coeff[1]
rho2[2]<-rho
i<-2
while (round(rho2[i-1],convergence)!=round(rho2[i],convergence)){
XB<-X[-1,]-rho*X[-n,]
YB<-Y[-1]-rho*Y[-n]
regCO<-lm(YB~XB-1)
ypCO<-regCO$coeff[1]+as.matrix(X[,-1])%*%regCO$coeff[-1]
e1<-ypCO-Y
e2<-e1[-1]
e3<-e1[-n]
regP<-lm(e2~e3-1)
rho<-summary(regP)$coeff[1]
i<-i+1
rho2[i]<-rho
}
regCO$number.interaction<-i-1
regCO$rho <- rho2[i-1]
regCO$DW <- c(lmtest::dwtest(reg)$statistic, lmtest::dwtest(reg)$p.value,
lmtest::dwtest(regCO)$statistic, lmtest::dwtest(regCO)$p.value)
dw_calc(regCO$residuals)
regF<-lm(YB ~ 1)
tF <- anova(regCO,regF)
regCO$Fs <- c(tF$F[2],tF$`Pr(>F)`[2])
# fitted.value
regCO$fitted.values <- model.matrix(reg) %*% (as.matrix(regCO$coeff))
# coeff
names(regCO$coefficients) <- colnames(X)
# st.err
regCO$std.error <- summary(regCO)$coeff[,2]
# t value
regCO$t.value <- summary(regCO)$coeff[,3]
# p value
regCO$p.value <- summary(regCO)$coeff[,4]
class(regCO) <- "orcutt"
# formula
regCO$call <- reg$call
# F statistics and p value
df1 <- dim(model.frame(reg))[2] - 1
df2 <- length(regCO$residuals) - df1 - 1
RSS <- sum((regCO$residuals)^2)
TSS <- sum((regCO$model[1] - mean(regCO$model[,1]))^2)
regCO$rse <- sqrt(RSS/df2)
regCO$r.squared <- 1 - (RSS/TSS)
regCO$adj.r.squared <- 1 - ((RSS/df2)/(TSS/(df1 + df2)))
regCO$gdl <- c(df1, df2)
#
regCO$rank <- df1
regCO$df.residual <- df2
regCO$assign <- regCO$assign[-(df1+1)]
regCO$residuals <- Y - regCO$fitted.values
regCO
}
我不太了解这个领域,但似乎更有可能有多种 definitions/estimation 方法用于此统计信息。回到 ?orcutt
中引用的书(Verbeek M. (2004) 现代计量经济学指南,John Wiley & Sons Ltd,ISBN:978-88-08 -17054-5),并搜索 Google books 得到
上面示例中 orcutt
计算出的值与书中给出的值一致。在本书的前面(第 108 页)它说
In [the Cochrane-Orcutt procedure], $\rho$ and $\beta$ are recursively estimated until convergence, i.e. having estimated $\beta$ by EGLS (by $\beta^*$), the residuals are recomputed and $\rho$ is estimated again using the residuals from the EGLS step. With this new estimate of $\rho$, EGLS is applied again and one obtains a new estimate of $\beta$ ...
换句话说,您上面给出的 $\rho$ 估计值似乎只对应于 Orcutt-Cochrane 程序的第一步。