一开始stl我们发现
x <- na.action(as.ts(x))
不久之后
period <- frequency(x)
if (period < 2 || n <= 2 * period)
stop("series is not periodic or has less than two periods")
也就是说，stl期望在（否则）之后x成为ts对象。让我们先检查一下。na.action(as.ts(x))period == 1na.omitna.exclude
显然，最后getAnywhere("na.omit.ts")我们发现
if (any(is.na(object)))
stop("time series contains internal NAs")
这很简单，什么都做不了（na.omit不排除NAs对象ts）。现在getAnywhere("na.exclude.default")排除NA观察，但返回一个类对象exclude：
attr(omit, "class") <- "exclude"
这是另一种情况。如上所述，stlexpects na.action(as.ts(x))to be ts，但na.exclude(as.ts(x))is of class exclude。
NAs因此，如果一个人对排除感到满意，那么例如
nottem[3] <- NA
frequency(nottem)
# [1] 12
na.new <- function(x) ts(na.exclude(x), frequency = 12)
stl(nottem, na.action = na.new, s.window = "per")
作品。通常，stl不适用于NA值（即 with na.action = na.pass），它在 Fortran 中崩溃得更深（请参阅此处的完整源代码）：
z <- .Fortran(C_stl, ...
替代方案na.new并不令人愉快：
na.contaguous- 查找时间序列对象中最长的连续非缺失值。na.approx, na.locffromzoo或其他一些插值函数。- 不确定这个，但可以在此处找到 Python 的另一个 Fortran 实现。经过一些修改后，可以使用 Python 可能从源代码安装 R，以防此模块确实允许缺失值。
正如我们在论文中看到的那样，没有一些简单的缺失值过程（比如在一开始就近似它们）可以在调用之前应用于时间序列stl。因此，考虑到原始实现相当冗长这一事实，我会考虑一些其他替代方案，而不是全新的实现。
更新：NAs当有可能na.approx来自时，在许多方面都是一个非常理想的选择zoo，所以让我们检查它的性能，即比较 的结果stl与完整数据集的结果，以及当有一定数量的 时的结果NAs，使用na.approx。我使用MAPE作为准确度的度量，但仅用于趋势，因为季节性分量和余数过零并且会扭曲结果。的位置NAs是随机选择的。
library(zoo)
library(plyr)
library(reshape)
library(ggplot2)
mape <- function(f, x) colMeans(abs(1 - f / x) * 100)
stlCheck <- function(data, p = 3, ...){
set.seed(20130201)
pos <- lapply(3^(0:p), function(x) sample(1:length(data), x))
datasetsNA <- lapply(pos, function(x) {data[x] <- NA; data})
original <- data.frame(stl(data, ...)$time.series, stringsAsFactors = FALSE)
original$id <- "Original"
datasetsNA <- lapply(datasetsNA, function(x)
data.frame(stl(x, na.action = na.approx, ...)$time.series,
id = paste(sum(is.na(x)), "NAs"),
stringsAsFactors = FALSE))
stlAll <- rbind.fill(c(list(original), datasetsNA))
stlAll$Date <- time(data)
stlAll <- melt(stlAll, id.var = c("id", "Date"))
results <- data.frame(trend = sapply(lapply(datasetsNA, '[', i = "trend"), mape, original[, "trend"]))
results$id <- paste(3^(0:p), "NAs")
results <- melt(results, id.var = "id")
results$x <- min(stlAll$Date) + diff(range(stlAll$Date)) / 4
results$y <- min(original[, "trend"]) + diff(range(original[, "trend"])) / (4 * p) * (0:p)
results$value <- round(results$value, 2)
ggplot(stlAll, aes(x = Date, y = value, colour = id, group = id)) + geom_line() +
facet_wrap(~ variable, scales = "free_y") + theme_bw() +
theme(legend.title = element_blank(), strip.background = element_rect(fill = "white")) +
labs(x = NULL, y = NULL) + scale_colour_brewer(palette = "Set1") +
lapply(unique(results$id), function(z)
geom_text(data = results, colour = "black", size = 3,
aes(x = x, y = y, label = paste0("MAPE (", id, "): ", value, "%"))))
}
nottem, 240 次观察
stlCheck(nottem, s.window = 4, t.window = 50, t.jump = 1)
co2, 468 次观察
stlCheck(log(co2), s.window = 21)
mdeaths, 72 个观测值
stlCheck(mdeaths, s.window = "per")
在视觉上，我们确实看到案例 1 和案例 3 的趋势存在一些差异。但考虑到样本量，这些差异在案例 1 中非常小，在案例 3 中也令人满意 (72)。