[R] Cochrane-Orcutt method

[John Janmaat]

John Fox wrote:
> 
> 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'm teaching an undergraduate econometrics class.  ML is a little more
than they are up to at this level, and the Cochrane-Orcott approach is
an easily understandable itterative procedure.  Thus, the reason is
historical.  I want a canned procedure to be available to verify results
of the student calculations, rather than running the risk that I do it
incorrectly.

Thanks,

John.

> 
> 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
> Hamilton, Ontario, Canada L8S 4M4
> email: jfox at mcmaster.ca
> phone: 905-525-9140x23604
> web: www.socsci.mcmaster.ca/jfox
> -----------------------------------------------------

-- 
==========================================
John Janmaat
Department of Economics
Acadia University, Wolfville, NS, B0P 1X0
(902)585-1461

All opinions stated are personal, unless
otherwise indicated.
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._