R/gbreakpoints.R

Defines functions confint.gbreakpoints AIC.gbreakpointsfull AIC.gbreakpoints logLik.gbreakpointsfull logLik.gbreakpoints plot.summary.gbreakpointsfull plot.gbreakpointsfull print.summary.gbreakpointsfull summary.gbreakpointsfull summary.gbreakpoints breakdates.gbreakpoints print.gbreakpoints breakpoints.gbreakpointsfull gbreakpoints nlogLik

Documented in AIC.gbreakpoints AIC.gbreakpointsfull breakdates.gbreakpoints breakpoints.gbreakpointsfull confint.gbreakpoints gbreakpoints logLik.gbreakpoints plot.gbreakpointsfull plot.summary.gbreakpointsfull print.gbreakpoints print.summary.gbreakpointsfull summary.gbreakpoints summary.gbreakpointsfull

nlogLik <- function(object, ...) -logLik(object, ...)

gbreakpoints <- function(formula, data, fit = lm, objfun = nlogLik,
                         h = 0.15, breaks = NULL, order.by = NULL,
			 ic = c("LWZ", "BIC"), hpc = c("none", "foreach"), ...)
{
  n <- NROW(data)
  k <- attr(objfun(fit(formula, data = data, ...)), "df")    
  if(h < 1) h <- floor(n*h)
  h <- round(h, digits = 0)
  if(h > floor(n/2))
    stop("minimum segment size must be smaller than half the number of observations")
  if(is.null(breaks)) breaks <- ceiling(n/h) - 2

  ## use order.by
  if(is.null(order.by)) {
    if(is.ts(data)) order.by <- time(data)
      else if(is.zoo(data)) order.by <- index(data)
      else order.by <- 1:n  
  }
  data <- data[ORDER(order.by), , drop = FALSE]
  order.by <- order.by[ORDER(order.by)]

  hpc <- match.arg(hpc)
  if(hpc == "foreach") {
    if(requireNamespace("foreach")) {
      `%dopar%` <- foreach::`%dopar%`
    } else {
      warning("High perfomance computing (hpc) support with 'foreach' package is not available, foreach is not installed.")    
      hpc <- "none"
    }
  }

  ## compute objfun() for all valid segmentes i, ..., j.
  
  #Z# for the default fit & objfun, this can be speeded up
  #Z# by re-using breakpoints(), otherwise this has to be done
  #Z# in a double for-loop and becomes *really* slow
  if(missing(fit) & missing(objfun)) {    
    #Z# obtain RSS-based breakpoints
    lm_bp <- breakpoints(formula, h = h, data = as.data.frame(data), hpc = hpc)
    
    RSS2obj.triang <- function(x) {
      rval <- x$RSS.triang
      n <- x$nobs
      df <- x$nreg + 1
      h <- x$nobs - length(rval)

      rss2nloglik <- function(i) {
        rss <- rval[[i]]
        ni <- seq_along(rss)
        loglik <- 0.5 * ni * (log(rss) + 1 - log(ni) + log(2 * pi))
        loglik
      }
      rval <- lapply(seq_along(rval), rss2nloglik)
      rval <- lapply(rval, function(x) x[-(1:h)])

      return(rval)
    }
    
    #Z# convert RSS to negative log-likelihood
    obj.triang <- RSS2obj.triang(lm_bp)

  } else {
    #Z# this becomes *really* slow...
    #Z# (but can be alleviated to some degree by hpc option)
    obj.triang <- list()
    if(hpc == "none") {
      for(i in 1:(n-h+1)) {
        obj_j <- rep(0, (n-i-h+2))
        for(j in (i+h-1):n) obj_j[(j-i-h+2)] <- as.numeric(objfun(fit(formula, data = data[(i:j),,drop = FALSE], ...)))
        obj.triang[[i]] <- obj_j
      }
    } else {
      obj.triang <- foreach::foreach(i = 1:(n-h+1)) %dopar% {
        obj_j <- rep(0, (n-i-h+2))
        for(j in (i+h-1):n) obj_j[(j-i-h+2)] <- as.numeric(objfun(fit(formula, data = data[(i:j),,drop = FALSE], ...)))
        obj_j
      }
    }
  }  
  #Z# For other special models (combinations of fit and objfun)
  #Z# it would be great to have external methods that generate
  #Z# something that can be easily transformed in an obj.triang
  #Z# object.

  ## function to extract the objective(i,j) from obj.triang
  objective <- function(i,j) {
    if((j < (i + h - 1)) || i < 1 || j > n) NA
      else obj.triang[[i]][j-i-h+2]
  }

  ## compute optimal previous partner if observation i is the mth break
  ## store results together with objective in obj.table

  ## breaks = 1
  index <- h:(n-h)
  obj.break <- sapply(index, function(i) objective(1,i))
  obj.table <- cbind(index, obj.break)
  rownames(obj.table) <- as.character(index)

  ## breaks >= 2

  extend.obj.table <- function(obj.table, breaks)
  {
    if((breaks*2) > ncol(obj.table)) {
      for(m in (ncol(obj.table)/2 + 1):breaks)
      {
        my.index <- (m*h):(n-h)
        my.obj.table <- obj.table[,c((m-1)*2 - 1, (m-1)*2)]
        my.obj.table <- cbind(my.obj.table, NA, NA)
        for(i in my.index)
        {
          pot.index <- ((m-1)*h):(i - h)
          obj.break <- sapply(pot.index, function(j) my.obj.table[as.character(j), 2] + objective(j+1,i))
          opt <- which.min(obj.break)
          my.obj.table[as.character(i), 3:4] <- c(pot.index[opt], obj.break[opt])
        }
        obj.table <- cbind(obj.table, my.obj.table[,3:4])
      }
      colnames(obj.table) <- as.vector(rbind(paste("break", 1:breaks, sep = ""),
                                             paste("objective", 1:breaks, sep = "")))
    }
    return(obj.table)
  }

  obj.table <- extend.obj.table(obj.table, breaks)

  ## extract optimal breaks

  extract.breaks <- function(obj.table, breaks)
  {
    if((breaks*2) > ncol(obj.table)) stop("compute obj.table with enough breaks before")
    index <- obj.table[, 1, drop = TRUE]
    obj.break <- sapply(index, function(i) obj.table[as.character(i),breaks*2] + objective(i + 1, n))
    opt <- index[which.min(obj.break)]
    if(breaks > 1) {
      for(i in ((breaks:2)*2 - 1))
        opt <- c(obj.table[as.character(opt[1]),i], opt)
    }
    names(opt) <- NULL
    return(opt)
  }

  opt <- extract.breaks(obj.table, breaks)

  RVAL <- list(breakpoints = opt,
               obj.table = obj.table,
	       obj.triang = obj.triang,
	       objective = objective,
	       extract.breaks = extract.breaks,
	       extend.obj.table = extend.obj.table,
	       nobs = n,
	       npar = k,
	       fit = fit,
	       fitname = deparse(substitute(fit)),
	       objfun = objfun,
	       objname = deparse(substitute(objfun)),
	       icname = match.arg(ic),
	       call = match.call(),
	       index = order.by)
  class(RVAL) <- c("gbreakpointsfull", "gbreakpoints", "breakpointsfull", "breakpoints")
  RVAL$breakpoints <- breakpoints(RVAL)$breakpoints
  return(RVAL)
}

