时间序列的数据分析:STL分解

Posted 2023-04-02 -派神-

之前已经完成了三篇关于时间序列的博客，还没有阅读过的读者请先阅读：

1. 时间序列的数据分析(一):主要成分

2. 时间序列的数据分析(二):数据趋势的计算

3. 时间序列的数据分析(三):经典时间序列分解 

六. STL分解

6.1 主要参数

STL（Seasonal and Trend decomposition using Loess）是一个非常通用和稳健强硬的分解时间序列的方法，其中Loess是一种估算非线性关系的方法。STL分解法由 R. B. Cleveland, Cleveland, McRae, & Terpenning (1990)提出。STL也是将时间序列分解成三个主要分量: 趋势、季节项和残差 。STL使用LOESS(locally estimated scatterplot smoothing) 来提取三个分量的平滑估计，在python中实现时间序列的STL分解主要是通过调用statsmodels类库的STL方法来实现的，该STL方法有四个主要的输入参数：

• endog：表示需要分解的数据集,它是STL方法的第一个参数，该数据集的类型可以是numpy的array,也可以是pandas的series 或者dataframe.
• period：表示季节性周期，如果endog的类型是numpy的array则需要指定period，如果是pandas的series 或dataframe则stl方法可以根据索引推断出period,因此无需指定peroid
• season： 表示季节性平滑器的长度，它必须是一个奇数，通常要>=7(默认)。
• trend：表示趋势平滑器的长度,通常要>period(或season)的1-1.5倍，并且它必须是一个奇数。默认值是最小的1-1.5倍的period,比如period=7则trend默认值是9，如果period=12则trend默认值是13

6.2 分解过程

下面我们使用statsmodels的STL方法对航空公司乘客数据进行分解并获取各个分量的结果：

from statsmodels.tsa.seasonal import STL
plt.rc("figure", figsize=(10, 6))

df=pd.read_csv("airline_Passengers.csv")
df['Period']=pd.to_datetime(df['Period'])
df.set_index('Period',inplace=True)

res = STL(df).fit()
res.plot()

df['trend']=res.trend
df['seasonal']=res.seasonal
df['resid']=res.resid

 

 这里的STL方法中我们只使用了第一个参数，其它均为默认参数,因为我们的数据集是dataframe，因此STL方法可以根据datetime的索引列推断出peroid，如果数据类型是numpy的array那就必须指定peroid。下面我们可以观察一下残差的分布以及它的均值，一般情况下如果残差呈现出以0为均值的近似正太分布(这不是必须的)那么说明我们使用了正确的分解方法。

print('residual mean:',df.resid.mean())
df.resid.hist();

 从上面的结果可知我们的残差近似正太分布并且均值在0的附近,这说明SLT分解是正确的。

 6.3 趋势性、季节性程度及季节项波峰的计算

时间序列数据可以被分解为:趋势(Trend)、季节性(seasonal)、残差(residual),其分解式一般可以表示为:

 其中T(t)表示t时刻的趋势值,S(t)表示t时刻的季节项值,R(t)表示t时刻的残差值。对于趋势性很强的数据，经季节调整后(删除季节项)的数据应比残差项的变动幅度更大。因此，会相对较小。但是，对于没有趋势或是趋势很弱的时间序列，两个方差应大致相同。因此，我们将趋势强度定义为：

 这可以给趋势强度的衡量标准，其值在 0-1 之间。因为有些情况下残差项的方差甚至比季节变换后的序列还大，我们令可取的最小值为0。

相似地，季节性的强度定义如下，其所用的数据为去除趋势后的数据而不是去除季节后的数据。

当季节强度接近 0 时表示该序列几乎没有季节性，当季节强度接近 1 时表示该序列的 远小于 。

在时间序列中季节性一般呈现周期性变化的规律,因此季节性周期中的波峰大体上也是固定的，因此我们只需要找到季节性周期中的最大值就可以确定波峰期。

下面我们来计算一下趋势程度、季节性程度以及季节性波峰期，首先我们需要在数据中删除趋势项和季节项并得到两个新列:detrend和deseasonal，其中detrend列表示, 而deseasonal表示:

#从数据中删除趋势项
df['detrend']=df['#Passengers']-df.trend
#从数据中删除季节项
df['deseasonal']=df['#Passengers']-df.seasonal

接下来我们套用公式来计算趋势和季节性程度：

trend_strength=max(0,1-df.resid.var()/df.deseasonal.var())
seasonal_strength=max(0,1-df.resid.var()/df.detrend.var())
print('trend_strength:',trend_strength)
print('seasonal_strength:',seasonal_strength)

 

 从结果中我们看到数据中的趋势和季节性程度都非常高(接近1)，趋势和季节性程度越高,那说明数据的可预测性越好。接下来我们来计算季节性波峰：

period=12
peak = (np.argmax(df.seasonal) + 1) % period
peak = period if peak == 0 else peak

print("peak:",peak)

 

波峰值为7，说明改每年的7月为波峰期,这个从数据趋势图中也能得到确认。

总结

