1bgtest <- function(formula, order = 1, order.by = NULL, type = c("Chisq", "F"), data = list(), fill = 0)
2{
3  dname <- paste(deparse(substitute(formula)))
4
5  if(!inherits(formula, "formula")) {
6    X <- if(is.matrix(formula$x))
7           formula$x
8         else model.matrix(terms(formula), model.frame(formula))
9    y <- if(is.vector(formula$y))
10           formula$y
11         else model.response(model.frame(formula))
12  } else {
13    mf <- model.frame(formula, data = data)
14    y <- model.response(mf)
15    X <- model.matrix(formula, data = data)
16  }
17
18  if(!is.null(order.by))
19  {
20    if(inherits(order.by, "formula")) {
21      z <- model.matrix(order.by, data = data)
22      z <- as.vector(z[,ncol(z)])
23    } else {
24      z <- order.by
25    }
26    X <- as.matrix(X[order(z),])
27    y <- y[order(z)]
28  }
29
30  n <- nrow(X)
31  k <- ncol(X)
32  order <- 1:order
33  m <- length(order)
34  resi <- lm.fit(X,y)$residuals
35
36  Z <- sapply(order, function(x) c(rep(fill, length.out = x), resi[1:(n-x)]))
37  if(any(na <- !complete.cases(Z))) {
38    X <- X[!na, , drop = FALSE]
39    Z <- Z[!na, , drop = FALSE]
40    y <- y[!na]
41    resi <- resi[!na]
42    n <- nrow(X)
43  }
44  auxfit <- lm.fit(cbind(X,Z), resi)
45
46  cf <- auxfit$coefficients
47  vc <- chol2inv(auxfit$qr$qr) * sum(auxfit$residuals^2) / auxfit$df.residual
48  names(cf) <- colnames(vc) <- rownames(vc) <- c(colnames(X), paste("lag(resid)", order, sep = "_"))
49
50  switch(match.arg(type),
51
52  "Chisq" = {
53    bg <- n * sum(auxfit$fitted.values^2)/sum(resi^2)
54    p.val <- pchisq(bg, m, lower.tail = FALSE)
55    df <- m
56    names(df) <- "df"
57  },
58
59  "F" = {
60    uresi <- auxfit$residuals
61    bg <- ((sum(resi^2) - sum(uresi^2))/m) / (sum(uresi^2) / (n-k-m))
62    df <- c(m, n-k-m)
63    names(df) <- c("df1", "df2")
64    p.val <- pf(bg, df1 = df[1], df2 = df[2], lower.tail = FALSE)
65  })
66
67  names(bg) <- "LM test"
68  RVAL <- list(statistic = bg, parameter = df,
69               method = paste("Breusch-Godfrey test for serial correlation of order up to", max(order)),
70               p.value = p.val,
71               data.name = dname,
72	       coefficients = cf,
73	       vcov = vc)
74  class(RVAL) <- c("bgtest", "htest")
75  return(RVAL)
76}
77
78vcov.bgtest <- function(object, ...) object$vcov
79df.residual.bgtest <- function(object, ...) if(length(df <- object$parameter) > 1L) df[2] else NULL
80
81