Nothing
utils::globalVariables(c("i", "%dopar%"))
breakpoints <- function(obj, ...)
{
UseMethod("breakpoints")
}
breakpoints.Fstats <- function(obj, ...)
{
RVAL <- list(breakpoints = obj$breakpoint,
RSS = obj$RSS,
nobs = obj$nobs,
nreg = obj$nreg,
call = match.call(),
datatsp = obj$datatsp)
class(RVAL) <- "breakpoints"
return(RVAL)
}
breakpoints.formula <- function(formula, h = 0.15, breaks = c("BIC", "LWZ", "RSS", "all"),
data = list(), hpc = c("none", "foreach"), ...)
{
mf <- model.frame(formula, data = data)
y <- model.response(mf)
modelterms <- terms(formula, data = data)
X <- model.matrix(modelterms, data = data)
RVAL <- breakpoints.matrix(X, y, h = h, breaks = breaks, hpc = hpc, ...)
n <- nrow(X)
# Not sure if this is any different from what the object returned already, but just in case...
if(is.ts(data)) {
if(NROW(data) == n) datatsp <- tsp(data)
else datatsp <- c(1/n, 1, n)
} else {
env <- environment(formula)
if(missing(data)) data <- env
orig.y <- eval(attr(terms(formula), "variables")[[2]], data, env)
if(is.ts(orig.y) & (NROW(orig.y) == n)) datatsp <- tsp(orig.y)
else datatsp <- c(1/n, 1, n)
}
RVAL$datatsp <- datatsp
return(RVAL)
}
breakpoints.matrix <- function(obj, y, h = 0.15, breaks = c("BIC", "LWZ", "RSS", "all"), hpc = c("none", "foreach"), ...)
{
X <- obj
n <- nrow(X)
k <- ncol(X)
breakstat <- NULL
intercept_only <- isTRUE(all.equal(as.vector(X), rep(1L, n)))
if(is.null(h)) h <- k + 1
if(h < 1) h <- floor(n*h)
if(h <= k)
stop("minimum segment size must be greater than the number of regressors")
if(h > floor(n/2))
stop("minimum segment size must be smaller than half the number of observations")
if (!is.numeric(breaks))
{
breakstat <- match.arg(breaks)
breaks <- ceiling(n/h) - 2
} else {
if (length(breaks) > 1)
stop("Argument 'breaks' takes a single number or method for optimal break estimation")
if (breaks %% 1 != 0)
stop("Please enter an integer number of breaks")
if(breaks < 1) {
breaks <- 1
warning("number of breaks must be at least 1")
}
if(breaks > ceiling(n/h) - 2) {
breaks0 <- breaks
breaks <- ceiling(n/h) - 2
warning(sprintf("requested number of breaks = %i too large, changed to %i", breaks0, breaks))
}
}
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 ith row of the RSS diagonal matrix, i.e,
## the recursive residuals for segments starting at i = 1:(n-h+1)
if (getOption("strucchange.use_armadillo", FALSE)) {
res = .sc_cpp_construct_rss_table(y,X,n,h,breaks,intercept_only,sqrt(.Machine$double.eps)/ncol(X),getOption("strucchange.armadillo_rcond_min",sqrt(.Machine$double.eps)))
RSS.table = res$RSS.table
dimnames(RSS.table) = list(as.character(h:(n-h)),
as.vector(rbind(paste("break", 1:breaks, sep = ""),paste("RSS", 1:breaks, sep = ""))))
RSS.triang = res$RSS.triang
RSS <- function(i, j) .sc_cpp_rss(RSS.triang, i, j)
extend.RSS.table <- function(RSS.table, breaks) {
if (2*breaks > ncol(RSS.table)) {
RSS.table = .sc_cpp_extend_rss_table(rss_table = RSS.table, rss_triang = RSS.triang, n = n, h=h, breaks = breaks)
dimnames(RSS.table) = list(as.character(h:(n-h)), as.vector(rbind(paste("break", 1:breaks, sep = ""),paste("RSS", 1:breaks, sep = ""))))
}
RSS.table
}
}
else {
RSSi <- function(i)
{
ssr <- if(intercept_only) {
(y[i:n] - cumsum(y[i:n])/(1L:(n-i+1L)))[-1L] * sqrt(1L + 1L/(1L:(n-i)))
} else {
recresid(X[i:n,,drop = FALSE],y[i:n], ...)
}
c(rep(NA, k), cumsum(ssr^2))
}
## employ HPC support if available/selected
RSS.triang <- if(hpc == "none") sapply(1:(n-h+1), RSSi) else foreach::foreach(i = 1:(n-h+1)) %dopar% RSSi(i)
## function to extract the RSS(i,j) from RSS.triang
RSS <- function(i,j) RSS.triang[[i]][j - i + 1]
## compute optimal previous partner if observation i is the mth break
## store results together with RSSs in RSS.table
## breaks = 1
index <- h:(n-h)
break.RSS <- sapply(index, function(i) RSS(1,i))
RSS.table <- cbind(index, break.RSS)
rownames(RSS.table) <- as.character(index)
## breaks >= 2
extend.RSS.table <- function(RSS.table, breaks)
{
if((breaks*2) > ncol(RSS.table)) {
for(m in (ncol(RSS.table)/2 + 1):breaks)
{
my.index <- (m*h):(n-h)
my.RSS.table <- RSS.table[,c((m-1)*2 - 1, (m-1)*2)]
my.RSS.table <- cbind(my.RSS.table, NA, NA)
for(i in my.index)
{
pot.index <- ((m-1)*h):(i - h)
break.RSS <- sapply(pot.index, function(j) my.RSS.table[as.character(j), 2] + RSS(j+1,i))
opt <- which.min(break.RSS)
my.RSS.table[as.character(i), 3:4] <- c(pot.index[opt], break.RSS[opt])
}
RSS.table <- cbind(RSS.table, my.RSS.table[,3:4])
}
colnames(RSS.table) <- as.vector(rbind(paste("break", 1:breaks, sep = ""),
paste("RSS", 1:breaks, sep = "")))
}
return(RSS.table)
}
RSS.table <- extend.RSS.table(RSS.table, breaks)
}
## extract optimal breaks
extract.breaks <- function(RSS.table, breaks)
{
if((breaks*2) > ncol(RSS.table)) stop("compute RSS.table with enough breaks before")
index <- RSS.table[, 1, drop = TRUE]
break.RSS <- sapply(index, function(i) RSS.table[as.character(i),breaks*2] + RSS(i + 1, n))
opt <- index[which.min(break.RSS)]
if(breaks > 1) {
for(i in ((breaks:2)*2 - 1))
opt <- c(RSS.table[as.character(opt[1]),i], opt)
}
names(opt) <- NULL
return(opt)
}
opt <- extract.breaks(RSS.table, breaks)
if(is.ts(y) && NROW(y) == n)
{
datatsp <- tsp(y)
} else {
datatsp <- c(1/n, 1, n)
}
RVAL <- list(breakpoints = opt,
RSS.table = RSS.table,
RSS.triang = RSS.triang,
RSS = RSS,
extract.breaks = extract.breaks,
extend.RSS.table = extend.RSS.table,
nobs = n,
nreg = k, y = y, X = X,
call = match.call(),
datatsp = datatsp)
class(RVAL) <- c("breakpointsfull", "breakpoints")
RVAL$breakpoints <- breakpoints(RVAL, breaks=breakstat)$breakpoints
return(RVAL)
}
breakpoints.breakpointsfull <- function(obj, breaks = c("BIC", "LWZ", "RSS", "all"), ...)
{
if (is.numeric(breaks))
{
if (length(breaks) > 1)
stop("This function is for extracting a single break")
if (breaks %% 1 != 0)
stop("Please enter an integer number of breaks")
} else
{
breakstat <- match.arg(breaks)
sbp <- summary(obj)
# Select optimal number of breaks by minimising a given statistic
# Note: we might want to handle cases where the difference is < 2
if (breakstat == "all")
breaks <- ncol(sbp$breakpoints)
else
breaks <- which.min(sbp$RSS[breakstat,]) - 1
}
if(breaks < 1)
{
breakpoints <- NA
RSS <- obj$RSS(1, obj$nobs)
} else {
RSS.tab <- obj$extend.RSS.table(obj$RSS.table, breaks)
breakpoints <- obj$extract.breaks(RSS.tab, breaks)
bp <- c(0, breakpoints, obj$nobs)
RSS <- sum(apply(cbind(bp[-length(bp)]+1,bp[-1]), 1,
function(x) obj$RSS(x[1], x[2])))
}
fvals = fitted(obj, breaks=breaks, bp=breakpoints)
mss = sum((fvals - mean(fvals))^2)
r.squared = mss/(mss + RSS)
RVAL <- list(breakpoints = breakpoints,
RSS = RSS,
nobs = obj$nobs,
nreg = obj$nreg,
call = match.call(),
datatsp = obj$datatsp,
r.squared = r.squared,
MSS = mss)
class(RVAL) <- "breakpoints"
return(RVAL)
}
print.breakpoints <- function(x, format.times = NULL, ...)
{
if(is.null(format.times)) format.times <- ((x$datatsp[3] > 1) & (x$datatsp[3] < x$nobs))
if(any(is.na(x$breakpoints))) lbp <- 0
else lbp <- length(x$breakpoints)
cat(paste("\n\t Optimal ", lbp + 1, "-segment partition: \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(breakdates(x, format.times = format.times),"\n")
}
breakdates <- function(obj, format.times = FALSE, ...)
{
UseMethod("breakdates")
}
breakdates.breakpoints <- function(obj, format.times = FALSE, breaks = NULL, ...)
{
if(inherits(obj, "breakpointsfull") && !is.null(breaks)) obj <- breakpoints(obj, breaks = breaks)
if(is.null(format.times)) format.times <- ((obj$datatsp[3] > 1) & (obj$datatsp[3] < obj$nobs))
format.time <- function(timevec, freq)
{
first <- floor(timevec + .001)
second <- floor(freq * (timevec - first) + 1 + .5 + .001)
RVAL <- cbind(first, second)
dummy <- function(x) paste(x[1], "(", x[2], ")", sep = "")
RVAL <- apply(RVAL, 1, dummy)
return(RVAL)
}
if(is.na(obj$breakpoints)[1])
breakdates <- NA
else {
breakdates <- (obj$breakpoints - 1)/obj$datatsp[3] + obj$datatsp[1]
if(format.times) breakdates <- format.time(breakdates, obj$datatsp[3])
}
return(breakdates)
}
breakfactor <- function(obj, breaks = NULL, labels = NULL, ...)
{
if("breakpointsfull" %in% class(obj)) obj <- breakpoints(obj, breaks = breaks)
breaks <- obj$breakpoints
if(all(is.na(breaks))) return(factor(rep("segment1", obj$nobs)))
nbreaks <- length(breaks)
fac <- rep(1:(nbreaks + 1), c(breaks[1], diff(c(breaks, obj$nobs))))
if(is.null(labels)) labels <- paste("segment", 1:(nbreaks+1), sep = "")
fac <- factor(fac, labels = labels, ...)
return(fac)
}
lines.breakpoints <- function(x, breaks = NULL, lty = 2, ...)
{
if("breakpointsfull" %in% class(x)) x <- breakpoints(x, breaks = breaks)
abline(v = breakdates(x), lty = lty, ...)
}
summary.breakpoints <- function(object, ...)
{
print(object)
cat(paste("\nRSS:", format(object$RSS),"MSS:", format(object$MSS), "\n"))
cat(paste("Multiple R-squared:", format(object$r.squared),"\n"))
}
summary.breakpointsfull <- function(object, breaks = NULL,
sort = TRUE, format.times = NULL, ...)
{
if(is.null(format.times)) format.times <- ((object$datatsp[3] > 1) & (object$datatsp[3] < object$nobs))
if(is.null(breaks)) breaks <- ncol(object$RSS.table)/2
n <- object$nobs
RSS <- c(object$RSS(1, n), rep(NA, breaks))
R.sq <- c(breakpoints(object, breaks = 0)$r.squared, rep(NA, breaks))
BIC <- c(n * (log(RSS[1]) + 1 - log(n) + log(2*pi)) + log(n) * (object$nreg + 1),
rep(NA, breaks))
names(RSS) <- as.character(0:breaks)
bp <- breakpoints(object, breaks = breaks)
bd <- breakdates(bp, format.times = format.times, breaks=breaks)
RSS[breaks + 1] <- bp$RSS
R.sq[breaks + 1] <- bp$r.squared
BIC[breaks + 1] <- AIC(bp, k = log(n))
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] <- breakdates(bpm, format.times = format.times)
RSS[m+1] <- bpm$RSS
R.sq[m+1] <- bpm$r.squared
BIC[m+1] <- AIC(bpm, k = log(n))
}} 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)
LWZ = LWZ.breakpointsfull(object)
RSS <- rbind(RSS, BIC, LWZ, R.sq)
rownames(RSS) <- c("RSS", "BIC", "LWZ", "R.sq")
RVAL <- list(breakpoints = bp,
breakdates = bd,
RSS = RSS,
call = object$call)
class(RVAL) <- "summary.breakpointsfull"
return(RVAL)
}
print.summary.breakpointsfull <- function(x, digits = max(2, getOption("digits") - 3), ...)
{
bp <- x$breakpoints
breaks <- ncol(bp)
bd <- x$breakdates
RSS <- x$RSS
bp[is.na(bp)] <- ""
bd[is.na(bd)] <- ""
rownames(bp) <- paste("m = ", rownames(bp), " ", sep = "")
rownames(bd) <- paste("m = ", rownames(bd), " ", sep = "")
RSS <- rbind(0:(ncol(RSS) - 1), format(RSS, digits = digits))
rownames(RSS) <- c("m","RSS", "BIC", "LWZ", "R.sq")
colnames(RSS) <- rep("", breaks + 1)
cat("\n\t Optimal (m+1)-segment partition: \n\n")
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(RSS, quote = FALSE)
}
plot.breakpointsfull <- function(x, breaks = NULL, ...)
{
rval <- summary(x, breaks = breaks)
plot(rval, ...)
invisible(rval)
}
plot.summary.breakpointsfull <- function(x, type = "b", col = c(1,4,5), legend = TRUE,
xlab = "Number of breakpoints", ylab = "", main = "BIC, LWZ and Residual Sum of Squares", ...)
{
breaks <- as.numeric(colnames(x$RSS))
RSS <- x$RSS["RSS",]
BIC <- x$RSS["BIC",]
LWZ <- x$RSS["LWZ",]
plot(breaks, BIC, ylab = "", ylim=c(min(c(BIC, LWZ)), max(c(BIC, LWZ))), xlab = xlab, main = main, type = type, col = col[1], ...)
points(breaks, LWZ, col=col[3], type=type)
onew <- getOption("new")
par(new = TRUE)
plot(breaks, RSS, type = type, axes = FALSE, col = col[2], xlab = "", ylab = "")
if(legend) legend("topright", c("BIC", "RSS", "LWZ"), lty = rep(1, 2), col = col, bty = "n")
axis(4)
par(new = onew)
invisible(x)
}
logLik.breakpoints <- function(object, ...)
{
n <- object$nobs
df <- (object$nreg + 1) * (length(object$breakpoints[!is.na(object$breakpoints)]) + 1)
logL <- -0.5 * n * (log(object$RSS) + 1 - log(n) + log(2 * pi))
attr(logL, "df") <- df
class(logL) <- "logLik"
return(logL)
}
logLik.breakpointsfull <- function(object, breaks = NULL, ...)
{
bp <- breakpoints(object, breaks = breaks)
logL <- logLik(bp)
return(logL)
}
AIC.breakpointsfull <- function(object, breaks = NULL, ..., k = 2)
{
if(is.null(breaks)) breaks <- 0:(ncol(object$RSS.table)/2)
RVAL <- NULL
for(m in breaks)
RVAL <- c(RVAL, AIC(breakpoints(object, breaks = m), k = k))
names(RVAL) <- breaks
return(RVAL)
}
LWZ <- function(object, ...)
{
UseMethod("LWZ")
}
LWZ.breakpointsfull <- function(object, ...)
{
return(AIC.breakpointsfull(object, ..., k=0.299 * log(object$nobs)^2.1))
}
LWZ.breakpoints <- function(object, ...)
{
return(AIC(object, k = 0.299 * log(object$nobs)^2.1))
}
pargmaxV <- function(x, xi = 1, phi1 = 1, phi2 = 1)
{
phi <- xi * (phi2/phi1)^2
G1 <- function(x, xi = 1, phi = 1)
{
x <- abs(x)
frac <- xi/phi
rval <- - exp(log(x)/2 - x/8 - log(2*pi)/2) -
(phi/xi * (phi + 2*xi)/(phi+xi)) * exp((frac * (1 + frac) * x/2) + pnorm(-(0.5 + frac) * sqrt(x), log.p = TRUE)) +
exp(log(x/2 - 2 + ((phi + 2 * xi)^2)/((phi + xi)*xi)) + pnorm(-sqrt(x)/2, log.p = TRUE))
rval
}
G2 <- function(x, xi = 1, phi = 1)
{
x <- abs(x)
frac <- xi^2/phi
rval <- 1 + sqrt(frac) * exp(log(x)/2 - (frac*x)/8 - log(2*pi)/2) +
(xi/phi * (2*phi + xi)/(phi + xi)) * exp(((phi + xi) * x/2) + pnorm(-(phi + xi/2)/sqrt(phi) * sqrt(x), log.p = TRUE)) -
exp(log(((2*phi + xi)^2)/((phi+xi)*phi) - 2 + frac*x/2) + pnorm(-sqrt(frac) * sqrt(x)/2 , log.p = TRUE))
rval
}
ifelse(x < 0, G1(x, xi = xi, phi = phi), G2(x, xi = xi, phi = phi))
}
confint.breakpointsfull <- function(object, parm = NULL, level = 0.95, breaks = NULL,
het.reg = TRUE, het.err = TRUE, vcov. = NULL, sandwich = TRUE, ...)
{
X <- object$X
y <- object$y
n <- object$nobs
a2 <- (1 - level)/2
if(!is.null(parm) & !is.null(breaks))
warning("`parm' and `breaks' are both specified: `breaks' is used")
else
if(!is.null(parm)) breaks <- parm
myfun <- function(x, level = 0.975, xi = 1, phi1 = 1, phi2 = 1)
(pargmaxV(x, xi = xi, phi1 = phi1, phi2 = phi2) - level)
myprod <- function(delta, mat) as.vector(crossprod(delta, mat) %*% delta)
bp <- breakpoints(object, breaks = breaks)$breakpoints
if(any(is.na(bp))) stop("cannot compute confidence interval when `breaks = 0'")
nbp <- length(bp)
upper <- rep(0, nbp)
lower <- rep(0, nbp)
bp <- c(0, bp, n)
res <- residuals(object, breaks = breaks)
sigma1 <- sigma2 <- sum(res^2)/n
Q1 <- Q2 <- crossprod(X)/n
if(is.null(vcov.))
Omega1 <- Omega2 <- sigma1 * Q1
else {
y.nb <- rowSums(X) + res
fm <- lm(y.nb ~ 0 + X)
if(sandwich) {
Omega1 <- Omega2 <- n * crossprod(Q1, vcov.(fm)) %*% Q1
} else {
Omega1 <- Omega2 <- vcov.(fm)
}
}
xi <- 1
X2 <- X[(bp[1]+1):bp[2],,drop = FALSE]
y2 <- y[(bp[1]+1):bp[2]]
fm2 <- lm(y2 ~ 0+ X2)
beta2 <- coef(fm2)
if(het.reg) Q2 <- crossprod(X2)/nrow(X2)
if(het.err) {
sigma2 <- sum(residuals(fm2)^2)/nrow(X2)
if(is.null(vcov.))
Omega2 <- sigma2 * Q2
else {
if(sandwich) Omega2 <- nrow(X2) * crossprod(Q2, vcov.(fm2)) %*% Q2
else Omega2 <- vcov.(fm2)
}
}
for(i in 2:(nbp+1))
{
X1 <- X2
y1 <- y2
beta1 <- beta2
sigma1 <- sigma2
Q1 <- Q2
Omega1 <- Omega2
X2 <- X[(bp[i]+1):bp[i+1],,drop = FALSE]
y2 <- y[(bp[i]+1):bp[i+1]]
fm2 <- lm(y2 ~ 0 + X2)
beta2 <- coef(fm2)
delta <- beta2 - beta1
if(het.reg) Q2 <- crossprod(X2)/nrow(X2)
if(het.err) {
sigma2 <- sum(residuals(fm2)^2)/nrow(X2)
if(is.null(vcov.))
Omega2 <- sigma2 * Q2
else {
if(sandwich) Omega2 <- nrow(X2) * crossprod(Q2, vcov.(fm2)) %*% Q2
else Omega2 <- vcov.(fm2)
}
}
Oprod1 <- myprod(delta, Omega1)
Oprod2 <- myprod(delta, Omega2)
Qprod1 <- myprod(delta, Q1)
Qprod2 <- myprod(delta, Q2)
if(het.reg) xi <- Qprod2/Qprod1
if(!is.null(vcov.)) phi1 <- sqrt(Oprod1/Qprod1)
else phi1 <- sqrt(sigma1)
if(!is.null(vcov.)) phi2 <- sqrt(Oprod2/Qprod2)
else phi2 <- sqrt(sigma2)
p0 <- pargmaxV(0, phi1 = phi1, phi2 = phi2, xi = xi)
if(is.nan(p0) || p0 < a2 || p0 > (1-a2)) {
warning(paste("Confidence interval", as.integer(i-1),
"cannot be computed: P(argmax V <= 0) =", round(p0, digits = 4)))
upper[i-1] <- NA
lower[i-1] <- NA
} else {
ub <- lb <- 0
while(pargmaxV(ub, phi1 = phi1, phi2 = phi2, xi = xi) < (1 - a2)) ub <- ub + 1000
while(pargmaxV(lb, phi1 = phi1, phi2 = phi2, xi = xi) > a2) lb <- lb - 1000
upper[i-1] <- uniroot(myfun, c(0, ub), level = (1-a2), xi = xi, phi1 = phi1, phi2 = phi2)$root
lower[i-1] <- uniroot(myfun, c(lb, 0), level = a2, xi = xi, phi1 = phi1, phi2 = phi2)$root
upper[i-1] <- upper[i-1] * phi1^2 / Qprod1
lower[i-1] <- lower[i-1] * phi1^2 / Qprod1
}
}
bp <- bp[-c(1, nbp+2)]
bp <- cbind(bp - ceiling(upper), bp, bp - floor(lower))
#V.BP# bp <- cbind(floor(bp - upper) - 1, bp, floor(bp - lower) + 1)
a2 <- round(a2 * 100, digits = 1)
colnames(bp) <- c(paste(a2, "%"), "breakpoints", paste(100 - a2, "%"))
rownames(bp) <- 1:nbp
RVAL <- list(confint = bp,
nobs = object$nobs,
nreg = object$nreg,
call = match.call(),
datatsp = object$datatsp)
class(RVAL) <- "confint.breakpoints"
return(RVAL)
}
breakdates.confint.breakpoints <- function(obj, format.times = FALSE, ...)
{
bp <- list(breakpoints = NA, nobs = obj$nobs, datatsp = obj$datatsp)
class(bp) <- "breakpoints"
RVAL <- obj$confint
for(i in 1:3) {
bp$breakpoints <- obj$confint[,i]
RVAL[,i] <- breakdates(bp, format.times = format.times, ...)
}
bp$breakpoints <- c(1, obj$nobs)
startend <- breakdates(bp, format.times = NULL, ...)
nbp <- nrow(obj$confint)
if(any(obj$confint < 1) | any(obj$confint > obj$nobs))
warning(paste("Confidence intervals outside data time interval\n\t from ",
startend[1], " to ", startend[2], " (", obj$nobs, " observations)", sep = ""), call. = FALSE)
if(any(obj$confint[-1,1] < obj$confint[-nbp,3]))
warning("Overlapping confidence intervals", call. = FALSE)
return(RVAL)
}
print.confint.breakpoints <- function(x, format.times = NULL, ...)
{
if(is.null(format.times)) format.times <- ((x$datatsp[3] > 1) & (x$datatsp[3] < x$nobs))
nbp <- nrow(x$confint)
cat("\n\t Confidence intervals for breakpoints")
cat(paste("\n\t of optimal ", nbp + 1, "-segment partition: \n\n", sep = ""))
cat("Call:\n")
print(x$call)
cat("\nBreakpoints at observation number:\n")
print(x$confint, quote = FALSE)
cat("\nCorresponding to breakdates:\n")
print(breakdates(x, format.times = format.times, ...), quote = FALSE)
}
lines.confint.breakpoints <- function(x, col = 2, angle = 90, length = 0.05,
code = 3, at = NULL, breakpoints = TRUE, ...)
{
nbp <- nrow(x$confint)
x <- breakdates(x)
if(breakpoints) abline(v = x[,2], lty = 2)
if(is.null(at)) {
at <- par("usr")[3:4]
at <- diff(at)/1.08 * 0.02 + at[1]
}
if(length(at) < nbp) at <- rep(at, length.out = nbp)
arrows(x[,1], at, x[,3], at, col = col, angle = angle, length = length, code = code, ...)
}
coef.breakpointsfull <- function(object, breaks = NULL, names = NULL, ...)
{
X <- object$X
y <- object$y
n <- object$nobs
bp <- obp <- breakpoints(object, breaks = breaks)$breakpoints
if(any(is.na(bp))) {
nbp <- 0
bp <- c(0, n)
} else {
nbp <- length(bp)
bp <- c(0, bp, n)
}
if(!is.null(names)) {
if(length(names) == 1) names <- paste(names, 1:(nbp+1))
else if(length(names) != (nbp+1)) names <- NULL
}
if(is.null(names)) {
bd1 <- structure(list(breakpoints = bp[-(nbp+2)] + 1, nobs = n, datatsp = object$datatsp),
class = "breakpoints")
bd2 <- structure(list(breakpoints = bp[-1], nobs = n, datatsp = object$datatsp),
class = "breakpoints")
bd1 <- breakdates(bd1, format.times = NULL)
bd2 <- breakdates(bd2, format.times = NULL)
names <- paste(bd1, "-", bd2)
}
rval <- NULL
for(i in 1:(nbp+1))
{
X2 <- X[(bp[i]+1):bp[i+1],,drop = FALSE]
y2 <- y[(bp[i]+1):bp[i+1]]
rval <- rbind(rval, lm.fit(X2, y2)$coef)
}
rownames(rval) <- names
return(rval)
}
fitted.breakpointsfull <- function(object, breaks = NULL, bp=NULL, ...)
{
X <- object$X
y <- object$y
n <- object$nobs
if (is.null(bp))
bp <- obp <- breakpoints(object, breaks = breaks)$breakpoints
else obp <- bp
if(any(is.na(bp))) {
nbp <- 0
bp <- c(0, n)
} else {
nbp <- length(bp)
bp <- c(0, bp, n)
}
rval <- NULL
for(i in 1:(nbp+1))
{
X2 <- X[(bp[i]+1):bp[i+1],,drop = FALSE]
y2 <- y[(bp[i]+1):bp[i+1]]
rval <- c(rval, lm.fit(X2, y2)$fitted.values)
}
rval <- ts(as.vector(rval))
tsp(rval) <- object$datatsp
return(rval)
}
residuals.breakpointsfull <- function(object, breaks = NULL, ...)
{
X <- object$X
y <- object$y
n <- object$nobs
bp <- obp <- breakpoints(object, breaks = breaks)$breakpoints
if(any(is.na(bp))) {
nbp <- 0
bp <- c(0, n)
} else {
nbp <- length(bp)
bp <- c(0, bp, n)
}
rval <- NULL
for(i in 1:(nbp+1))
{
X2 <- X[(bp[i]+1):bp[i+1],,drop = FALSE]
y2 <- y[(bp[i]+1):bp[i+1]]
rval <- c(rval, lm.fit(X2, y2)$residuals)
}
rval <- ts(as.vector(rval))
tsp(rval) <- object$datatsp
return(rval)
}
vcov.breakpointsfull <- function(object, breaks = NULL, names = NULL, het.reg = TRUE,
het.err = TRUE, vcov. = NULL, sandwich = TRUE, ...)
{
X <- object$X
y <- object$y
n <- object$nobs
bp <- breakpoints(object, breaks = breaks)$breakpoints
if(any(is.na(bp))) {
nbp <- 0
bp <- c(0, n)
} else {
nbp <- length(bp)
bp <- c(0, bp, n)
}
if(!is.null(names)) {
if(length(names) == 1) names <- paste(names, 1:(nbp+1))
else if(length(names) != (nbp+1)) names <- NULL
}
if(is.null(names)) {
bd1 <- structure(list(breakpoints = bp[-(nbp+2)] + 1, nobs = n, datatsp = object$datatsp),
class = "breakpoints")
bd2 <- structure(list(breakpoints = bp[-1], nobs = n, datatsp = object$datatsp),
class = "breakpoints")
bd1 <- breakdates(bd1, format.times = NULL)
bd2 <- breakdates(bd2, format.times = NULL)
names <- paste(bd1, "-", bd2)
}
res <- residuals(object, breaks = breaks)
sigma2 <- sum(res^2)/n
Q2 <- crossprod(X)/n
if(is.null(vcov.))
Omega2 <- sigma2 * solve(Q2) / n
else {
y.nb <- rowSums(X) + res
fm <- lm(y.nb ~ 0 + X)
if(sandwich) {
Omega2 <- vcov.(fm)
} else {
modelv <- summary(fm)$cov.unscaled
Omega2 <- n * modelv %*% vcov.(fm) %*% modelv
}
}
rownames(Omega2) <- colnames(Omega2) <- colnames(X)
rval <- list()
for(i in 1:(nbp+1))
{
X2 <- X[(bp[i]+1):bp[i+1],,drop = FALSE]
y2 <- y[(bp[i]+1):bp[i+1]]
fm2 <- lm(y2 ~ 0 + X2)
if(het.reg) Q2 <- crossprod(X2)/nrow(X2)
if(het.err) {
sigma2 <- sum(residuals(fm2)^2)/nrow(X2)
if(is.null(vcov.))
Omega2 <- sigma2 * solve(Q2) / nrow(X2)
else {
if(sandwich) {
Omega2 <- vcov.(fm2)
} else {
modelv <- summary(fm2)$cov.unscaled
Omega2 <- n * modelv %*% vcov.(fm2) %*% modelv
}
}
rownames(Omega2) <- colnames(Omega2) <- colnames(X)
}
rval[[i]] <- Omega2
}
names(rval) <- names
return(rval)
}
df.residual.breakpointsfull <- function(object, ...)
{
rval <- table(breakfactor(object, ...)) - object$nreg
names(rval) <- rownames(coef(object, ...))
return(rval)
}
magnitude <- function(object, ...)
{
UseMethod("magnitude")
}
# Returns a vector of magnitudes of change
magnitude.breakpointsfull <- function(object, interval = 0.1, breaks = NULL, component = "trend", ...)
{
X <- object$X[,!colnames(object$X) %in% "(Intercept)", drop=FALSE] # Do not take intercept
y <- object$y
# Also filter out the intercept from the components, in case users are lazy and put in all model names()
component <- component[!component %in% "(Intercept)"]
if (interval <= 0 || interval > length(y))
stop("Requested interval for magnitude computation out of valid range")
if (interval < 1)
interval <- floor(length(y)*interval) # Convert to number of samples #TODO: Add handling of time
bp <- breakpoints(object, breaks=breaks)$breakpoints
nrbp <- length(bp)
if (nrbp < 2 && is.na(bp))
stop("There are no breakpoints to calculate magnitudes for!")
if (!any(colnames(object$X) %in% component))
stop(paste("The specified component", component, "is missing"))
co <- coef(object, breaks=breaks)
Mag <- matrix(NA, nrbp, 6)
for (i in 1:nrbp) {
interval_start <- max(bp[i]-interval, 1)
interval_end <- min(bp[i]+interval, nrow(X))
# Fitted components over the interval range from the breakpoint
fit_prev <- co[i, "(Intercept)"]
fit_next <- co[i+1, "(Intercept)"]
for (comp in component) {
fit_prev <- X[interval_start:interval_end,comp] * co[i, comp] + fit_prev
fit_next <- X[interval_start:interval_end,comp] * co[i+1, comp] + fit_next
}
# First fitted values before and after the break; for legacy reasons
Mag[i, 1] <- co[i, "(Intercept)"]
Mag[i, 2] <- co[i+1, "(Intercept)"]
for (comp in component) {
Mag[i, 1] <- X[bp[i], comp] * co[i, comp] + Mag[i, 1]
Mag[i, 2] <- X[bp[i]+1,comp] * co[i+1, comp] + Mag[i, 2]
}
Mag[i, 3] <- Mag[i, 2] - Mag[i, 1]
Mag[i, 4] <- sqrt(mean((fit_next - fit_prev)^2))
Mag[i, 5] <- mean(abs(fit_next - fit_prev))
Mag[i, 6] <- mean(fit_next - fit_prev)
colnames(Mag) = c("before", "after", "diff", "RMSD", "MAD", "MD")
}
index <- which.max(abs(Mag[, 3]))
m.x <- rep(bp[index], 2)
m.y <- c(Mag[index, 1], Mag[index, 2]) #Magnitude position
Magnitude <- Mag[index, 3] # Magnitude of biggest change
Time <- bp[index]
Result <- list(Mag=Mag, m.x=m.x, m.y=m.y, Magnitude=Magnitude, Time=Time)
class(Result) <- "magnitude"
return(Result)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.