今天我们主要介绍了STL的分解的主要参数,和分解的过程，并观察了分解以后残差的分布和均值并确认了残差服从以0为均值的近似正太分布,这说明STL分解是正确的。其次我们还介绍了趋势程度、季节性程度以及季节性波峰的计算方法，这有助于确定数据是否具有良好的可预测性。

参考资料

statsmodels.tsa.seasonal.STL — statsmodels 

Seasonal-Trend decomposition using LOESS (STL) — statsmodels

https://www.scb.se/contentassets/ca21efb41fee47d293bbee5bf7be7fb3/stl-a-seasonal-trend-decomposition-procedure-based-on-loess.pdf

用于异常检测的具有缺失值的时间序列的 STL 分解

【中文标题】用于异常检测的具有缺失值的时间序列的 STL 分解

【英文标题】：STL decomposition of time series with missing values for anomaly detection

【发布时间】：2012-08-17 00:09:57

【问题描述】：

我正在尝试检测时间序列气候数据中的异常值，其中一些观测值缺失。在网上搜索我发现了许多可用的方法。其中，stl 分解似乎很有吸引力，因为它去除了趋势和季节性成分并研究了其余部分。阅读STL: A Seasonal-Trend Decomposition Procedure Based on Loess，stl 在确定分配可变性的设置方面似乎很灵活，不受异常值的影响，并且尽管有缺失值也可以应用。但是，尝试在 R 中应用它，经过四年的观察并根据 http://stat.ethz.ch/R-manual/R-patched/library/stats/html/stl.html 定义所有参数，我遇到了错误：

时间序列包含内部 NAs

当na.action = na.omit，并且

序列不是周期性的或少于两个周期

当na.action = na.exclude。

我已经仔细检查了频率是否正确定义。我在博客中看到了相关问题，但没有找到任何可以解决此问题的建议。是否可以在缺少值的系列中应用 stl？我非常不愿意插入它们，因为我不想引入（并因此检测......）伪影。出于同样的原因，我不知道改用 ARIMA 方法有多么可取（如果缺失值仍然是一个问题）。

如果您知道在缺失值的系列中应用 stl 的方法，或者您认为我的选择在方法上不合理，或者您有更好的建议，请分享。我是该领域的新手，被大量（看似......）相关信息所淹没。

【问题讨论】：

你有多少缺失值？您可以尝试 zoo 包中的 na.approx 函数来插入丢失的数据。

一些额外的想法：显然，现在包 stlplus 可以处理 ts 分解中的缺失值。 Rbeast 可能是另一种选择；它是贝叶斯算法，不像stl只分解时间序列，Rbeast同时进行时间序列分解和变化点检测，允许缺失值。这是一个例子：library(Rbeast); co2[ sample(1:length(co2), 200) ]=NA; plot(beast(co2)).

【参考方案1】：

stl开头我们find

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 在 na.action(as.ts(x)) 之后成为 ts 对象（否则 period == 1）。让我们先检查na.omit 和na.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"

这是一种不同的情况。如上所述，stl 期望 na.action(as.ts(x)) 为 ts，但 na.exclude(as.ts(x)) 属于 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 值（即使用 na.action = na.pass），它在 Fortran 中崩溃更深（参见完整源代码 here）：

z <- .Fortran(C_stl, ...

na.new 的替代品并不令人愉快：

na.contaguous - 查找时间序列对象中最长的连续非缺失值。 na.approx、na.locf 来自 zoo 或其他一些插值函数。 不确定这个，但是可以找到另一个用于 Python here 的 Fortran 实现。可以使用 Python 进行一些修改后从源代码安装 R，以防此模块确实允许缺失值。

正如我们在paper 中看到的那样，在调用stl 之前，没有一些简单的缺失值程序（比如在一开始就近似它们）可以应用于时间序列。因此，考虑到original implementation 相当冗长这一事实，我会考虑一些其他替代方案，而不是全新的实现。

更新：当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 中看到了一些趋势差异。但考虑到样本量 (72)，这些差异在案例 1 中非常小，在案例 3 中也令人满意。

【讨论】：

对，看来stl 在 R 中的实现是基于问题中链接的论文第 4 节（第 20 页）中讨论的 Fortran 实现，它避开了设计目标 4（第 5 页） ) 有利于速度（见第 21 页底部）。我希望有一个允许缺失值的解决方案（S 显然存在一个实现）。尽管我怀疑在大多数情况下，在 stl 之前通过 na.approx 传递输入之间的差异与处理缺失值的 stl 实现并没有太大差异，这确实是我的后备选项。

@user1935457，我更新了我的答案。显然，R 中 stl 的实现缺少一些部分，但正如我在回答中提到的，可能不值得重写它，特别是如果它真的不会像你说的那样有太大差异。只是考虑了一些具有完整数据和na.approx 的测试，我将用结果更新我的答案..

【参考方案2】：

意识到这是一个老问题，但我想我会更新，因为 R 中有一个更新的 stl 包，称为 stlplus。 Here is its homepage on github。您可以使用 install.packages("stlplus") 从 CRAN 安装它，或者使用 devtools::install_github("hafen/stlplus") 直接从 github 安装它。