## wrapper to lm() that also has an estimating
## function for the sigma
fw_lm <- function(formula, data = list()) {
  rval <- lm(formula, data = data)
  class(rval) <- c("fw_lm", "lm")
  return(rval)
}
estfun.fw_lm <- function(obj) {
  res <- residuals(obj)
  ef <- estfun.lm(obj)
  sigma2 <- mean(res^2)
  rval <- cbind(ef, res^2 - sigma2)
  colnames(rval) <- c(colnames(ef), "(Variance)")
  return(rval)
}

## simple monitoring function...this could be more elaborate
## but should be sufficient for our current purposes
## Note: critical value is computed from Bonferroni correction for alpha = 0.05 and T = 3
## based on Table III in Zeileis et al. (2005, JAE) and then interpolation via
## round(approxfun(c(0.01, 0.005), c(2.3, 2.478))(0.008513), digits = 3)
fw_monitor <- function(formula, data, monitor, start = NULL, end = NULL, critval = 2.353, hac = TRUE)
{
  if(!(is.null(start) & is.null(end))) data <- window(data, start = start, end = end)

  ## historical data
  monitor <- end(window(data, end = monitor))
  n <- NROW(window(data, end = monitor))
  hist.fm <- fw_lm(formula, data = as.data.frame(window(data, end = monitor)))
  sigma2 <- sum(residuals(hist.fm)^2)/n

  ## full data
  mf <- model.frame(formula, data = as.data.frame(data))
  y <- model.response(mf)
  X <- model.matrix(formula, data = as.data.frame(data))
  attr(X, "assign") <- NULL

  ## estimating functions
  y.pred <- predict(hist.fm, as.data.frame(X))
  mscore <- as.vector(y - y.pred) * X
  mscore <- cbind(mscore, ((y - y.pred)^2 - sigma2))

  ## cumulative estimating functions
  mscore <- mscore/sqrt(n)
  hist.score <- mscore[1:n,]
  J12 <- if(hac) kernHAC(hist.fm, sandwich = FALSE) else crossprod(hist.score)
  J12 <- root.matrix(J12)

  mscore <- apply(as.matrix(mscore), 2, cumsum)
  mscore <- zoo(t(chol2inv(chol(J12)) %*% t(mscore)), order.by = index(data))
  colnames(mscore) <- c(colnames(X), "(Variance)")
  
  rval <- list(process = mscore,
               n = n,
	       formula = formula,
	       data = data,
	       monitor = monitor,
	       critval = critval,
	       J12 = J12)
  class(rval) <- "fw_monitor"
  return(rval)
}

plot.fw_monitor <- function(x,
  ylim = NULL, xlab = "Time", ylab = "Empirical fluctuation process",
  main = "Monitoring of Frankel-Wei regression", ...)
{
  n <- x$n
  proc <- x$process
  maxabs <- zoo(apply(abs(coredata(proc)), 1, max), order.by = index(proc))
  dummy <- zoo((n+1):NROW(proc)/n, index(proc)[-(1:n)]) * x$critval

  if(is.null(ylim)) ylim <- c(0, max(c(max(maxabs), max(dummy))))
  plot(maxabs, ylim = ylim, xlab = xlab, ylab = ylab, main = main, ...)
  abline(h = 0)
  abline(v = x$monitor, lty = 2)
  lines(dummy, col = 2)
  invisible(x)
}

print.fw_monitor <- function(x, ...)
{
  n <- x$n
  proc <- x$process
  maxabs <- apply(abs(coredata(proc)), 1, max)[-(1:n)]
  dummy <- ((n+1):NROW(proc)/n) * x$critval
  breakpoint <- if(any(maxabs > dummy)) min(which(maxabs > dummy)) + n else NA
  form <- as.character(x$formula)
  form <- paste(form[2], form[1], form[3])

  cat("Monitoring of Frankel-Wei regression\n")
  cat("--------------------------------------------------------\n")
  cat(paste("Formula:", form, "\n"))
  cat(paste("History period:", start(proc), "to", x$monitor, "\n"))
  cat(paste("Break detected:", ifelse(is.na(breakpoint), "none", as.character(index(proc)[breakpoint])), "\n"))
  cat("--------------------------------------------------------\n")
  
  invisible(x)
}

breakpoints.fw_monitor <- function(obj, ...)
{
  n <- obj$n
  proc <- obj$process
  maxabs <- apply(abs(coredata(proc)), 1, max)[-(1:n)]
  dummy <- ((n+1):NROW(proc)/n) * obj$critval
  breakpoint <- if(any(maxabs > dummy)) min(which(maxabs > dummy)) + n else NA
  return(breakpoint)
}