Search code examples
rstatisticstime-serieseconomics

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:

The orcutt package, version 2.2 uses a special procedure to calculate the DW statistic for its CO models. The MWE uses the orcutt package's example to show that its Durbin-Watson statistic is not based on the residuals of the OC estimation.

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

The orcutt package reports 1.54884 but the actual DW is 1.006431 for the new residuals. The reported value, 1.54884, comes from the last round of the convergence procedure (see Hildreth-Lu). See below for a thorough explanation:

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
}
Solution

  • I don't know this area well, but it seems much more likely that there are multiple definitions/estimation methods for this statistic. Going back to the book cited in ?orcutt (Verbeek M. (2004) A guide to modern econometrics, John Wiley & Sons Ltd, ISBN:978-88-08-17054-5), and searching on Google books gives

    enter image description here

    The value shown as computed by orcutt in your example above agrees with the value given in the book. Earlier in the book (p. 108) it says

    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$ ...

    In other words, it seems as though the estimate of $\rho$ that you give above corresponds only to the first step of the Orcutt-Cochrane procedure.