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: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED]
> Sent: June-17-08 1:17 PM
> To: John Fox
> Cc: r-help@r-project.org; [EMAIL PROTECTED]
> 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" <[EMAIL PROTECTED]>
>
> 17/06/2008 17:51 To
> <[EMAIL PROTECTED]>
> cc
> <r-help@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: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED]
> On
> > Behalf Of [EMAIL PROTECTED]
> > Sent: June-17-08 11:13 AM
> > To: r-help@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.
______________________________________________
R-help@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.