breakpoints.gbreakpointsfull <- function(obj, breaks = NULL, ...)
{
  if(is.null(breaks)) breaks <- obj$icname
  if(is.character(breaks))
  {
    breaks <- match.arg(breaks, c("BIC", "LWZ"))
    IC <- summary(obj)$objective[breaks,]
    breaks <- if(any(is.na(IC))) length(IC)-1 else which.min(IC)-1
  }
  if(breaks < 1)
  {
    breakpoints <- NA
    objective <- obj$objective(1, obj$nobs)
  } else {
    obj.tab <- obj$extend.obj.table(obj$obj.table, breaks)
    breakpoints <- obj$extract.breaks(obj.tab, breaks)
    bp <- c(0, breakpoints, obj$nobs)
    objective <- sum(apply(cbind(bp[-length(bp)]+1,bp[-1]), 1,
                     function(x) obj$objective(x[1], x[2])))
  }
  RVAL <- list(breakpoints = breakpoints,
               objective = objective,
               nobs = obj$nobs,
	       npar = obj$npar,
	       call = match.call(),
	       fitname = obj$fitname,
	       objname = obj$objname,
               index = obj$index)
  class(RVAL) <- c("gbreakpoints", "breakpoints")
  return(RVAL)
}


print.gbreakpoints <- function(x, format.times = NULL, ...)
{
  if(any(is.na(x$breakpoints))) lbp <- 0
    else lbp <- length(x$breakpoints)
  cat(paste("\n\t Optimal ", lbp + 1, "-segment partition for `", x$fitname,"' fit: \n\n", sep = ""))
  cat("Call:\n")
  print(x$call)
  cat("\nBreakpoints at observation number:\n")
  cat(x$breakpoints,"\n")
  cat("\nCorresponding to breakdates:\n")
  cat(as.character(breakdates(x, format.times = format.times)),"\n")
}

