[R] Problems with Cochrane-Orcutt procedures
John Fox
jfox at mcmaster.ca
Tue Jun 17 20:30:18 CEST 2008
Dear Tolga,
I'm afraid that your data work fine for me:
> regrCMSlm <- lm(regrCMS[,1] ~ regrCMS[,2])
> cochrane.orcutt.lm(regrCMSlm)
$coefficients
(Intercept) regrCMS[, 2]
23.5679065 -0.1784187
$cov
(Intercept) regrCMS[, 2]
(Intercept) 60.449366 -2.14491371
regrCMS[, 2] -2.144914 0.07676436
$sigma
[1] 1.081036
$rho
[1] 0.7208652
attr(,"class")
[1] "cochrane.orcutt"
>
I'm using R 2.7.0 on Windows Vista; I put your data in a standard data
frame; and, as I said, the method dispatch no longer works, so I called
Cochrane.orcutt.lm() directly.
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 2:13 PM
> To: John Fox
> Cc: r-help at r-project.org
> Subject: Re: [R] Problems with Cochrane-Orcutt procedures
>
> Sure, of course. And thanks for looking John.
>
> Here is the entire data:
>
> > regrCMS
> a b
> 09/20/07 26.084 28.40
> 09/21/07 22.458 28.90
> 09/24/07 21.297 29.25
> 09/25/07 21.733 29.40
> 09/26/07 21.319 28.75
> 09/27/07 22.507 28.85
> 09/28/07 19.571 28.90
> 10/01/07 21.961 29.00
> 10/02/07 21.729 28.45
> 10/03/07 21.241 28.00
> 10/04/07 21.253 28.55
> 10/05/07 21.401 27.25
> 10/10/07 20.923 25.70
> 10/11/07 20.401 25.25
> 10/12/07 20.122 24.35
> 10/15/07 19.225 26.35
> 10/16/07 18.759 27.15
> 10/17/07 18.413 26.85
> 10/18/07 18.960 27.75
> 10/19/07 18.643 28.65
> 10/22/07 17.829 28.35
> 10/23/07 17.835 28.60
> 10/24/07 17.335 29.15
> 10/25/07 17.971 29.05
> 10/26/07 16.205 28.60
> 10/29/07 16.722 28.40
> 10/30/07 16.772 28.40
> 10/31/07 17.127 27.05
> 11/01/07 15.328 27.85
> 11/02/07 17.088 28.15
> > regrCMSlm<-lm(regrCMS[,1]~regrCMS[,2])
> > summary(regrCMSlm)
>
> Call:
> lm(formula = regrCMS[, 1] ~ regrCMS[, 2])
>
> Residuals:
> 09/20/07 10/01/07 10/12/07 10/23/07 11/02/07
> 6.4365 2.2112 0.3183 -2.2462 -2.5355
>
> Coefficients:
> Estimate Std. Error t value Pr(>|t|)
> (Intercept) 16.92553 10.21538 1.657 0.109
> regrCMS[, 2] 0.09584 0.36477 0.263 0.795
>
> Residual standard error: 2.436 on 28 degrees of freedom
> Multiple R-squared: 0.00246, Adjusted R-squared: -0.03317
> F-statistic: 0.06904 on 1 and 28 DF, p-value: 0.7947
>
> > cochrane.orcutt(regrCMSlm)
> 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]>
> debug: e <- residuals(mod)
> 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')
> >
>
>
>
>
> "John Fox" <jfox at mcmaster.ca>
> 17/06/2008 18:49
>
> To
> <tolga.i.uzuner at jpmorgan.com>
> cc
> <r-help at r-project.org>
> Subject
> RE: [R] Problems with Cochrane-Orcutt procedures
>
>
>
>
>
>
> 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.
>
>
>
>
>
>
> 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.
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
More information about the R-help
mailing list