[R] Problems with Cochrane-Orcutt procedures

Tue Jun 17 19:49:22 CEST 2008

Dear Tolga,

That's a little more information, but because the code seems to work for me
on other data (though no longer the message dispatch), I can't say what
produces the error. I guess that if you can't debug this yourself, you'll
have to share the data (generally a good idea in any event).

> mod <- lm(Hartnagel[,5] ~ Hartnagel[,7])
> mod

Call:
lm(formula = Hartnagel[, 5] ~ Hartnagel[, 7])

Coefficients:
   (Intercept)  Hartnagel[, 7]  
      -85.1117          0.2126  

> cochrane.orcutt.lm(mod)
$coefficients
   (Intercept) Hartnagel[, 7] 
   57.05592426     0.04015223 

$cov
               (Intercept) Hartnagel[, 7]
(Intercept)    1226.739528   -1.381235045
Hartnagel[, 7]   -1.381235    0.001713204

$sigma
[1] 14.00277

$rho
[1] 0.7835837

attr(,"class")
[1] "cochrane.orcutt"

Regards,
 John

------------------------------
John Fox, Professor
Department of Sociology
McMaster University
Hamilton, Ontario, Canada
web: socserv.mcmaster.ca/jfox

> -----Original Message-----
> From: tolga.i.uzuner at jpmorgan.com [mailto:tolga.i.uzuner at jpmorgan.com]
> Sent: June-17-08 1:17 PM
> To: John Fox
> Cc: r-help at r-project.org; tolga.i.uzuner at jpmorgan.com
> Subject: RE: [R] Problems with Cochrane-Orcutt procedures
> 
> 
> Sure, I can imagine GLS is a much better way to deal with this.
> 
> I guess I was looking at this because I did try GLS but got exactly the
same
> results as LM and I just wanted to be sure.
> 
> I did "debug" the code in https://stat.ethz.ch/pipermail/r-help/2002-
> January/017774.html and the offending line is:
> 
> ...
> mod <- lm(y ~ X - 1)
> ...
> 
> in the Orcutt procedure. Essentially, X is twice the length of y, as
below:
> 
> > cochrane.orcutt(regrCMSlm)
> debugging in: cochrane.orcutt(regrCMSlm)
> debug: {
>     UseMethod("cochrane.orcutt")
> }
> Browse[1]>
> debug: UseMethod("cochrane.orcutt")
> Browse[1]>
> debugging in: cochrane.orcutt.lm(regrCMSlm)
> debug: {
>     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
> }
> Browse[1]>
> debug: X <- model.matrix(mod)
> Browse[1]>
> debug: y <- model.response(model.frame(mod))
> Browse[1]> length(X)
> [1] 378 #<<<<<<<<<<<<<<<<------------------------------
> Browse[1]>
> debug: e <- residuals(mod)
> Browse[1]> length(y)
> [1] 189 #<<<<<<<<<<<<<<<<------------------------------
> Browse[1]>
> debug: n <- length(e)
> Browse[1]>
> debug: names <- colnames(X)
> Browse[1]>
> debug: rho <- sum(e[1:(n - 1)] * e[2:n])/sum(e^2)
> Browse[1]>
> debug: y <- y[2:n] - rho * y[1:(n - 1)]
> Browse[1]>
> debug: X <- X[2:n, ] - rho * X[1:(n - 1), ]
> Browse[1]>
> debug: mod <- lm(y ~ X - 1)
> Browse[1]>
> Error in model.frame.default(formula = y ~ X - 1, drop.unused.levels =
TRUE)
> :
>   variable lengths differ (found for 'X')
> 
> Tolga
> 
> 
> 
> 
> "John Fox" <jfox at mcmaster.ca>
> 
> 17/06/2008 17:51 To
> <tolga.i.uzuner at jpmorgan.com>
> cc
> <r-help at r-project.org>
> Subject
> RE: [R] Problems with Cochrane-Orcutt procedures
> 
> 
> 
> 
> 
> 
> Dear Tolga,
> 
> I'm afraid that I don't see an error. (I expect in any event that the
> Cochrane-Orcott and Prais estimators are now only of historical interest.)
> 
> Regards,
> John
> 
> ------------------------------
> John Fox, Professor
> Department of Sociology
> McMaster University
> Hamilton, Ontario, Canada
> web: socserv.mcmaster.ca/jfox
> 
> 
> > -----Original Message-----
> > From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
> On
> > Behalf Of tolga.i.uzuner at jpmorgan.com
> > Sent: June-17-08 11:13 AM
> > To: r-help at r-project.org
> > Subject: [R] Problems with Cochrane-Orcutt procedures
> >
> > Hi John,
> >
> > Hi Folks/Prof. Fox,
> >
> > I found some code John Fox had written sometime back on the
> Cochrane-Orcutt
> > and Prais procedures here:
> > https://stat.ethz.ch/pipermail/r-help/2002-January/017774.html
> >
> > I thought I would try it out and get the following errors below. Was
> > wondering if anyone had any immediate opinions why this might be ?
> >
> > The linear model is the object regrCMSlm .
> >
> > Thanks,
> > Tolga
> >
> > > regrCMSlm
> >
> > Call:
> > lm(formula = regrCMS[, 1] ~ regrCMS[, 2])
> >
> > Coefficients:
> >  (Intercept)  regrCMS[, 2]
> >      25.7067       -0.3409
> >
> > > summary(regrCMSlm)
> >
> > Call:
> > lm(formula = regrCMS[, 1] ~ regrCMS[, 2])
> >
> > Residuals:
> > 09/20/07 11/28/07 02/01/08 04/09/08 06/16/08
> >  10.0593   0.3588  -1.1459   0.1340  -9.8520
> >
> > Coefficients:
> >              Estimate Std. Error t value Pr(>|t|)
> > (Intercept)  25.70673    0.85300   30.14   <2e-16 ***
> > regrCMS[, 2] -0.34092    0.02205  -15.46   <2e-16 ***
> > ---
> > Signif. codes:  0 b
> 
> 
> 
> 
> 
> 
> ________________________________
> 
> Generally, this communication is for informational purposes only and it is
> not intended as an offer or solicitation for the purchase or sale of any
> financial instrument or as an official confirmation of any transaction. In
> the event you are receiving the offering materials attached below related
to
> your interest in hedge funds or private equity, this communication may be
> intended as an offer or solicitation for the purchase or sale of such
> fund(s). All market prices, data and other information are not warranted
as
> to completeness or accuracy and are subject to change without notice. Any
> comments or statements made herein do not necessarily reflect those of
> JPMorgan Chase & Co., its subsidiaries and affiliates. This transmission
may
> contain information that is privileged, confidential, legally privileged,
> and/or exempt from disclosure under applicable law. If you are not the
> intended recipient, you are hereby notified that any disclosure, copying,
> distribution, or use of the information contained herein (including any
> reliance thereon) is STRICTLY PROHIBITED. Although this transmission and
any
> attachments are believed to be free of any virus or other defect that
might
> affect any computer system into which it is received and opened, it is the
> responsibility of the recipient to ensure that it is virus free and no
> responsibility is accepted by JPMorgan Chase & Co., its subsidiaries and
> affiliates, as applicable, for any loss or damage arising in any way from
its
> use. If you received this transmission in error, please immediately
contact
> the sender and destroy the material in its entirety, whether in electronic
or
> hard copy format. Thank you. Please refer to
> http://www.jpmorgan.com/pages/disclosures for disclosures relating to UK
> legal entities.