breakdates.gbreakpoints <- function(obj, format.times = NULL, breaks = NULL, ...)
{
  if(is.null(format.times)) format.times <- !identical(obj$index, 1:obj$nobs)
  if("gbreakpointsfull" %in% class(obj)) obj <- breakpoints(obj, breaks = breaks)

  if(is.na(obj$breakpoints)[1])
    NA
  else {
    if(format.times) obj$index[obj$breakpoints] else obj$breakpoints/obj$nobs
  }
}

summary.gbreakpoints <- function(object, ...)
{
  print(object)
  cat(paste("\nObjective function `", object$objname, "': ", format(object$objective),"\n", sep = ""))
}

summary.gbreakpointsfull <- function(object, breaks = NULL,
  sort = TRUE, format.times = NULL, ...)
{
  if(is.null(format.times)) format.times <- !identical(object$index, 1:object$nobs)
  if(is.null(breaks)) breaks <- ncol(object$obj.table)/2
  n <- object$nobs
  objective <- c(object$objective(1, n), rep(NA, breaks))
  BIC <- c(AIC(breakpoints(object, breaks = 0), k = log(n)), rep(NA, breaks))
  LWZ <- c(AIC(breakpoints(object, breaks = 0), k = 0.299 * log(n)^2.1), rep(NA, breaks))
  names(objective) <- as.character(0:breaks)
  bp <- breakpoints(object, breaks = breaks)
  bd <- as.character(breakdates(bp, format.times = format.times))
  objective[breaks + 1] <- bp$objective
  BIC[breaks + 1] <- AIC(bp, k = log(n))
  LWZ[breaks + 1] <- AIC(bp, k = 0.299 * log(n)^2.1)
  bp <- bp$breakpoints
  if(breaks > 1) {
  for(m in (breaks-1):1)
  {
    bp <- rbind(NA, bp)
    bd <- rbind(NA, bd)
    bpm <- breakpoints(object, breaks = m)
    if(sort) {
      pos <- apply(outer(bpm$breakpoints, bp[nrow(bp),],
                   FUN = function(x,y) abs(x - y)), 1, which.min)
      if(length(pos) > unique(length(pos))) {
        warning("sorting not possible", call. = FALSE)
	sort <- FALSE
      }
    }
    if(!sort) pos <- 1:m
    bp[1,pos] <- bpm$breakpoints
    bd[1,pos] <- as.character(breakdates(bpm, format.times = format.times))
    objective[m+1] <- bpm$objective
    BIC[m+1] <- AIC(bpm, k = log(n))
    LWZ[m+1] <- AIC(bpm, k = 0.299 * log(n)^2.1)
  }} else {
    bp <- as.matrix(bp)
    bd <- as.matrix(bd)
  }
  rownames(bp) <- as.character(1:breaks)
  colnames(bp) <- rep("", breaks)
  rownames(bd) <- as.character(1:breaks)
  colnames(bd) <- rep("", breaks)
  objective <- rbind(objective, BIC, LWZ)
  rownames(objective) <- c(object$objname, "BIC", "LWZ")
  RVAL <- list(breakpoints = bp,
               breakdates = bd,
	       objective = objective,
	       objname = object$objname,
	       fitname = object$fitname,
	       icname = object$icname,
	       call = object$call)
  class(RVAL) <- "summary.gbreakpointsfull"
  return(RVAL)
}

