# [R] Cochrane-Orcutt method

John Fox jfox at mcmaster.ca
Thu Jan 10 00:30:27 CET 2002

Dear John,

At 03:49 PM 1/9/2002 -0400, you wrote:

>Is there a package that implements the Cochrane-Orcutt itterative
>procedure for dealing with autocorrelation in a regression model?

I'm not aware of an implementation in a package, but I have two-step
Cochrane-Orcott and Prais-Winsten estimators for my own use, which I've
attached to this note.

That said, I'm not sure why you'd want to use these methods, other than for
historical reasons. The gls function in the nlme package will fit linear
models with autocorrelated errors by ML.

I hope that this is of use to you,
John

# these functions take a linear-model object as input

cochrane.orcutt <- function(mod, ...){
UseMethod("cochrane.orcott")
}

cochrane.orcutt.lm <- function(mod){
X <- model.matrix(mod)
y <- model.response(model.frame(mod))
e <- residuals(mod)
n <- length(e)
names <- colnames(X)
rho <- sum(e[1:(n-1)]*e[2:n])/sum(e^2)
y <- y[2:n] - rho * y[1:(n-1)]
X <- X[2:n,] - rho * X[1:(n-1),]
mod <- lm(y ~ X - 1)
result <- list()
result\$coefficients <- coef(mod)
names(result\$coefficients) <- names
summary <- summary(mod, corr = F)
result\$cov <- (summary\$sigma^2) * summary\$cov.unscaled
dimnames(result\$cov) <- list(names, names)
result\$sigma <- summary\$sigma
result\$rho <- rho
class(result) <- 'cochrane.orcutt'
result
}

prais.winsten <- function(mod, ...){
UseMethod("prais.winsten")
}

prais.winsten.lm <- function(mod){
X <- model.matrix(mod)
y <- model.response(model.frame(mod))
e <- residuals(mod)
n <- length(e)
names <- colnames(X)
rho <- sum(e[1:(n-1)]*e[2:n])/sum(e^2)
y <- c(y[1] * (1 - rho^2)^0.5, y[2:n] - rho * y[1:(n-1)])
X <- rbind(X[1,] * (1 - rho^2)^0.5, X[2:n,] - rho * X[1:(n-1),])
mod <- lm(y ~ X - 1)
result <- list()
result\$coefficients <- coef(mod)
names(result\$coefficients) <- names
summary <- summary(mod, corr = F)
result\$cov <- (summary\$sigma^2) * summary\$cov.unscaled
dimnames(result\$cov) <- list(names, names)
result\$sigma <- summary\$sigma
result\$rho <- rho
class(result) <- 'prais.winsten'
result
}

-----------------------------------------------------
John Fox
Department of Sociology
McMaster University
email: jfox at mcmaster.ca
phone: 905-525-9140x23604
web: www.socsci.mcmaster.ca/jfox
-----------------------------------------------------

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._