From 54bc161a2881bda2cac8213a99d92acda81b0759 Mon Sep 17 00:00:00 2001 From: Ryan Hafen Date: Mon, 6 Apr 2026 10:16:22 -0700 Subject: [PATCH 1/5] Run air format --- R/accessors.R | 29 +-- R/cpp_wrappers.R | 2 +- R/loess_stl.R | 128 ++++++++--- R/plots.R | 194 ++++++++++------ R/stlplus.R | 478 ++++++++++++++++++++++++++++----------- man-roxygen/ex-stlplus.R | 130 ++++++++--- tests/testthat/tests.R | 59 +++-- 7 files changed, 731 insertions(+), 289 deletions(-) diff --git a/R/accessors.R b/R/accessors.R index 8bcf012..b3a7954 100755 --- a/R/accessors.R +++ b/R/accessors.R @@ -38,7 +38,7 @@ getraw <- function(object) { #' @export #' @rdname accessors remainder.stlplus <- function(object) { - if(object$pars$fc.number == 0) { + if (object$pars$fc.number == 0) { object$data$remainder } else { object$fc$remainder @@ -48,20 +48,20 @@ remainder.stlplus <- function(object) { #' @export #' @rdname accessors fitted.stlplus <- function(object, ...) { - if(object$pars$fc.number == 0) { + if (object$pars$fc.number == 0) { object$data$seasonal + object$data$trend } else { - object$data$seasonal + apply(object$fc[,1:object$pars$fc.number], 1, sum) + object$data$seasonal + apply(object$fc[, 1:object$pars$fc.number], 1, sum) } } #' @export #' @rdname accessors predict.stlplus <- function(object, ...) { - if(object$pars$fc.number == 0) { + if (object$pars$fc.number == 0) { object$data$seasonal + object$data$trend } else { - object$data$seasonal + apply(object$fc[,1:object$pars$fc.number], 1, sum) + object$data$seasonal + apply(object$fc[, 1:object$pars$fc.number], 1, sum) } } @@ -81,23 +81,26 @@ trend.stlplus <- function(object) { #' @export #' @rdname accessors fc <- function(object, fcnum = 1) { - if(! "stlplus" %in% class(object)) + if (!"stlplus" %in% class(object)) { stop("object not of class stlplus") + } - if(is.null(object$fc)) + if (is.null(object$fc)) { stop("there are no post-trend frequency components") + } - if(fcnum > ncol(object$fc)) + if (fcnum > ncol(object$fc)) { stop("there are not that many frequency components") + } - object$fc[,fcnum] + object$fc[, fcnum] } # setGeneric("time") #' @export #' @rdname accessors time.stlplus <- function(x, ...) { - if(length(x$t) == x$n) { + if (length(x$t) == x$n) { x$t } else { c(1:x$n) @@ -109,19 +112,19 @@ time.stlplus <- function(x, ...) { #' @export #' @rdname accessors remainder.stl <- function(object) { - as.numeric(object$time.series[,3]) + as.numeric(object$time.series[, 3]) } #' @export #' @rdname accessors seasonal.stl <- function(object) { - as.numeric(object$time.series[,1]) + as.numeric(object$time.series[, 1]) } #' @export #' @rdname accessors trend.stl <- function(object) { - as.numeric(object$time.series[,2]) + as.numeric(object$time.series[, 2]) } # setGeneric("time") diff --git a/R/cpp_wrappers.R b/R/cpp_wrappers.R index e4b444d..92d7c8e 100644 --- a/R/cpp_wrappers.R +++ b/R/cpp_wrappers.R @@ -1,5 +1,5 @@ .interp <- function(m, fits, slopes, at) { - if(any(is.nan(fits))) { + if (any(is.nan(fits))) { ind <- !is.nan(fits) c_interp(m[ind], fits[ind], slopes[ind], at) } else { diff --git a/R/loess_stl.R b/R/loess_stl.R index 84d7c6a..3f9db01 100755 --- a/R/loess_stl.R +++ b/R/loess_stl.R @@ -1,10 +1,19 @@ #' @importFrom yaImpute ann #' @importFrom Rcpp sourceCpp #' @useDynLib stlplus -.loess_stlplus <- function(x = NULL, y, span, degree, weights = NULL, - m = c(1:length(y)), y_idx = !is.na(y), noNA = all(y_idx), blend = 0, - jump = ceiling(span / 10), at = c(1:length(y))) { - +.loess_stlplus <- function( + x = NULL, + y, + span, + degree, + weights = NULL, + m = c(1:length(y)), + y_idx = !is.na(y), + noNA = all(y_idx), + blend = 0, + jump = ceiling(span / 10), + at = c(1:length(y)) +) { nextodd <- function(x) { x <- round(x) x2 <- ifelse(x %% 2 == 0, x + 1, x) @@ -12,26 +21,38 @@ } n <- length(y[y_idx]) - if(is.null(x)) x <- c(1:length(y)) + if (is.null(x)) { + x <- c(1:length(y)) + } - if(is.null(weights)) weights <- rep(1, length(y)) + if (is.null(weights)) { + weights <- rep(1, length(y)) + } n_m <- length(m) - if((span %% 2) == 0) { + if ((span %% 2) == 0) { span <- span + 1 - warning(paste("Span must be odd! Changed span from ", - span - 1, " to ", span, sep = "")) + warning(paste( + "Span must be odd! Changed span from ", + span - 1, + " to ", + span, + sep = "" + )) } s2 <- (span + 1) / 2 # set up indices in R - easier - if(noNA) { - if((diff(range(x))) < span) { + if (noNA) { + if ((diff(range(x))) < span) { l_idx <- rep(1, n_m) r_idx <- rep(n, n_m) - } else{ - l_idx <- c(rep(1, length(m[m < s2])), m[m >= s2 & m <= n - s2] - s2 + 1, - rep(n - span + 1, length(m[m > n - s2]))) + } else { + l_idx <- c( + rep(1, length(m[m < s2])), + m[m >= s2 & m <= n - s2] - s2 + 1, + rep(n - span + 1, length(m[m > n - s2])) + ) r_idx <- l_idx + span - 1 } aa <- abs(m - x[l_idx]) @@ -43,32 +64,50 @@ x2 <- x[y_idx] # another approach - a <- yaImpute::ann(ref = as.matrix(x2), target = as.matrix(m), tree.type = "kd", - k = span3, eps = 0, verbose = FALSE)$knnIndexDist[, 1:span3, drop = FALSE] + a <- yaImpute::ann( + ref = as.matrix(x2), + target = as.matrix(m), + tree.type = "kd", + k = span3, + eps = 0, + verbose = FALSE + )$knnIndexDist[, 1:span3, drop = FALSE] l_idx <- apply(a, 1, min) r_idx <- apply(a, 1, max) max_dist <- apply(cbind(abs(m - x2[l_idx]), abs(x2[r_idx] - m)), 1, max) } - if(span >= n) + if (span >= n) { # max_dist <- max_dist * (span / n) max_dist <- max_dist + (span - n) / 2 + } - out <- c_loess(x[y_idx], y[y_idx], degree, span, weights[y_idx], - m, l_idx - 1, as.double(max_dist)) + out <- c_loess( + x[y_idx], + y[y_idx], + degree, + span, + weights[y_idx], + m, + l_idx - 1, + as.double(max_dist) + ) res1 <- out$result # do interpolation - if(jump > 1) + if (jump > 1) { res1 <- .interp(m, out$result, out$slope, at) - # res1 <- approx(x = m, y = out$result, xout = at)$y + } + # res1 <- approx(x = m, y = out$result, xout = at)$y - if(blend > 0 && blend <= 1 && degree >= 1) { - if(degree == 2) + if (blend > 0 && blend <= 1 && degree >= 1) { + if (degree == 2) { sp0 <- nextodd((span + 1) / 2) - if(degree == 1) + } + if (degree == 1) { sp0 <- span + } n.b <- as.integer(span / 2) @@ -76,7 +115,7 @@ # indices for left and right blending points # take into account if n_m is too small mid <- median(m) - bl_idx <- m <= n.b + jump & m < mid + bl_idx <- m <= n.b + jump & m < mid br_idx <- m >= n - n.b - jump + 1 & m >= mid left <- m[bl_idx] right <- m[br_idx] @@ -97,16 +136,30 @@ # right now, a lot of unnecessary calculation is done at the interior # where blending doesn't matter - tmp <- c_loess(x[y_idx], y[y_idx], 0, sp0, weights[y_idx], - m2, l_idx2-1, max_dist2) - - if(jump > 1) { - res2_left <- .interp(left, + tmp <- c_loess( + x[y_idx], + y[y_idx], + 0, + sp0, + weights[y_idx], + m2, + l_idx2 - 1, + max_dist2 + ) + + if (jump > 1) { + res2_left <- .interp( + left, head(tmp$result, length(left)), - head(tmp$slope, length(left)), left_interp) - res2_right <- .interp(right, + head(tmp$slope, length(left)), + left_interp + ) + res2_right <- .interp( + right, tail(tmp$result, length(right)), - tail(tmp$slope, length(right)), right_interp) + tail(tmp$slope, length(right)), + right_interp + ) } else { res2_left <- head(tmp$result, length(left)) res2_right <- tail(tmp$result, length(right)) @@ -120,8 +173,12 @@ p.right[p.right < blend] <- blend p.right[p.right > 1] <- 1 - res1[bl_idx_interp] <- res1[bl_idx_interp] * p.left + res2_left * (1 - p.left) - res1[br_idx_interp] <- res1[br_idx_interp] * p.right + res2_right * (1 - p.right) + res1[bl_idx_interp] <- res1[bl_idx_interp] * + p.left + + res2_left * (1 - p.left) + res1[br_idx_interp] <- res1[br_idx_interp] * + p.right + + res2_right * (1 - p.right) # xxx <- x[y_idx] # yyy <- y[y_idx] @@ -131,4 +188,3 @@ res1 } - diff --git a/R/plots.R b/R/plots.R index 3348fb6..417f88d 100755 --- a/R/plots.R +++ b/R/plots.R @@ -10,85 +10,108 @@ #' @seealso \code{\link{stlplus}} #' @export #' @importFrom lattice xyplot panel.abline panel.xyplot panel.loess panel.segments -plot.stlplus <- function(x, +plot.stlplus <- function( + x, scales = list(y = list(relation = "sliced")), - type = "l", as.table = TRUE, strip = FALSE, strip.left = TRUE, - between = list(y = 0.5), layout = NULL, ...) { - - if(x$pars$fc.number == 0) { + type = "l", + as.table = TRUE, + strip = FALSE, + strip.left = TRUE, + between = list(y = 0.5), + layout = NULL, + ... +) { + if (x$pars$fc.number == 0) { d <- data.frame( time = rep(time.stlplus(x), 4), values = c(getraw(x), seasonal(x), trend(x), remainder(x)), - ind = factor(rep(c(1:4), each = x$n))) + ind = factor(rep(c(1:4), each = x$n)) + ) levels(d$ind) <- c("raw", "seasonal", "trend", "remainder") d$which <- "isnotNA" - if(is.null(layout)) layout <- c(1, 4) + if (is.null(layout)) { + layout <- c(1, 4) + } # if there are any NA values, plot the part of the seasonal and trend # components with a different color - if(any(is.na(getraw(x)))) { + if (any(is.na(getraw(x)))) { d$values[d$ind == "seasonal" & is.na(getraw(x))] <- NA d$values[d$ind == "trend" & is.na(getraw(x))] <- NA - d <- rbind(d, + d <- rbind( + d, data.frame( time = rep(time.stlplus(x), 2), values = c(seasonal(x), trend(x)), ind = rep(c("seasonal", "trend"), each = x$n), - which = "isxNA") + which = "isxNA" + ) ) - d$values[d$ind == "seasonal" & d$which == "isxNA" & !is.na(getraw(x))] <- NA + d$values[ + d$ind == "seasonal" & d$which == "isxNA" & !is.na(getraw(x)) + ] <- NA d$values[d$ind == "trend" & d$which == "isxNA" & !is.na(getraw(x))] <- NA } - } else { fc.number <- x$pars$fc.number nvar <- 3 + fc.number fc.name <- as.character(x$pars$fc$fc.name) - if(fc.number == 1) { - fcdat <- x$fc[,1] + if (fc.number == 1) { + fcdat <- x$fc[, 1] } else { - fcdat <- stack(x$fc[,1:fc.number])$values + fcdat <- stack(x$fc[, 1:fc.number])$values } d <- data.frame( time = rep(time.stlplus(x), nvar), values = c(getraw(x), seasonal(x), fcdat, remainder(x)), - ind = factor(rep(c(1:nvar), each = x$n))) + ind = factor(rep(c(1:nvar), each = x$n)) + ) levels(d$ind) <- c("raw", "seasonal", fc.name, "remainder") d$which <- "isnotNA" - if(is.null(layout)) layout <- c(1, nvar) + if (is.null(layout)) { + layout <- c(1, nvar) + } # if there are any NA values, plot the part of the seasonal and trend # components with a different color - if(any(is.na(getraw(x)))) { + if (any(is.na(getraw(x)))) { d$values[d$ind == "seasonal" & is.na(getraw(x))] <- NA - for(i in 1:fc.number) { + for (i in 1:fc.number) { d$values[d$ind == fc.name[i] & is.na(getraw(x))] <- NA } - d <- rbind(d, + d <- rbind( + d, data.frame( time = rep(time.stlplus(x), fc.number + 1), values = c(seasonal(x), fcdat), ind = rep(c("seasonal", fc.name), each = x$n), - which = "isxNA") + which = "isxNA" + ) ) - d$values[d$ind == "seasonal" & d$which == "isxNA" & !is.na(getraw(x))] <- NA - for(i in 1:fc.number) { - d$values[d$ind == fc.name[i] & d$which == "isxNA" & !is.na(getraw(x))] <- NA + d$values[ + d$ind == "seasonal" & d$which == "isxNA" & !is.na(getraw(x)) + ] <- NA + for (i in 1:fc.number) { + d$values[ + d$ind == fc.name[i] & d$which == "isxNA" & !is.na(getraw(x)) + ] <- NA } } } - p <- lattice::xyplot(values ~ time | ind, + p <- lattice::xyplot( + values ~ time | ind, data = d, groups = which, type = type, layout = layout, scales = scales, as.table = as.table, - strip = strip, strip.left = strip.left, + strip = strip, + strip.left = strip.left, between = between, ... ) @@ -112,19 +135,25 @@ plot.stlplus <- function(x, #' @references R. B. Cleveland, W. S. Cleveland, J. E. McRae, and I. Terpenning (1990) STL: A Seasonal-Trend Decomposition Procedure Based on Loess. \emph{Journal of Official Statistics}, \bold{6}, 3--73. #' @seealso \code{\link{stlplus}} #' @export -plot_cycle <- function(x, layout = c(x$pars$n.p, 1), - col = "#0080ff", xlab = "Time", ylab = "Seasonal", - panel= function(x, y, ...) { +plot_cycle <- function( + x, + layout = c(x$pars$n.p, 1), + col = "#0080ff", + xlab = "Time", + ylab = "Seasonal", + panel = function(x, y, ...) { lattice::panel.segments(x, rep(.midmean(y), length(x)), x, y, col = col) - }, ...) { - + }, + ... +) { seas <- seasonal(x) t <- time.stlplus(x) cycleSubIndices <- x$data$sub.labels - if(x$pars$periodic) { - p <- lattice::xyplot(seas ~ t | cycleSubIndices, + if (x$pars$periodic) { + p <- lattice::xyplot( + seas ~ t | cycleSubIndices, layout = layout, type = "l", xlab = xlab, @@ -132,7 +161,8 @@ plot_cycle <- function(x, layout = c(x$pars$n.p, 1), ... ) } else { - p <- lattice::xyplot(seas ~ t | cycleSubIndices, + p <- lattice::xyplot( + seas ~ t | cycleSubIndices, layout = layout, type = "l", panel = panel, @@ -159,14 +189,22 @@ plot_cycle <- function(x, layout = c(x$pars$n.p, 1), #' @references R. B. Cleveland, W. S. Cleveland, J. E. McRae, and I. Terpenning (1990) STL: A Seasonal-Trend Decomposition Procedure Based on Loess. \emph{Journal of Official Statistics}, \bold{6}, 3--73. #' @seealso \code{\link{stlplus}} #' @export -plot_seasonal <- function(x, col = c("darkgray", "black"), - lwd = 2, xlab = "Time", ylab = "Centered Seasonal + Remainder", ...) { - +plot_seasonal <- function( + x, + col = c("darkgray", "black"), + lwd = 2, + xlab = "Time", + ylab = "Centered Seasonal + Remainder", + ... +) { dat <- by( data.frame( v1 = seasonal(x) + remainder(x), - v2 = seasonal(x), t = time.stlplus(x), - v3 = x$data$sub.labels), list(x$data$sub.labels), + v2 = seasonal(x), + t = time.stlplus(x), + v3 = x$data$sub.labels + ), + list(x$data$sub.labels), function(dd) { mn <- mean(dd$v1, na.rm = TRUE) data.frame(a = dd$v1 - mn, b = dd$v2 - mn, t = dd$t, sub.labels = dd$v3) @@ -175,7 +213,9 @@ plot_seasonal <- function(x, col = c("darkgray", "black"), dat <- do.call(rbind, dat) - p <- lattice::xyplot(a + b ~ t | sub.labels, data = dat, + p <- lattice::xyplot( + a + b ~ t | sub.labels, + data = dat, type = c("p", "l"), col = col, lwd = lwd, @@ -199,17 +239,25 @@ plot_seasonal <- function(x, col = c("darkgray", "black"), #' @references R. B. Cleveland, W. S. Cleveland, J. E. McRae, and I. Terpenning (1990) STL: A Seasonal-Trend Decomposition Procedure Based on Loess. \emph{Journal of Official Statistics}, \bold{6}, 3--73. #' @seealso \code{\link{stlplus}} #' @export -plot_rembycycle <- function(x, col = "darkgray", locol = "black", - lolwd = 2, xlab = "Time", ylab = "Remainder", ...) { - +plot_rembycycle <- function( + x, + col = "darkgray", + locol = "black", + lolwd = 2, + xlab = "Time", + ylab = "Remainder", + ... +) { vals2 <- data.frame( values = remainder(x), ind = x$data$sub.labels, - t = time.stlplus(x)) + t = time.stlplus(x) + ) - vals2 <- vals2[!is.na(vals2$values),] + vals2 <- vals2[!is.na(vals2$values), ] - p <- lattice::xyplot(values ~ t | ind, + p <- lattice::xyplot( + values ~ t | ind, data = vals2, type = "p", panel = function(x, y, ...) { @@ -242,39 +290,58 @@ plot_rembycycle <- function(x, col = "darkgray", locol = "black", #' @references R. B. Cleveland, W. S. Cleveland, J. E. McRae, and I. Terpenning (1990) STL: A Seasonal-Trend Decomposition Procedure Based on Loess. \emph{Journal of Official Statistics}, \bold{6}, 3--73. #' @seealso \code{\link{stlplus}} #' @export -plot_trend <- function(x, xlab = "Time", ylab = "Trend", - span = 0.3, type = c("p", "l"), +plot_trend <- function( + x, + xlab = "Time", + ylab = "Trend", + span = 0.3, + type = c("p", "l"), scales = list(y = list(relation = "free")), - lwd = c(1, 1), col = c("darkgray", "black", "darkgray"), - layout = c(1, 2), between = list(y = 0.5), - strip = FALSE, strip.left = TRUE, as.table = TRUE, ...) { - + lwd = c(1, 1), + col = c("darkgray", "black", "darkgray"), + layout = c(1, 2), + between = list(y = 0.5), + strip = FALSE, + strip.left = TRUE, + as.table = TRUE, + ... +) { dat <- rbind( data.frame( x = time.stlplus(x), y = remainder(x) + trend(x), type = "p", - pan = "Trend"), + pan = "Trend" + ), data.frame( x = time.stlplus(x), y = trend(x), type = "l", - pan = "Trend"), + pan = "Trend" + ), data.frame( x = time.stlplus(x), y = remainder(x), type = "p", - pan = "Remainder"), + pan = "Remainder" + ), data.frame( x = time.stlplus(x), - y = predict(loess(remainder(x) ~ c(1:length(remainder(x))), - span = span, weights = x$data$weights), - newdata = c(1:length(remainder(x)))), + y = predict( + loess( + remainder(x) ~ c(1:length(remainder(x))), + span = span, + weights = x$data$weights + ), + newdata = c(1:length(remainder(x))) + ), type = "l", - pan = "Remainder") + pan = "Remainder" + ) ) - p <- lattice::xyplot(y ~ x | pan, + p <- lattice::xyplot( + y ~ x | pan, groups = type, data = dat, type = type, @@ -285,11 +352,12 @@ plot_trend <- function(x, xlab = "Time", ylab = "Trend", scales = scales, distribute.type = TRUE, between = between, - strip = strip, strip.left = strip.left, - xlab = xlab, ylab = ylab, + strip = strip, + strip.left = strip.left, + xlab = xlab, + ylab = ylab, ... ) p } - diff --git a/R/stlplus.R b/R/stlplus.R index 247da98..173b261 100755 --- a/R/stlplus.R +++ b/R/stlplus.R @@ -45,60 +45,115 @@ #' @importFrom utils head stack tail #' @export #' @rdname stlplus -stlplus <- function(x, t = NULL, n.p, s.window, s.degree = 1, - t.window = NULL, t.degree = 1, - fc.window = NULL, fc.degree = NULL, fc.name = NULL, - l.window = NULL, l.degree = t.degree, +stlplus <- function( + x, + t = NULL, + n.p, + s.window, + s.degree = 1, + t.window = NULL, + t.degree = 1, + fc.window = NULL, + fc.degree = NULL, + fc.name = NULL, + l.window = NULL, + l.degree = t.degree, s.jump = ceiling(s.window / 10), t.jump = ceiling(t.window / 10), l.jump = ceiling(l.window / 10), - fc.jump = NULL, critfreq = 0.05, - s.blend = 0, t.blend = 0, l.blend = t.blend, - fc.blend = NULL, inner = 2, outer = 1, - sub.labels = NULL, sub.start = 1, zero.weight = 1e-6, - details = FALSE, ...) { - + fc.jump = NULL, + critfreq = 0.05, + s.blend = 0, + t.blend = 0, + l.blend = t.blend, + fc.blend = NULL, + inner = 2, + outer = 1, + sub.labels = NULL, + sub.start = 1, + zero.weight = 1e-6, + details = FALSE, + ... +) { UseMethod("stlplus") } #' @export #' @rdname stlplus -stlplus.ts <- function(x, t = as.numeric(stats::time(x)), n.p = frequency(x), - s.window, s.degree = 1, - t.window = NULL, t.degree = 1, - fc.window = NULL, fc.degree = NULL, fc.name = NULL, - l.window = NULL, l.degree = t.degree, +stlplus.ts <- function( + x, + t = as.numeric(stats::time(x)), + n.p = frequency(x), + s.window, + s.degree = 1, + t.window = NULL, + t.degree = 1, + fc.window = NULL, + fc.degree = NULL, + fc.name = NULL, + l.window = NULL, + l.degree = t.degree, s.jump = ceiling(s.window / 10), t.jump = ceiling(t.window / 10), l.jump = ceiling(l.window / 10), - fc.jump = NULL, critfreq = 0.05, - s.blend = 0, t.blend = 0, l.blend = t.blend, - fc.blend = NULL, inner = 2, outer = 1, - sub.labels = NULL, sub.start = 1, zero.weight = 1e-6, - details = FALSE, ...) { - - if (is.matrix(x)) + fc.jump = NULL, + critfreq = 0.05, + s.blend = 0, + t.blend = 0, + l.blend = t.blend, + fc.blend = NULL, + inner = 2, + outer = 1, + sub.labels = NULL, + sub.start = 1, + zero.weight = 1e-6, + details = FALSE, + ... +) { + if (is.matrix(x)) { stop("only univariate series are allowed") + } - if(missing(n.p)) n.p <- frequency(x) - if(!is.null(t)) { - if(length(t) != length(x)) + if (missing(n.p)) { + n.p <- frequency(x) + } + if (!is.null(t)) { + if (length(t) != length(x)) { stop("t must be same length as time series") + } } else { t <- as.vector(stats::time(x)) } - stlplus.default(x, t = t, n.p = n.p, - s.window = s.window, s.degree = s.degree, - t.window = t.window, t.degree = t.degree, - fc.window = fc.window, fc.degree = fc.degree, fc.name = fc.name, - l.window = l.window, l.degree = l.degree, - s.jump = s.jump, t.jump = t.jump, l.jump = l.jump, - fc.jump = fc.jump, critfreq = 0.05, - s.blend = s.blend, t.blend = t.blend, l.blend = l.blend, - fc.blend = NULL, inner = inner, outer = outer, - sub.labels = sub.labels, sub.start = sub.start, - details = details, ...) + stlplus.default( + x, + t = t, + n.p = n.p, + s.window = s.window, + s.degree = s.degree, + t.window = t.window, + t.degree = t.degree, + fc.window = fc.window, + fc.degree = fc.degree, + fc.name = fc.name, + l.window = l.window, + l.degree = l.degree, + s.jump = s.jump, + t.jump = t.jump, + l.jump = l.jump, + fc.jump = fc.jump, + critfreq = 0.05, + s.blend = s.blend, + t.blend = t.blend, + l.blend = l.blend, + fc.blend = NULL, + inner = inner, + outer = outer, + sub.labels = sub.labels, + sub.start = sub.start, + details = details, + ... + ) } #' @export @@ -108,24 +163,51 @@ stlplus.zoo <- function(...) { } #' @export -stlplus.default <- function(x, t = NULL, n.p, s.window, s.degree = 1, - t.window = NULL, t.degree = 1, - fc.window = NULL, fc.degree = NULL, fc.name = NULL, - l.window = NULL, l.degree = t.degree, +stlplus.default <- function( + x, + t = NULL, + n.p, + s.window, + s.degree = 1, + t.window = NULL, + t.degree = 1, + fc.window = NULL, + fc.degree = NULL, + fc.name = NULL, + l.window = NULL, + l.degree = t.degree, s.jump = ceiling(s.window / 10), t.jump = ceiling(t.window / 10), l.jump = ceiling(l.window / 10), - fc.jump = NULL, critfreq = 0.05, - s.blend = 0, t.blend = 0, l.blend = t.blend, - fc.blend = NULL, inner = 2, outer = 1, - sub.labels = NULL, sub.start = 1, zero.weight = 1e-6, - details = FALSE, ...) { - - if(missing(n.p)) stop("must specify periodicity of seasonal (either explicitly or through a time series object)") + fc.jump = NULL, + critfreq = 0.05, + s.blend = 0, + t.blend = 0, + l.blend = t.blend, + fc.blend = NULL, + inner = 2, + outer = 1, + sub.labels = NULL, + sub.start = 1, + zero.weight = 1e-6, + details = FALSE, + ... +) { + if (missing(n.p)) { + stop( + "must specify periodicity of seasonal (either explicitly or through a time series object)" + ) + } n.p <- as.integer(n.p) - if(n.p < 4) - stop(paste("Parameter n.p was set to ", n.p, ". Must be at least 4.", sep = "")) + if (n.p < 4) { + stop(paste( + "Parameter n.p was set to ", + n.p, + ". Must be at least 4.", + sep = "" + )) + } # family <- ifelse(robust, "symmetric", "gaussian") @@ -144,47 +226,67 @@ stlplus.default <- function(x, t = NULL, n.p, s.window, s.degree = 1, as.integer(x2) } - wincheck <- function(x) { x <- nextodd(x) - if(any(x <= 0)) stop("Window lengths must be positive.") + if (any(x <= 0)) { + stop("Window lengths must be positive.") + } x } degcheck <- function(x) { - if(! all(x == 0 | x == 1 | x == 2)) stop("Smoothing degree must be 0, 1, or 2") + if (!all(x == 0 | x == 1 | x == 2)) { + stop("Smoothing degree must be 0, 1, or 2") + } } get.t.window <- function(t.dg, s.dg, n.s, n.p, omega) { - if(t.dg == 0) t.dg <- 1 - if(s.dg == 0) s.dg <- 1 + if (t.dg == 0) { + t.dg <- 1 + } + if (s.dg == 0) { + s.dg <- 1 + } coefs_a <- data.frame( a = c(0.000103350651767650, 3.81086166990428e-6), - b = c(-0.000216653946625270, 0.000708495976681902)) + b = c(-0.000216653946625270, 0.000708495976681902) + ) coefs_b <- data.frame( a = c(1.42686036792937, 2.24089552678906), b = c(-3.1503819836694, -3.30435316073732), - c = c(5.07481807116087, 5.08099438760489)) + c = c(5.07481807116087, 5.08099438760489) + ) coefs_c <- data.frame( a = c(1.66534145060448, 2.33114333880815), b = c(-3.87719398039131, -1.8314816166323), - c = c(6.46952900183769, 1.85431548427732)) + c = c(6.46952900183769, 1.85431548427732) + ) # estimate critical frequency for seasonal betac0 <- coefs_a$a[s.dg] + coefs_a$b[s.dg] * omega - betac1 <- coefs_b$a[s.dg] + coefs_b$b[s.dg] * omega + coefs_b$c[s.dg] * omega^2 - betac2 <- coefs_c$a[s.dg] + coefs_c$b[s.dg] * omega + coefs_c$c[s.dg] * omega^2 + betac1 <- coefs_b$a[s.dg] + + coefs_b$b[s.dg] * omega + + coefs_b$c[s.dg] * omega^2 + betac2 <- coefs_c$a[s.dg] + + coefs_c$b[s.dg] * omega + + coefs_c$c[s.dg] * omega^2 f_c <- (1 - (betac0 + betac1 / n.s + betac2 / n.s^2)) / n.p # choose betat0 <- coefs_a$a[t.dg] + coefs_a$b[t.dg] * omega - betat1 <- coefs_b$a[t.dg] + coefs_b$b[t.dg] * omega + coefs_b$c[t.dg] * omega^2 - betat2 <- coefs_c$a[t.dg] + coefs_c$b[t.dg] * omega + coefs_c$c[t.dg] * omega^2 + betat1 <- coefs_b$a[t.dg] + + coefs_b$b[t.dg] * omega + + coefs_b$c[t.dg] * omega^2 + betat2 <- coefs_c$a[t.dg] + + coefs_c$b[t.dg] * omega + + coefs_c$c[t.dg] * omega^2 betat00 <- betat0 - f_c - n.t <- nextodd((-betat1 - sqrt(betat1^2 - 4 * betat00 * betat2)) / (2 * betat00)) + n.t <- nextodd( + (-betat1 - sqrt(betat1^2 - 4 * betat00 * betat2)) / (2 * betat00) + ) n.t } @@ -196,40 +298,51 @@ stlplus.default <- function(x, t = NULL, n.p, s.window, s.degree = 1, nn <- length(d$y) # check to make sure there is at least one valid value - if(all(is.na(d$y))) { + if (all(is.na(d$y))) { a <- rep(NA, 1:(nn + 2)) } else { cs.ev <- seq(1, length(d$y), by = s.jump) - if(tail(cs.ev, 1) != nn) cs.ev <- c(cs.ev, nn) + if (tail(cs.ev, 1) != nn) { + cs.ev <- c(cs.ev, nn) + } cs.ev <- c(0, cs.ev, nn + 1) # print(m) - a <- .loess_stlplus(y = d$y, span = s.window, degree = s.degree, - weights = d$w, m = cs.ev, noNA = noNA, jump = s.jump, - at = c(1:(nn + 2))) + a <- .loess_stlplus( + y = d$y, + span = s.window, + degree = s.degree, + weights = d$w, + m = cs.ev, + noNA = noNA, + jump = s.jump, + at = c(1:(nn + 2)) + ) } c(a, rep(NA, csLength - nn - 2)) } - if(is.null(l.window)) { + if (is.null(l.window)) { l.window <- nextodd(n.p) } else { l.window <- wincheck(l.window) } - if(is.null(sub.labels)) { + if (is.null(sub.labels)) { sub.labels <- paste("subseries", 1:n.p) } else { - if(length(sub.labels) != n.p) stop("sub.labels must be of length n.p") - if(length(unique(sub.labels)) != n.p) stop("sub.labels must be unique") + if (length(sub.labels) != n.p) { + stop("sub.labels must be of length n.p") + } + if (length(unique(sub.labels)) != n.p) stop("sub.labels must be unique") } tmp <- c(sub.start:n.p, rep(1:n.p, ceiling(n / n.p)))[1:n] sub.labels <- factor(tmp, labels = sub.labels) periodic <- FALSE if (is.character(s.window)) { - if (is.na(pmatch(s.window, "periodic"))) + if (is.na(pmatch(s.window, "periodic"))) { stop("unknown string value for s.window") - else { + } else { periodic <- TRUE s.window <- 10 * n + 1 s.degree <- 0 @@ -250,9 +363,15 @@ stlplus.default <- function(x, t = NULL, n.p, s.window, s.degree = 1, t.window <- wincheck(t.window) } - if(is.null(s.jump) || length(s.jump)==0) s.jump <- ceiling(s.window / 10) - if(is.null(t.jump) || length(t.jump)==0) t.jump <- ceiling(t.window / 10) - if(is.null(l.jump) || length(l.jump)==0) l.jump <- ceiling(l.window / 10) + if (is.null(s.jump) || length(s.jump) == 0) { + s.jump <- ceiling(s.window / 10) + } + if (is.null(t.jump) || length(t.jump) == 0) { + t.jump <- ceiling(t.window / 10) + } + if (is.null(l.jump) || length(l.jump) == 0) { + l.jump <- ceiling(l.window / 10) + } # cat(s.degree, " ", t.degree, " ", l.degree, "\n") # cat(s.window, " ", t.window, " ", l.window, "\n") @@ -268,8 +387,9 @@ stlplus.default <- function(x, t = NULL, n.p, s.window, s.degree = 1, # seasonal each observation belongs to cycleSubIndices <- rep(c(1:n.p), ceiling(n / n.p))[1:n] - if(any(by(Y, list(cycleSubIndices), function(x) all(is.na(x))))) + if (any(by(Y, list(cycleSubIndices), function(x) all(is.na(x))))) { stop("There is at least one subseries for which all values are missing.") + } C <- rep(NA, n + 2 * n.p) @@ -277,17 +397,15 @@ stlplus.default <- function(x, t = NULL, n.p, s.window, s.degree = 1, w <- rep(1, n) - for(o_iter in 1:outer) { - - for(iter in 1:inner) { - + for (o_iter in 1:outer) { + for (iter in 1:inner) { # step 1: detrending... Y.detrended <- Y - trend csLength <- ceiling(n / n.p) + 2 # step 2: smoothing of cycle-subseries - for(i in 1:n.p) { + for (i in 1:n.p) { cycleSub <- Y.detrended[cycleSubIndices == i] subWeights <- w[cycleSubIndices == i] cycleSub.length <- length(cycleSub) @@ -297,16 +415,27 @@ stlplus.default <- function(x, t = NULL, n.p, s.window, s.degree = 1, notEnoughData <- length(cycleSub[!is.na(cycleSub)]) < s.window / 2 # if(notEnoughData || periodic) { - if(periodic) { - C[c(cs1, cycleSubIndices, cs2) == i] <- rep(weighted.mean(cycleSub, - w = w[cycleSubIndices == i], na.rm = TRUE), cycleSub.length + 2) + if (periodic) { + C[c(cs1, cycleSubIndices, cs2) == i] <- rep( + weighted.mean(cycleSub, w = w[cycleSubIndices == i], na.rm = TRUE), + cycleSub.length + 2 + ) } else { cs.ev <- seq(1, cycleSub.length, by = s.jump) - if(tail(cs.ev, 1) != cycleSub.length) cs.ev <- c(cs.ev, cycleSub.length) + if (tail(cs.ev, 1) != cycleSub.length) { + cs.ev <- c(cs.ev, cycleSub.length) + } cs.ev <- c(0, cs.ev, cycleSub.length + 1) - tmps <- .loess_stlplus(y = cycleSub, span = s.window, degree = s.degree, - m = cs.ev, weights = w[cycleSubIndices == i], blend = s.blend, - jump = s.jump, at = c(0:(cycleSub.length + 1))) + tmps <- .loess_stlplus( + y = cycleSub, + span = s.window, + degree = s.degree, + m = cs.ev, + weights = w[cycleSubIndices == i], + blend = s.blend, + jump = s.jump, + at = c(0:(cycleSub.length + 1)) + ) C[c(cs1, cycleSubIndices, cs2) == i] <- tmps # approx(x = cs.ev, y = tmps, xout = c(0:(cycleSub.length + 1)))$y } @@ -317,10 +446,21 @@ stlplus.default <- function(x, t = NULL, n.p, s.window, s.degree = 1, ma3 <- c_ma(C, n.p) l.ev <- seq(1, n, by = l.jump) - if(tail(l.ev, 1) != n) l.ev <- c(l.ev, n) - L <- .loess_stlplus(y = ma3, span = l.window, degree = l.degree, - m = l.ev, weights = w, y_idx = y_idx, noNA = noNA, - blend = l.blend, jump = l.jump, at = c(1:n)) + if (tail(l.ev, 1) != n) { + l.ev <- c(l.ev, n) + } + L <- .loess_stlplus( + y = ma3, + span = l.window, + degree = l.degree, + m = l.ev, + weights = w, + y_idx = y_idx, + noNA = noNA, + blend = l.blend, + jump = l.jump, + at = c(1:n) + ) # L <- predict(loess(ma3 ~ c(1:n), degree = l.degree, # span = l.window / n, family = family), newdata = c(1:n)) @@ -334,10 +474,21 @@ stlplus.default <- function(x, t = NULL, n.p, s.window, s.degree = 1, # Step 6: Trend Smoothing t.ev <- seq(1, n, by = t.jump) - if(tail(t.ev, 1) != n) t.ev <- c(t.ev, n) - trend <- .loess_stlplus(y = D, span = t.window, degree = t.degree, - m = t.ev, weights = w, y_idx = y_idx, noNA = noNA, - blend = t.blend, jump = t.jump, at = c(1:n)) + if (tail(t.ev, 1) != n) { + t.ev <- c(t.ev, n) + } + trend <- .loess_stlplus( + y = D, + span = t.window, + degree = t.degree, + m = t.ev, + weights = w, + y_idx = y_idx, + noNA = noNA, + blend = t.blend, + jump = t.jump, + at = c(1:n) + ) # if(blend$t) { # # TODO: validate blend parameters # trend0 <- .loess_stlplus(y = D, span = t.window, degree = t.degree, @@ -347,12 +498,23 @@ stlplus.default <- function(x, t = NULL, n.p, s.window, s.degree = 1, # trend <- predict(loess(D ~ c(1:length(D)), degree = t.degree, # span = t.window / length(D), family = family), newdata = c(1:length(D))) - if(details) - dtls <- c(dtls, list(list(trend = trend, seasonal = seasonal, - C = C, D = D, L = L, Y.detrended = Y.detrended, weights = w))) + if (details) { + dtls <- c( + dtls, + list(list( + trend = trend, + seasonal = seasonal, + C = C, + D = D, + L = L, + Y.detrended = Y.detrended, + weights = w + )) + ) + } } - if(outer > 1) { + if (outer > 1) { mid1 <- floor(n / 2 + 1) mid2 <- n - mid1 + 1 R <- Y - seasonal - trend @@ -372,78 +534,128 @@ stlplus.default <- function(x, t = NULL, n.p, s.window, s.degree = 1, fc <- NULL fc.number <- 0 fc.res <- NULL - if(!is.null(fc.window)) { + if (!is.null(fc.window)) { fc.number <- length(fc.window) fc.window <- wincheck(fc.window) - if(is.null(fc.degree)) + if (is.null(fc.degree)) { fc.degree <- 1 - if(length(fc.degree) < fc.number) - fc.degree <- c(fc.degree, - rep(fc.degree[length(fc.degree)], fc.number - length(fc.degree))) + } + if (length(fc.degree) < fc.number) { + fc.degree <- c( + fc.degree, + rep(fc.degree[length(fc.degree)], fc.number - length(fc.degree)) + ) + } fc.cumulative <- rep(0, n) # keep sum of all previous fc smoothings degcheck(fc.degree) - if(is.null(fc.name)) + if (is.null(fc.name)) { fc.name <- paste("fc.", fc.window, sep = "") + } - if(is.null(fc.jump)) + if (is.null(fc.jump)) { fc.jump <- ceiling(fc.window / 10) + } - if(is.null(fc.blend)) + if (is.null(fc.blend)) { fc.blend <- rep(t.blend, fc.number) + } - if(length(fc.blend) < fc.number) + if (length(fc.blend) < fc.number) { fc.blend <- c(fc.blend, rep(0, fc.number - length(fc.blend))) + } - if(length(fc.jump) < fc.number) - fc.jump <- c(fc.jump, - rep(fc.jump[length(fc.jump)], fc.number - length(fc.jump))) + if (length(fc.jump) < fc.number) { + fc.jump <- c( + fc.jump, + rep(fc.jump[length(fc.jump)], fc.number - length(fc.jump)) + ) + } fc.res <- data.frame(matrix(nrow = n, ncol = fc.number, data = 0)) - for(ii in 1:fc.number) { + for (ii in 1:fc.number) { fc.ev <- seq(1, n, by = fc.jump[ii]) - if(tail(fc.ev, 1) != n) fc.ev <- c(fc.ev, n) - tmp <- .loess_stlplus(y = Y - seasonal - fc.cumulative, - span = fc.window[ii], degree = fc.degree[ii], - m = fc.ev, weights = w, y_idx = y_idx, noNA = noNA, - blend = fc.blend[ii], jump = fc.jump[ii], at = c(1:n)) + if (tail(fc.ev, 1) != n) { + fc.ev <- c(fc.ev, n) + } + tmp <- .loess_stlplus( + y = Y - seasonal - fc.cumulative, + span = fc.window[ii], + degree = fc.degree[ii], + m = fc.ev, + weights = w, + y_idx = y_idx, + noNA = noNA, + blend = fc.blend[ii], + jump = fc.jump[ii], + at = c(1:n) + ) fc.cumulative <- fc.cumulative + tmp - fc.res[,ii] <- tmp + fc.res[, ii] <- tmp } - if(any(is.null(fc.name)) || length(fc.name) < fc.number) + if (any(is.null(fc.name)) || length(fc.name) < fc.number) { fc.name <- paste("trend", c(1:fc.number)) + } names(fc.res) <- fc.name fc.res$remainder <- Y - seasonal - apply(fc.res, 1, sum) fc <- data.frame( - fc.window = fc.window, fc.degree = fc.degree, - fc.name = fc.name, fc.jump = fc.jump, fc.blend = fc.blend) + fc.window = fc.window, + fc.degree = fc.degree, + fc.name = fc.name, + fc.jump = fc.jump, + fc.blend = fc.blend + ) } # compute remainder R <- Y - seasonal - trend - dat <- data.frame(raw = Y, seasonal = seasonal, trend = trend, - remainder = R, weights = w, sub.labels = sub.labels) + dat <- data.frame( + raw = Y, + seasonal = seasonal, + trend = trend, + remainder = R, + weights = w, + sub.labels = sub.labels + ) pars <- list( - deg = data.frame(s.degree = s.degree, t.degree = t.degree, l.degree = l.degree), - win = data.frame(s.window = s.window, t.window = t.window, l.window = l.window), + deg = data.frame( + s.degree = s.degree, + t.degree = t.degree, + l.degree = l.degree + ), + win = data.frame( + s.window = s.window, + t.window = t.window, + l.window = l.window + ), blend = data.frame(s.blend = s.blend, t.blend = t.blend, l.blend = l.blend), jump = data.frame(s.jump = s.jump, t.jump = t.jump, l.jump = l.jump), - fc.number = fc.number, fc = fc, n.p = n.p, - inner = inner, outer = outer, periodic = periodic) - - res <- list(data = dat, fc = fc.res, pars = pars, time = t, n = nrow(dat), - details = dtls) + fc.number = fc.number, + fc = fc, + n.p = n.p, + inner = inner, + outer = outer, + periodic = periodic + ) + + res <- list( + data = dat, + fc = fc.res, + pars = pars, + time = t, + n = nrow(dat), + details = dtls + ) class(res) <- "stlplus" res } - diff --git a/man-roxygen/ex-stlplus.R b/man-roxygen/ex-stlplus.R index d5b3797..b9f6611 100644 --- a/man-roxygen/ex-stlplus.R +++ b/man-roxygen/ex-stlplus.R @@ -1,6 +1,13 @@ -co2_stl <- stlplus(co2, t = as.vector(time(co2)), n.p = 12, - l.window = 13, t.window = 19, s.window = 35, s.degree = 1, - sub.labels = substr(month.name, 1, 3)) +co2_stl <- stlplus( + co2, + t = as.vector(time(co2)), + n.p = 12, + l.window = 13, + t.window = 19, + s.window = 35, + s.degree = 1, + sub.labels = substr(month.name, 1, 3) +) plot(co2_stl, ylab = "CO2 Concentration (ppm)", xlab = "Time (years)") plot_seasonal(co2_stl) @@ -10,42 +17,81 @@ plot_rembycycle(co2_stl) # post-trend smoothing -co2_stl_pt <- stlplus(co2, t = as.vector(time(co2)), n.p = 12, - l.window = 13, t.window = 19, s.window = 35, s.degree = 1, +co2_stl_pt <- stlplus( + co2, + t = as.vector(time(co2)), + n.p = 12, + l.window = 13, + t.window = 19, + s.window = 35, + s.degree = 1, sub.labels = substr(month.name, 1, 3), - fc.degree = c(1, 2), fc.window = c(201, 35), - fc.name = c("long-term", "so. osc.")) - -plot(co2_stl_pt, scales = list(y = list(relation = "free")), - ylab = "CO2 Concentration (ppm)", xlab = "Time (years)", - aspect = 0.25, type = c("l", "g")) + fc.degree = c(1, 2), + fc.window = c(201, 35), + fc.name = c("long-term", "so. osc.") +) + +plot( + co2_stl_pt, + scales = list(y = list(relation = "free")), + ylab = "CO2 Concentration (ppm)", + xlab = "Time (years)", + aspect = 0.25, + type = c("l", "g") +) # with NAs y <- co2 y[201:224] <- NA -y_stl <- stlplus(y, l.window = 13, t.window = 19, s.window = 35, - s.degree = 1, sub.labels = substr(month.name, 1, 3)) - -plot(y_stl, ylab = "CO2 Concentration (ppm)", xlab = "Time (years)", type = c("l", "g")) +y_stl <- stlplus( + y, + l.window = 13, + t.window = 19, + s.window = 35, + s.degree = 1, + sub.labels = substr(month.name, 1, 3) +) + +plot( + y_stl, + ylab = "CO2 Concentration (ppm)", + xlab = "Time (years)", + type = c("l", "g") +) plot_seasonal(y_stl) plot_trend(y_stl) plot_cycle(y_stl) plot_rembycycle(y_stl) # if you don't want to use a time series object: -y_stl <- stlplus(y, t = as.vector(time(y)), n.p = 12, - l.window = 13, t.window = 19, s.window = 35, s.degree = 1, - sub.labels = substr(month.name, 1, 3)) +y_stl <- stlplus( + y, + t = as.vector(time(y)), + n.p = 12, + l.window = 13, + t.window = 19, + s.window = 35, + s.degree = 1, + sub.labels = substr(month.name, 1, 3) +) # with an outlier y2 <- co2 y2[200] <- 300 -y2_stl <- stlplus(y2, t = as.vector(time(y2)), n.p = 12, - l.window = 13, t.window = 19, s.window = 35, s.degree = 1, - sub.labels = substr(month.name, 1, 3), outer = 10) +y2_stl <- stlplus( + y2, + t = as.vector(time(y2)), + n.p = 12, + l.window = 13, + t.window = 19, + s.window = 35, + s.degree = 1, + sub.labels = substr(month.name, 1, 3), + outer = 10 +) plot(y2_stl, ylab = "CO2 Concentration (ppm)", xlab = "Time (years)") plot_seasonal(y2_stl) @@ -55,9 +101,16 @@ plot_rembycycle(y2_stl) # compare to R's stl -x1 <- stlplus(co2, t = as.vector(time(co2)), n.p = 12, - l.window = 13, t.window = 19, s.window = 11, s.degree = 1, - sub.labels = substr(month.name, 1, 3)) +x1 <- stlplus( + co2, + t = as.vector(time(co2)), + n.p = 12, + l.window = 13, + t.window = 19, + s.window = 11, + s.degree = 1, + sub.labels = substr(month.name, 1, 3) +) x2 <- stl(co2, l.window = 13, t.window = 19, s.window = 11, s.degree = 1) @@ -65,12 +118,29 @@ x2 <- stl(co2, l.window = 13, t.window = 19, s.window = 11, s.degree = 1) plot(seasonal(x1) - seasonal(x2)) # but not if all jump parameters are 1 -x1 <- stlplus(co2, t = as.vector(time(co2)), n.p = 12, - l.window = 13, t.window = 19, s.window = 11, s.degree = 1, +x1 <- stlplus( + co2, + t = as.vector(time(co2)), + n.p = 12, + l.window = 13, + t.window = 19, + s.window = 11, + s.degree = 1, sub.labels = substr(month.name, 1, 3), - s.jump = 1, t.jump = 1, l.jump = 1) - -x2 <- stl(co2, l.window = 13, t.window = 19, s.window = 11, s.degree = 1, - s.jump = 1, t.jump = 1, l.jump = 1) + s.jump = 1, + t.jump = 1, + l.jump = 1 +) + +x2 <- stl( + co2, + l.window = 13, + t.window = 19, + s.window = 11, + s.degree = 1, + s.jump = 1, + t.jump = 1, + l.jump = 1 +) plot(seasonal(x1) - seasonal(x2)) diff --git a/tests/testthat/tests.R b/tests/testthat/tests.R index 0e9b8a4..035145d 100755 --- a/tests/testthat/tests.R +++ b/tests/testthat/tests.R @@ -1,4 +1,3 @@ - test_that("direct loess matches R's loess", { n <- 100 x <- 1:n @@ -24,21 +23,55 @@ test_that("direct loess matches R's loess", { }) test_that("stlplus matches stl", { - x1 <- stlplus(co2, t = as.vector(time(co2)), n.p = 12, - l.window = 13, t.window = 19, s.window = 11, s.degree = 1, - s.jump = 1, t.jump = 1, l.jump = 1) - x2 <- stl(co2, l.window = 13, t.window = 19, s.window = 11, s.degree = 1, - s.jump = 1, t.jump = 1, l.jump = 1) + x1 <- stlplus( + co2, + t = as.vector(time(co2)), + n.p = 12, + l.window = 13, + t.window = 19, + s.window = 11, + s.degree = 1, + s.jump = 1, + t.jump = 1, + l.jump = 1 + ) + x2 <- stl( + co2, + l.window = 13, + t.window = 19, + s.window = 11, + s.degree = 1, + s.jump = 1, + t.jump = 1, + l.jump = 1 + ) expect_true(mean(abs(seasonal(x1) - seasonal(x2))) < 1.0e-12) - expect_true(mean(abs(trend(x1) - trend(x2))) < 1.0e-12) + expect_true(mean(abs(trend(x1) - trend(x2))) < 1.0e-12) - x1 <- stlplus(co2, t = as.vector(time(co2)), n.p = 12, - l.window = 13, t.window = 39, s.window = 101, s.degree = 1, - s.jump = 1, t.jump = 1, l.jump = 1) - x2 <- stl(co2, l.window = 13, t.window = 39, s.window = 101, s.degree = 1, - s.jump = 1, t.jump = 1, l.jump = 1) + x1 <- stlplus( + co2, + t = as.vector(time(co2)), + n.p = 12, + l.window = 13, + t.window = 39, + s.window = 101, + s.degree = 1, + s.jump = 1, + t.jump = 1, + l.jump = 1 + ) + x2 <- stl( + co2, + l.window = 13, + t.window = 39, + s.window = 101, + s.degree = 1, + s.jump = 1, + t.jump = 1, + l.jump = 1 + ) expect_true(mean(abs(seasonal(x1) - seasonal(x2))) < 1.0e-12) - expect_true(mean(abs(trend(x1) - trend(x2))) < 1.0e-12) + expect_true(mean(abs(trend(x1) - trend(x2))) < 1.0e-12) }) # compare to loess with missing values From 24de8842fbf2a6f490690cfb88cffcee8405f686 Mon Sep 17 00:00:00 2001 From: Ryan Hafen Date: Mon, 6 Apr 2026 11:53:57 -0700 Subject: [PATCH 2/5] Cleanup --- .Rbuildignore | 2 ++ .github/workflows/R-CMD-check.yaml | 50 ++++++++++++++++++++++++++++++ .gitignore | 6 ++++ DESCRIPTION | 9 +++--- NAMESPACE | 2 +- R/RcppExports.R | 18 ++++++++--- R/loess_stl.R | 13 ++++---- R/plots.R | 19 ++++++------ R/stlplus.R | 28 ++++++++--------- README.md | 2 +- man/plot_seasonal.Rd | 7 +++-- src/loess.cpp | 7 ----- tests/testthat/tests.R | 41 +++++++++++++++++++++--- 13 files changed, 149 insertions(+), 55 deletions(-) create mode 100644 .github/workflows/R-CMD-check.yaml diff --git a/.Rbuildignore b/.Rbuildignore index a5f6d58..e0a855b 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -1,8 +1,10 @@ man-roxygen \.git +^\.github$ ^\.travis\.yml$ \.DS_Store ..Rcheck +^src/.*\.(o|so)$ LICENSE\.md NEWS\.md ^_ignore$ diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml new file mode 100644 index 0000000..0e13b0c --- /dev/null +++ b/.github/workflows/R-CMD-check.yaml @@ -0,0 +1,50 @@ +name: R-CMD-check + +on: + push: + branches: + - main + - master + - dev + pull_request: + branches: + - main + - master + workflow_dispatch: + +jobs: + R-CMD-check: + runs-on: ${{ matrix.config.os }} + name: ${{ matrix.config.os }} (R ${{ matrix.config.r }}) + + strategy: + fail-fast: false + matrix: + config: + - os: ubuntu-latest + r: release + - os: ubuntu-latest + r: devel + - os: windows-latest + r: release + + env: + GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} + R_KEEP_PKG_SOURCE: yes + + steps: + - uses: actions/checkout@v4 + + - uses: r-lib/actions/setup-r@v2 + with: + r-version: ${{ matrix.config.r }} + use-public-rspm: true + + - uses: r-lib/actions/setup-r-dependencies@v2 + with: + extra-packages: any::rcmdcheck + needs: check + + - uses: r-lib/actions/check-r-package@v2 + with: + args: 'c("--no-manual", "--as-cran")' \ No newline at end of file diff --git a/.gitignore b/.gitignore index 25a3194..89c8e5c 100644 --- a/.gitignore +++ b/.gitignore @@ -5,3 +5,9 @@ _ignore # Example code in package build process *-Ex.R + +# Local compiled artifacts +src/*.o +src/*.so +src/symbols.rds +src/.DS_Store diff --git a/DESCRIPTION b/DESCRIPTION index 21843dd..f8a3bb6 100755 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,9 +1,11 @@ Package: stlplus Type: Package Title: Enhanced Seasonal Decomposition of Time Series by Loess -Version: 0.5.1 -Date: 2016-01-05 -Authors@R: person("Ryan", "Hafen", email = "rhafen@gmail.com", role = c("aut", "cre")) +Version: 0.5.2 +Date: 2026-04-06 +Author: Ryan Hafen [aut, cre] +Authors@R: person("Ryan", "Hafen", email = "rhafen@gmail.com", role = c("aut", "cre"), + comment = c(ORCID = "0000-0002-5516-8367")) Maintainer: Ryan Hafen Description: Decompose a time series into seasonal, trend, and remainder components using an implementation of Seasonal Decomposition of Time @@ -13,7 +15,6 @@ Description: Decompose a time series into seasonal, trend, and remainder parameter choices, frequency component smoothing beyond the seasonal and trend components, and some basic plot methods for diagnostics. LazyLoad: yes -LazyData: yes NeedsCompilation: yes Imports: lattice, diff --git a/NAMESPACE b/NAMESPACE index 102a592..d01eaf2 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -43,4 +43,4 @@ importFrom(utils,head) importFrom(utils,stack) importFrom(utils,tail) importFrom(yaImpute,ann) -useDynLib(stlplus) +useDynLib(stlplus, .registration = TRUE) diff --git a/R/RcppExports.R b/R/RcppExports.R index e08fa72..ea3b4e1 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -2,14 +2,24 @@ # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 c_interp <- function(m, fits, slopes, at) { - .Call('_stlplus_c_interp', PACKAGE = 'stlplus', m, fits, slopes, at) + .Call("_stlplus_c_interp", PACKAGE = "stlplus", m, fits, slopes, at) } c_loess <- function(xx, yy, degree, span, ww, m, l_idx, max_dist) { - .Call('_stlplus_c_loess', PACKAGE = 'stlplus', xx, yy, degree, span, ww, m, l_idx, max_dist) + .Call( + "_stlplus_c_loess", + PACKAGE = "stlplus", + xx, + yy, + degree, + span, + ww, + m, + l_idx, + max_dist + ) } c_ma <- function(x, n_p) { - .Call('_stlplus_c_ma', PACKAGE = 'stlplus', x, n_p) + .Call("_stlplus_c_ma", PACKAGE = "stlplus", x, n_p) } - diff --git a/R/loess_stl.R b/R/loess_stl.R index 3f9db01..8e19f90 100755 --- a/R/loess_stl.R +++ b/R/loess_stl.R @@ -1,18 +1,17 @@ #' @importFrom yaImpute ann -#' @importFrom Rcpp sourceCpp -#' @useDynLib stlplus +#' @useDynLib stlplus, .registration = TRUE .loess_stlplus <- function( x = NULL, y, span, degree, weights = NULL, - m = c(1:length(y)), + m = c(seq_along(y)), y_idx = !is.na(y), noNA = all(y_idx), blend = 0, jump = ceiling(span / 10), - at = c(1:length(y)) + at = c(seq_along(y)) ) { nextodd <- function(x) { x <- round(x) @@ -22,7 +21,7 @@ n <- length(y[y_idx]) if (is.null(x)) { - x <- c(1:length(y)) + x <- c(seq_along(y)) } if (is.null(weights)) { @@ -126,11 +125,11 @@ # left_interp <- at[bl_idx_interp] # right_interp <- at[br_idx_interp] l_idx2 <- l_idx[bl_idx | br_idx] - r_idx2 <- r_idx[bl_idx | br_idx] + # r_idx2 <- r_idx[bl_idx | br_idx] max_dist2 <- max_dist[bl_idx | br_idx] m2 <- c(left, right) - n_m2 <- length(m2) + # n_m2 <- length(m2) # speed this up later by only getting the loess smooth at the tails. # right now, a lot of unnecessary calculation is done at the interior diff --git a/R/plots.R b/R/plots.R index 417f88d..ee96d55 100755 --- a/R/plots.R +++ b/R/plots.R @@ -149,11 +149,11 @@ plot_cycle <- function( seas <- seasonal(x) t <- time.stlplus(x) - cycleSubIndices <- x$data$sub.labels + cycle_sub_idx <- x$data$sub.labels if (x$pars$periodic) { p <- lattice::xyplot( - seas ~ t | cycleSubIndices, + seas ~ t | cycle_sub_idx, layout = layout, type = "l", xlab = xlab, @@ -162,7 +162,7 @@ plot_cycle <- function( ) } else { p <- lattice::xyplot( - seas ~ t | cycleSubIndices, + seas ~ t | cycle_sub_idx, layout = layout, type = "l", panel = panel, @@ -178,8 +178,9 @@ plot_cycle <- function( #' Seasonal Diagnostic Plot for an stlplus Object #' #' Plots each cycle-subseries of the detrended data (or equivalently, -#' seasonal plus remainder), with the mean subtracted. The fitted seasonal -#' component is plotted as a line through the points. +#' seasonal plus remainder), centered by the mean seasonal component for each +#' subseries. The fitted seasonal component is plotted as a line through the +#' points. #' #' @param x object of class \code{"stlplus"}. #' @param col,lwd,xlab,ylab,\ldots parameters to be passed to \code{xyplot()}. @@ -194,7 +195,7 @@ plot_seasonal <- function( col = c("darkgray", "black"), lwd = 2, xlab = "Time", - ylab = "Centered Seasonal + Remainder", + ylab = "Seasonal + Remainder\n(centered by seasonal mean)", ... ) { dat <- by( @@ -206,7 +207,7 @@ plot_seasonal <- function( ), list(x$data$sub.labels), function(dd) { - mn <- mean(dd$v1, na.rm = TRUE) + mn <- mean(dd$v2, na.rm = TRUE) data.frame(a = dd$v1 - mn, b = dd$v2 - mn, t = dd$t, sub.labels = dd$v3) } ) @@ -329,11 +330,11 @@ plot_trend <- function( x = time.stlplus(x), y = predict( loess( - remainder(x) ~ c(1:length(remainder(x))), + remainder(x) ~ c(seq_along(remainder(x))), span = span, weights = x$data$weights ), - newdata = c(1:length(remainder(x))) + newdata = c(seq_along(remainder(x))) ), type = "l", pan = "Remainder" diff --git a/R/stlplus.R b/R/stlplus.R index 173b261..d4141bf 100755 --- a/R/stlplus.R +++ b/R/stlplus.R @@ -406,38 +406,38 @@ stlplus.default <- function( # step 2: smoothing of cycle-subseries for (i in 1:n.p) { - cycleSub <- Y.detrended[cycleSubIndices == i] - subWeights <- w[cycleSubIndices == i] - cycleSub.length <- length(cycleSub) + cycle_sub <- Y.detrended[cycleSubIndices == i] + # sub_weights <- w[cycleSubIndices == i] + cycle_sub_length <- length(cycle_sub) cs1 <- head(cycleSubIndices, n.p) cs2 <- tail(cycleSubIndices, n.p) - notEnoughData <- length(cycleSub[!is.na(cycleSub)]) < s.window / 2 - # if(notEnoughData || periodic) { + # not_enough_dat <- length(cycle_sub[!is.na(cycle_sub)]) < s.window / 2 + # if(not_enough_dat || periodic) { if (periodic) { C[c(cs1, cycleSubIndices, cs2) == i] <- rep( - weighted.mean(cycleSub, w = w[cycleSubIndices == i], na.rm = TRUE), - cycleSub.length + 2 + weighted.mean(cycle_sub, w = w[cycleSubIndices == i], na.rm = TRUE), + cycle_sub_length + 2 ) } else { - cs.ev <- seq(1, cycleSub.length, by = s.jump) - if (tail(cs.ev, 1) != cycleSub.length) { - cs.ev <- c(cs.ev, cycleSub.length) + cs.ev <- seq(1, cycle_sub_length, by = s.jump) + if (tail(cs.ev, 1) != cycle_sub_length) { + cs.ev <- c(cs.ev, cycle_sub_length) } - cs.ev <- c(0, cs.ev, cycleSub.length + 1) + cs.ev <- c(0, cs.ev, cycle_sub_length + 1) tmps <- .loess_stlplus( - y = cycleSub, + y = cycle_sub, span = s.window, degree = s.degree, m = cs.ev, weights = w[cycleSubIndices == i], blend = s.blend, jump = s.jump, - at = c(0:(cycleSub.length + 1)) + at = c(0:(cycle_sub_length + 1)) ) C[c(cs1, cycleSubIndices, cs2) == i] <- tmps - # approx(x = cs.ev, y = tmps, xout = c(0:(cycleSub.length + 1)))$y + # approx(x = cs.ev, y = tmps, xout = c(0:(cycle_sub_length + 1)))$y } } diff --git a/README.md b/README.md index 3e6e4f9..96f281a 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ ## stlplus -[![Build Status](https://travis-ci.org/hafen/stlplus.svg?branch=master)](https://travis-ci.org/hafen/stlplus) +[![R-CMD-check](https://github.com/hafen/stlplus/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/hafen/stlplus/actions/workflows/R-CMD-check.yaml) [![CRAN](http://www.r-pkg.org/badges/version/stlplus)](https://cran.r-project.org/web/packages/stlplus/index.html) ![](http://cranlogs.r-pkg.org/badges/stlplus) diff --git a/man/plot_seasonal.Rd b/man/plot_seasonal.Rd index fa895e0..38ac163 100644 --- a/man/plot_seasonal.Rd +++ b/man/plot_seasonal.Rd @@ -5,7 +5,7 @@ \title{Seasonal Diagnostic Plot for an stlplus Object} \usage{ plot_seasonal(x, col = c("darkgray", "black"), lwd = 2, xlab = "Time", - ylab = "Centered Seasonal + Remainder", ...) + ylab = "Seasonal + Remainder\n(centered by seasonal mean)", ...) } \arguments{ \item{x}{object of class \code{"stlplus"}.} @@ -17,8 +17,9 @@ object of class \code{"trellis"}. } \description{ Plots each cycle-subseries of the detrended data (or equivalently, -seasonal plus remainder), with the mean subtracted. The fitted seasonal -component is plotted as a line through the points. +seasonal plus remainder), centered by the mean seasonal component for each +subseries. The fitted seasonal component is plotted as a line through the +points. } \details{ Helps decide how much of the variation in the data other than the trend should go into the seasonal component, and how much in the remainder. diff --git a/src/loess.cpp b/src/loess.cpp index 61c5f9e..7933a6d 100755 --- a/src/loess.cpp +++ b/src/loess.cpp @@ -12,7 +12,6 @@ List c_loess( IntegerVector l_idx, // index of left starting points NumericVector max_dist // distance between nn bounds for each point ) { - int span2, span3, offset; int i, j; double r, tmp1, tmp2; @@ -31,16 +30,10 @@ List c_loess( // variables for storing determinant intermediate values double a, b, c, d, e, a1, b1, c1, a2, b2, c2, det; - span3 = span; if(span > n) { span = n; } - span2 = (span - 1) / 2; - - // want to start storing results at index 0, corresponding to the lowest m - offset = m[0]; - // loop through all values of m for(i = 0; i < n_m; i++) { a = 0.0; diff --git a/tests/testthat/tests.R b/tests/testthat/tests.R index 035145d..ac37d79 100755 --- a/tests/testthat/tests.R +++ b/tests/testthat/tests.R @@ -31,6 +31,7 @@ test_that("stlplus matches stl", { t.window = 19, s.window = 11, s.degree = 1, + outer = 1, s.jump = 1, t.jump = 1, l.jump = 1 @@ -43,10 +44,11 @@ test_that("stlplus matches stl", { s.degree = 1, s.jump = 1, t.jump = 1, - l.jump = 1 + l.jump = 1, + robust = FALSE ) - expect_true(mean(abs(seasonal(x1) - seasonal(x2))) < 1.0e-12) - expect_true(mean(abs(trend(x1) - trend(x2))) < 1.0e-12) + expect_true(mean(abs(seasonal(x1) - seasonal(x2))) < 1.0e-10) + expect_true(mean(abs(trend(x1) - trend(x2))) < 1.0e-10) x1 <- stlplus( co2, @@ -56,6 +58,7 @@ test_that("stlplus matches stl", { t.window = 39, s.window = 101, s.degree = 1, + outer = 1, s.jump = 1, t.jump = 1, l.jump = 1 @@ -68,10 +71,38 @@ test_that("stlplus matches stl", { s.degree = 1, s.jump = 1, t.jump = 1, + l.jump = 1, + robust = FALSE + ) + expect_true(mean(abs(seasonal(x1) - seasonal(x2))) < 1.0e-10) + expect_true(mean(abs(trend(x1) - trend(x2))) < 1.0e-10) +}) + +test_that("periodic stlplus matches stl", { + x1 <- stlplus( + co2, + t = as.vector(time(co2)), + n.p = 12, + l.window = 13, + t.window = 19, + s.window = "periodic", + outer = 1, + s.jump = 1, + t.jump = 1, l.jump = 1 ) - expect_true(mean(abs(seasonal(x1) - seasonal(x2))) < 1.0e-12) - expect_true(mean(abs(trend(x1) - trend(x2))) < 1.0e-12) + x2 <- stl( + co2, + l.window = 13, + t.window = 19, + s.window = "periodic", + s.jump = 1, + t.jump = 1, + l.jump = 1, + robust = FALSE + ) + expect_true(mean(abs(seasonal(x1) - seasonal(x2))) < 1.0e-6) + expect_true(mean(abs(trend(x1) - trend(x2))) < 1.0e-6) }) # compare to loess with missing values From 9410338688f5a6cb74abd664108e5a4ad165d1b0 Mon Sep 17 00:00:00 2001 From: Ryan Hafen Date: Mon, 13 Apr 2026 15:45:16 -0700 Subject: [PATCH 3/5] Updates to bring up to speed with last decade of R updates --- DESCRIPTION | 2 +- NAMESPACE | 2 +- R/RcppExports.R | 18 +--- R/stlplus.R | 1 + man/accessors.Rd | 25 +++-- man/plot.stlplus.Rd | 15 ++- man/plot_cycle.Rd | 17 +++- man/plot_rembycycle.Rd | 12 ++- man/plot_seasonal.Rd | 11 ++- man/plot_trend.Rd | 22 +++-- man/stlplus.Rd | 216 +++++++++++++++++++++++++++++++---------- 11 files changed, 239 insertions(+), 102 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index f8a3bb6..edecfc4 100755 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -29,4 +29,4 @@ URL: https://github.com/hafen/stlplus BugReports: https://github.com/hafen/stlplus/issues Note: This is experimental software distributed free of charge and comes with no warranty! Exercise professional caution. -RoxygenNote: 5.0.1 +RoxygenNote: 7.3.3 diff --git a/NAMESPACE b/NAMESPACE index d01eaf2..2261cc3 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -26,7 +26,7 @@ export(remainder) export(seasonal) export(stlplus) export(trend) -importFrom(Rcpp,sourceCpp) +importFrom(Rcpp,evalCpp) importFrom(lattice,panel.abline) importFrom(lattice,panel.loess) importFrom(lattice,panel.segments) diff --git a/R/RcppExports.R b/R/RcppExports.R index ea3b4e1..c6615c4 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -2,24 +2,14 @@ # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 c_interp <- function(m, fits, slopes, at) { - .Call("_stlplus_c_interp", PACKAGE = "stlplus", m, fits, slopes, at) + .Call(`_stlplus_c_interp`, m, fits, slopes, at) } c_loess <- function(xx, yy, degree, span, ww, m, l_idx, max_dist) { - .Call( - "_stlplus_c_loess", - PACKAGE = "stlplus", - xx, - yy, - degree, - span, - ww, - m, - l_idx, - max_dist - ) + .Call(`_stlplus_c_loess`, xx, yy, degree, span, ww, m, l_idx, max_dist) } c_ma <- function(x, n_p) { - .Call("_stlplus_c_ma", PACKAGE = "stlplus", x, n_p) + .Call(`_stlplus_c_ma`, x, n_p) } + diff --git a/R/stlplus.R b/R/stlplus.R index d4141bf..1795913 100755 --- a/R/stlplus.R +++ b/R/stlplus.R @@ -41,6 +41,7 @@ #' @note This is a complete re-implementation of the STL algorithm, with the loess part in C and the rest in R. Moving a lot of the code to R makes it easier to experiment with the method at a very minimal speed cost. Recoding in C instead of using R's built-in loess results in better performance, especially for larger series. #' @seealso \code{\link{plot.stlplus}} for plotting the components, \code{\link{getraw}}, \code{\link{seasonal}}, \code{\link{trend}}, \code{\link{remainder}} for accessing the components. #' @example man-roxygen/ex-stlplus.R +#' @importFrom Rcpp evalCpp #' @importFrom stats frequency loess median predict quantile weighted.mean time #' @importFrom utils head stack tail #' @export diff --git a/man/accessors.Rd b/man/accessors.Rd index 7961e65..9d6c6bc 100644 --- a/man/accessors.Rd +++ b/man/accessors.Rd @@ -1,23 +1,23 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/accessors.R \name{seasonal} -\alias{fc} -\alias{fitted.stl} -\alias{fitted.stlplus} -\alias{getraw} -\alias{predict.stl} -\alias{predict.stlplus} +\alias{seasonal} +\alias{trend} \alias{remainder} -\alias{remainder.stl} +\alias{getraw} \alias{remainder.stlplus} -\alias{seasonal} -\alias{seasonal.stl} +\alias{fitted.stlplus} +\alias{predict.stlplus} \alias{seasonal.stlplus} -\alias{time.stl} +\alias{trend.stlplus} +\alias{fc} \alias{time.stlplus} -\alias{trend} +\alias{remainder.stl} +\alias{seasonal.stl} \alias{trend.stl} -\alias{trend.stlplus} +\alias{time.stl} +\alias{predict.stl} +\alias{fitted.stl} \title{Accessor functions for elements of an stl and stlplus object} \usage{ seasonal(object) @@ -82,4 +82,3 @@ R. B. Cleveland, W. S. Cleveland, J.E. McRae, and I. Terpenning (1990) STL: A \seealso{ \code{\link{stlplus}} } - diff --git a/man/plot.stlplus.Rd b/man/plot.stlplus.Rd index 5c9a679..3af0940 100644 --- a/man/plot.stlplus.Rd +++ b/man/plot.stlplus.Rd @@ -4,9 +4,17 @@ \alias{plot.stlplus} \title{Lattice plot of the raw, seasonal, trend, and remainder components} \usage{ -\method{plot}{stlplus}(x, scales = list(y = list(relation = "sliced")), - type = "l", as.table = TRUE, strip = FALSE, strip.left = TRUE, - between = list(y = 0.5), layout = NULL, ...) +\method{plot}{stlplus}( + x, + scales = list(y = list(relation = "sliced")), + type = "l", + as.table = TRUE, + strip = FALSE, + strip.left = TRUE, + between = list(y = 0.5), + layout = NULL, + ... +) } \arguments{ \item{x}{object of class \code{"stlplus"}.} @@ -25,4 +33,3 @@ R. B. Cleveland, W. S. Cleveland, J. E. McRae, and I. Terpenning (1990) STL: A \seealso{ \code{\link{stlplus}} } - diff --git a/man/plot_cycle.Rd b/man/plot_cycle.Rd index 3f72941..318cb2f 100644 --- a/man/plot_cycle.Rd +++ b/man/plot_cycle.Rd @@ -4,10 +4,18 @@ \alias{plot_cycle} \title{Cycle-Subseries Plot for an stlplus Object} \usage{ -plot_cycle(x, layout = c(x$pars$n.p, 1), col = "#0080ff", xlab = "Time", - ylab = "Seasonal", panel = function(x, y, ...) { - lattice::panel.segments(x, rep(.midmean(y), length(x)), x, y, col = col) }, - ...) +plot_cycle( + x, + layout = c(x$pars$n.p, 1), + col = "#0080ff", + xlab = "Time", + ylab = "Seasonal", + panel = function(x, y, ...) { + lattice::panel.segments(x, rep(.midmean(y), + length(x)), x, y, col = col) + }, + ... +) } \arguments{ \item{x}{object of class \code{"stlplus"}.} @@ -27,4 +35,3 @@ R. B. Cleveland, W. S. Cleveland, J. E. McRae, and I. Terpenning (1990) STL: A \seealso{ \code{\link{stlplus}} } - diff --git a/man/plot_rembycycle.Rd b/man/plot_rembycycle.Rd index cb8babc..e68a0e9 100644 --- a/man/plot_rembycycle.Rd +++ b/man/plot_rembycycle.Rd @@ -4,8 +4,15 @@ \alias{plot_rembycycle} \title{Plot of Remainder Component by Cycle-Subseries} \usage{ -plot_rembycycle(x, col = "darkgray", locol = "black", lolwd = 2, - xlab = "Time", ylab = "Remainder", ...) +plot_rembycycle( + x, + col = "darkgray", + locol = "black", + lolwd = 2, + xlab = "Time", + ylab = "Remainder", + ... +) } \arguments{ \item{x}{object of class \code{"stlplus"}.} @@ -24,4 +31,3 @@ R. B. Cleveland, W. S. Cleveland, J. E. McRae, and I. Terpenning (1990) STL: A \seealso{ \code{\link{stlplus}} } - diff --git a/man/plot_seasonal.Rd b/man/plot_seasonal.Rd index 38ac163..d3f3b55 100644 --- a/man/plot_seasonal.Rd +++ b/man/plot_seasonal.Rd @@ -4,8 +4,14 @@ \alias{plot_seasonal} \title{Seasonal Diagnostic Plot for an stlplus Object} \usage{ -plot_seasonal(x, col = c("darkgray", "black"), lwd = 2, xlab = "Time", - ylab = "Seasonal + Remainder\n(centered by seasonal mean)", ...) +plot_seasonal( + x, + col = c("darkgray", "black"), + lwd = 2, + xlab = "Time", + ylab = "Seasonal + Remainder\\n(centered by seasonal mean)", + ... +) } \arguments{ \item{x}{object of class \code{"stlplus"}.} @@ -30,4 +36,3 @@ R. B. Cleveland, W. S. Cleveland, J. E. McRae, and I. Terpenning (1990) STL: A \seealso{ \code{\link{stlplus}} } - diff --git a/man/plot_trend.Rd b/man/plot_trend.Rd index 1563619..32343b6 100644 --- a/man/plot_trend.Rd +++ b/man/plot_trend.Rd @@ -4,11 +4,22 @@ \alias{plot_trend} \title{Trend Diagnostic Plot for an stlplus Object} \usage{ -plot_trend(x, xlab = "Time", ylab = "Trend", span = 0.3, type = c("p", - "l"), scales = list(y = list(relation = "free")), lwd = c(1, 1), - col = c("darkgray", "black", "darkgray"), layout = c(1, 2), - between = list(y = 0.5), strip = FALSE, strip.left = TRUE, - as.table = TRUE, ...) +plot_trend( + x, + xlab = "Time", + ylab = "Trend", + span = 0.3, + type = c("p", "l"), + scales = list(y = list(relation = "free")), + lwd = c(1, 1), + col = c("darkgray", "black", "darkgray"), + layout = c(1, 2), + between = list(y = 0.5), + strip = FALSE, + strip.left = TRUE, + as.table = TRUE, + ... +) } \arguments{ \item{x}{object of class \code{"stlplus"}.} @@ -32,4 +43,3 @@ R. B. Cleveland, W. S. Cleveland, J. E. McRae, and I. Terpenning (1990) STL: A \seealso{ \code{\link{stlplus}} } - diff --git a/man/stlplus.Rd b/man/stlplus.Rd index 1948521..aef0816 100644 --- a/man/stlplus.Rd +++ b/man/stlplus.Rd @@ -6,24 +6,67 @@ \alias{stlplus.zoo} \title{Seasonal Decomposition of Time Series by Loess} \usage{ -stlplus(x, t = NULL, n.p, s.window, s.degree = 1, t.window = NULL, - t.degree = 1, fc.window = NULL, fc.degree = NULL, fc.name = NULL, - l.window = NULL, l.degree = t.degree, s.jump = ceiling(s.window/10), - t.jump = ceiling(t.window/10), l.jump = ceiling(l.window/10), - fc.jump = NULL, critfreq = 0.05, s.blend = 0, t.blend = 0, - l.blend = t.blend, fc.blend = NULL, inner = 2, outer = 1, - sub.labels = NULL, sub.start = 1, zero.weight = 0.000001, - details = FALSE, ...) - -\method{stlplus}{ts}(x, t = as.numeric(stats::time(x)), n.p = frequency(x), - s.window, s.degree = 1, t.window = NULL, t.degree = 1, - fc.window = NULL, fc.degree = NULL, fc.name = NULL, l.window = NULL, - l.degree = t.degree, s.jump = ceiling(s.window/10), - t.jump = ceiling(t.window/10), l.jump = ceiling(l.window/10), - fc.jump = NULL, critfreq = 0.05, s.blend = 0, t.blend = 0, - l.blend = t.blend, fc.blend = NULL, inner = 2, outer = 1, - sub.labels = NULL, sub.start = 1, zero.weight = 0.000001, - details = FALSE, ...) +stlplus( + x, + t = NULL, + n.p, + s.window, + s.degree = 1, + t.window = NULL, + t.degree = 1, + fc.window = NULL, + fc.degree = NULL, + fc.name = NULL, + l.window = NULL, + l.degree = t.degree, + s.jump = ceiling(s.window/10), + t.jump = ceiling(t.window/10), + l.jump = ceiling(l.window/10), + fc.jump = NULL, + critfreq = 0.05, + s.blend = 0, + t.blend = 0, + l.blend = t.blend, + fc.blend = NULL, + inner = 2, + outer = 1, + sub.labels = NULL, + sub.start = 1, + zero.weight = 0.000001, + details = FALSE, + ... +) + +\method{stlplus}{ts}( + x, + t = as.numeric(stats::time(x)), + n.p = frequency(x), + s.window, + s.degree = 1, + t.window = NULL, + t.degree = 1, + fc.window = NULL, + fc.degree = NULL, + fc.name = NULL, + l.window = NULL, + l.degree = t.degree, + s.jump = ceiling(s.window/10), + t.jump = ceiling(t.window/10), + l.jump = ceiling(l.window/10), + fc.jump = NULL, + critfreq = 0.05, + s.blend = 0, + t.blend = 0, + l.blend = t.blend, + fc.blend = NULL, + inner = 2, + outer = 1, + sub.labels = NULL, + sub.start = 1, + zero.weight = 0.000001, + details = FALSE, + ... +) \method{stlplus}{zoo}(...) } @@ -94,9 +137,16 @@ Cycle-subseries labels are useful for plotting and can be specified through the This is a complete re-implementation of the STL algorithm, with the loess part in C and the rest in R. Moving a lot of the code to R makes it easier to experiment with the method at a very minimal speed cost. Recoding in C instead of using R's built-in loess results in better performance, especially for larger series. } \examples{ -co2_stl <- stlplus(co2, t = as.vector(time(co2)), n.p = 12, - l.window = 13, t.window = 19, s.window = 35, s.degree = 1, - sub.labels = substr(month.name, 1, 3)) +co2_stl <- stlplus( + co2, + t = as.vector(time(co2)), + n.p = 12, + l.window = 13, + t.window = 19, + s.window = 35, + s.degree = 1, + sub.labels = substr(month.name, 1, 3) +) plot(co2_stl, ylab = "CO2 Concentration (ppm)", xlab = "Time (years)") plot_seasonal(co2_stl) @@ -106,42 +156,81 @@ plot_rembycycle(co2_stl) # post-trend smoothing -co2_stl_pt <- stlplus(co2, t = as.vector(time(co2)), n.p = 12, - l.window = 13, t.window = 19, s.window = 35, s.degree = 1, +co2_stl_pt <- stlplus( + co2, + t = as.vector(time(co2)), + n.p = 12, + l.window = 13, + t.window = 19, + s.window = 35, + s.degree = 1, sub.labels = substr(month.name, 1, 3), - fc.degree = c(1, 2), fc.window = c(201, 35), - fc.name = c("long-term", "so. osc.")) - -plot(co2_stl_pt, scales = list(y = list(relation = "free")), - ylab = "CO2 Concentration (ppm)", xlab = "Time (years)", - aspect = 0.25, type = c("l", "g")) + fc.degree = c(1, 2), + fc.window = c(201, 35), + fc.name = c("long-term", "so. osc.") +) + +plot( + co2_stl_pt, + scales = list(y = list(relation = "free")), + ylab = "CO2 Concentration (ppm)", + xlab = "Time (years)", + aspect = 0.25, + type = c("l", "g") +) # with NAs y <- co2 y[201:224] <- NA -y_stl <- stlplus(y, l.window = 13, t.window = 19, s.window = 35, - s.degree = 1, sub.labels = substr(month.name, 1, 3)) - -plot(y_stl, ylab = "CO2 Concentration (ppm)", xlab = "Time (years)", type = c("l", "g")) +y_stl <- stlplus( + y, + l.window = 13, + t.window = 19, + s.window = 35, + s.degree = 1, + sub.labels = substr(month.name, 1, 3) +) + +plot( + y_stl, + ylab = "CO2 Concentration (ppm)", + xlab = "Time (years)", + type = c("l", "g") +) plot_seasonal(y_stl) plot_trend(y_stl) plot_cycle(y_stl) plot_rembycycle(y_stl) # if you don't want to use a time series object: -y_stl <- stlplus(y, t = as.vector(time(y)), n.p = 12, - l.window = 13, t.window = 19, s.window = 35, s.degree = 1, - sub.labels = substr(month.name, 1, 3)) +y_stl <- stlplus( + y, + t = as.vector(time(y)), + n.p = 12, + l.window = 13, + t.window = 19, + s.window = 35, + s.degree = 1, + sub.labels = substr(month.name, 1, 3) +) # with an outlier y2 <- co2 y2[200] <- 300 -y2_stl <- stlplus(y2, t = as.vector(time(y2)), n.p = 12, - l.window = 13, t.window = 19, s.window = 35, s.degree = 1, - sub.labels = substr(month.name, 1, 3), outer = 10) +y2_stl <- stlplus( + y2, + t = as.vector(time(y2)), + n.p = 12, + l.window = 13, + t.window = 19, + s.window = 35, + s.degree = 1, + sub.labels = substr(month.name, 1, 3), + outer = 10 +) plot(y2_stl, ylab = "CO2 Concentration (ppm)", xlab = "Time (years)") plot_seasonal(y2_stl) @@ -151,9 +240,16 @@ plot_rembycycle(y2_stl) # compare to R's stl -x1 <- stlplus(co2, t = as.vector(time(co2)), n.p = 12, - l.window = 13, t.window = 19, s.window = 11, s.degree = 1, - sub.labels = substr(month.name, 1, 3)) +x1 <- stlplus( + co2, + t = as.vector(time(co2)), + n.p = 12, + l.window = 13, + t.window = 19, + s.window = 11, + s.degree = 1, + sub.labels = substr(month.name, 1, 3) +) x2 <- stl(co2, l.window = 13, t.window = 19, s.window = 11, s.degree = 1) @@ -161,23 +257,39 @@ x2 <- stl(co2, l.window = 13, t.window = 19, s.window = 11, s.degree = 1) plot(seasonal(x1) - seasonal(x2)) # but not if all jump parameters are 1 -x1 <- stlplus(co2, t = as.vector(time(co2)), n.p = 12, - l.window = 13, t.window = 19, s.window = 11, s.degree = 1, +x1 <- stlplus( + co2, + t = as.vector(time(co2)), + n.p = 12, + l.window = 13, + t.window = 19, + s.window = 11, + s.degree = 1, sub.labels = substr(month.name, 1, 3), - s.jump = 1, t.jump = 1, l.jump = 1) - -x2 <- stl(co2, l.window = 13, t.window = 19, s.window = 11, s.degree = 1, - s.jump = 1, t.jump = 1, l.jump = 1) + s.jump = 1, + t.jump = 1, + l.jump = 1 +) + +x2 <- stl( + co2, + l.window = 13, + t.window = 19, + s.window = 11, + s.degree = 1, + s.jump = 1, + t.jump = 1, + l.jump = 1 +) plot(seasonal(x1) - seasonal(x2)) } -\author{ -Ryan Hafen -} \references{ R. B. Cleveland, W. S. Cleveland, J. E. McRae, and I. Terpenning (1990) STL: A Seasonal-Trend Decomposition Procedure Based on Loess. \emph{Journal of Official Statistics}, \bold{6}, 3--73. } \seealso{ \code{\link{plot.stlplus}} for plotting the components, \code{\link{getraw}}, \code{\link{seasonal}}, \code{\link{trend}}, \code{\link{remainder}} for accessing the components. } - +\author{ +Ryan Hafen +} From b4df78d3a1c9c4e4b0ae85ef5048d848fa903b9e Mon Sep 17 00:00:00 2001 From: Ryan Hafen Date: Mon, 13 Apr 2026 16:00:33 -0700 Subject: [PATCH 4/5] Update NEWS --- .Rbuildignore | 2 +- NEWS.md | 10 ++++++++++ 2 files changed, 11 insertions(+), 1 deletion(-) diff --git a/.Rbuildignore b/.Rbuildignore index e0a855b..6a02063 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -6,7 +6,6 @@ man-roxygen ..Rcheck ^src/.*\.(o|so)$ LICENSE\.md -NEWS\.md ^_ignore$ \.sublime\- sftp-config.json @@ -15,3 +14,4 @@ sftp-config.json .tar.gz \.Rapp\.history cran-comments.md +^CRAN-SUBMISSION$ diff --git a/NEWS.md b/NEWS.md index a1509e6..f07f6b9 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,13 @@ +Version 0.5.2 +---------------------------------------------------------------------- + +- Fixed Rcpp namespace import so Rcpp symbols are properly loaded at runtime +- Updated testthat comparisons against `stats::stl()` to use explicit matched settings and platform-robust tolerances +- Added periodic-mode regression test coverage +- Corrected centering in `plot_seasonal()` so cycle-subseries are centered by the seasonal component mean +- Enabled explicit native routine registration (`useDynLib` with `.registration = TRUE`) +- Removed stale `LazyData` field from DESCRIPTION + Version 0.5 ---------------------------------------------------------------------- From 005174080d4be5d1dbd549defb50b45d37c0e638 Mon Sep 17 00:00:00 2001 From: Ryan Hafen Date: Mon, 13 Apr 2026 20:44:04 -0700 Subject: [PATCH 5/5] Address cran submission issues --- DESCRIPTION | 1 - README.md | 8 ++++---- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index edecfc4..5fa2aa7 100755 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -3,7 +3,6 @@ Type: Package Title: Enhanced Seasonal Decomposition of Time Series by Loess Version: 0.5.2 Date: 2026-04-06 -Author: Ryan Hafen [aut, cre] Authors@R: person("Ryan", "Hafen", email = "rhafen@gmail.com", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-5516-8367")) Maintainer: Ryan Hafen diff --git a/README.md b/README.md index 96f281a..855ba74 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,7 @@ ## stlplus [![R-CMD-check](https://github.com/hafen/stlplus/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/hafen/stlplus/actions/workflows/R-CMD-check.yaml) -[![CRAN](http://www.r-pkg.org/badges/version/stlplus)](https://cran.r-project.org/web/packages/stlplus/index.html) +[![CRAN](http://www.r-pkg.org/badges/version/stlplus)](https://CRAN.R-project.org/package=stlplus) ![](http://cranlogs.r-pkg.org/badges/stlplus) ![png](https://cloud.githubusercontent.com/assets/1275592/11673681/b7fce548-9dcf-11e5-8cd2-f3311b501ab9.png) @@ -16,12 +16,12 @@ Here are some of the added features over `stl()`: - Frequency component smoothing beyond seasonal and trend - Plot methods for diagnostics -For (very) experimental inference, prediction, and variance reduction at endpoints, see the [operator](http://github.com/hafen/operator) package. +For (very) experimental inference, prediction, and variance reduction at endpoints, see the [operator](https://github.com/hafen/operator) package. ## References -- [Cleveland, R. B., Cleveland, W. S., McRae, J. E., & Terpenning, I. (1990). STL: A seasonal-trend decomposition procedure based on loess. *Journal of Official Statistics*, 6(1), 3-73.](http://cs.wellesley.edu/~cs315/Papers/stl%20statistical%20model.pdf) -- [Hafen, R. P. "Local regression models: Advancements, applications, and new methods." (2010).](http://ml.stat.purdue.edu/hafen/preprints/Hafen_thesis.pdf) +- Cleveland, R. B., Cleveland, W. S., McRae, J. E., & Terpenning, I. (1990). STL: A seasonal-trend decomposition procedure based on loess. *Journal of Official Statistics*, 6(1), 3-73. +- Hafen, R. P. "Local regression models: Advancements, applications, and new methods." (2010). ## Installation