print.summary.gbreakpointsfull <- function(x, ...)
{
  bp <- x$breakpoints
  breaks <- ncol(bp)
  bd <- x$breakdates
  objective <- x$objective
  bp[is.na(bp)] <- ""
  bd[is.na(bd)] <- ""
  rownames(bp) <- paste("m = ", rownames(bp), "  ", sep = "")
  rownames(bd) <- paste("m = ", rownames(bd), "  ", sep = "")
  objective <- rbind(0:(ncol(objective) - 1), format(objective))
  rownames(objective) <- c("m", x$objname, "BIC", "LWZ")
  colnames(objective) <- rep("", breaks + 1)

  cat("\n\t Optimal (m+1)-segment partition for `", x$fitname,"' : \n\n", sep = "")
  cat("Call:\n")
  print(x$call)
  cat("\nBreakpoints at observation number:\n")
  print(bp, quote = FALSE)
  cat("\nCorresponding to breakdates:\n")
  print(bd, quote = FALSE)
  cat("\nFit:\n")
  print(objective, quote = FALSE)
  invisible(x)
}

plot.gbreakpointsfull <- function(x, breaks = NULL, ...)
{
  plot(summary(x, breaks = breaks), ...)
}

plot.summary.gbreakpointsfull <- function(x, type = "b", col = c(1, 4), lty = c(1, 1),
  legend = TRUE, main = NULL, xlab = "Number of breakpoints", ylab = NULL, ...)
{
  breaks <- as.numeric(colnames(x$objective))
  obj <- x$objective[x$objname,]
  icname <- x$icname
  IC <- x$objective[icname,]
  
  col <- rep(col, length.out = 2)
  lty = rep(lty, length.out = 2)
  
  if(any(is.na(IC))) {
    if(is.null(main)) main <- ""
    if(is.null(ylab)) ylab <- x$objname
    plot(breaks, obj, ylab = ylab, xlab = xlab, main = main,
         type = type, col = col[1], lty = lty[1], ...)
  } else {
    if(is.null(main)) main <- paste(icname, "and Negative Log-Likelihood")
    if(is.null(ylab)) ylab <- ""
    plot(breaks, IC, ylab = ylab, xlab = xlab, main = main,
	 type = type, col = col[1], lty = lty[1], ...)
    par(new = TRUE)
    plot(breaks, obj, type = type, axes = FALSE, col = col[2], lty = lty[2],
         xlab = "", ylab = "")
    if(legend) legend("topright", legend = c(icname, "neg. Log-Lik."),
                      lty = lty, col = col, bty = "n")
    axis(4)
    par(new = FALSE)
  }
}


logLik.gbreakpoints <- function(object, ...)
{
  if(is.null(object$npar)) return(NA)
  n <- object$nobs
  nb <- length(object$breakpoints[!is.na(object$breakpoints)])
  df <- object$npar * (nb + 1) + nb
  logL <- -object$objective
  attr(logL, "df") <- df
  class(logL) <- "logLik"
  return(logL)
}

logLik.gbreakpointsfull <- function(object, breaks = NULL, ...)
{
  bp <- breakpoints(object, breaks = length(object$breakpoints))
  logL <- logLik(bp)
  return(logL)
}

AIC.gbreakpoints <- function(object, ..., k = 2)
{
  if(is.null(object$npar)) return(NA)
  NextMethod()
}

AIC.gbreakpointsfull <- function(object, breaks = NULL, ..., k = 2)
{
  if(is.null(object$npar)) return(NA)
  if(is.null(breaks)) breaks <- 0:(ncol(object$obj.table)/2)
  RVAL <- NULL
  for(m in breaks)
    RVAL <- c(RVAL, AIC(breakpoints(object, breaks = m), k = k))
  names(RVAL) <- breaks
  return(RVAL)
}

confint.gbreakpoints <- function(object, parm, level = 0.95, ...)
{
  stop(paste("no confidence intervals available for", dQuote("gbreakpoints"), "objects"))
}

Try the fxregime package in your browser

Any scripts or data that you put into this service are public.

fxregime documentation built on May 2, 2020, 9:05 a.m.