无需使用stl或decompose提取时间序列季节性效应的方法

你收到报错的核心原因是当前构造的time_series对象频率设置为1，代表年度数据，R的decompose()和stl()默认会认为序列不存在小于1年的传统季节性周期，且没有至少2个完整周期可用于分解，因此触发报错。
注：年度数据本身不存在常规的“年内季节性效应”，你要提取的实际是多年重复的周期性波动规律，以下是可落地的实现方法：

• 方法1：手动指定周期后调用原生分解函数
  如果通过行业先验知识（比如煤业生产周期、经济周期规律）已经确定序列的波动周期长度，直接重新构造时间序列对象，将frequency参数设置为对应周期值即可。比如假设周期为10年，构造代码如下：

time_series <- ts(bicoal, start = 1920, frequency = 10)

构造完成后就可以正常调用decompose(time_series)或者stl(time_series, s.window = "periodic")完成分解，提取周期性（即你所说的季节性）效应。注意frequency的取值需要满足序列总长度 >= 2*frequency，你的数据总长度为49，只要周期设置为24及以下都满足要求。

• 方法2：先通过频谱分析确定周期再分解
  如果不知道合理的周期长度，可以先通过功率谱分析识别序列的潜在周期：

# 绘制功率谱图，峰值对应频率的倒数就是潜在周期
spec_pgram(time_series, spans = c(3,3), main = "序列功率谱图")

找到功率谱最高的几个峰值，计算对应周期后再按照方法1的步骤设置frequency做分解即可。

• 方法3：手动拟合分解序列
  也可以通过回归的方式手动拆分趋势、周期和残差项：

# 第一步：用局部加权回归拟合趋势项
trend_fit <- loess(bicoal ~ seq_along(bicoal), span = 0.7)
trend <- predict(trend_fit)
# 第二步：得到去趋势后的序列
detrend_series <- bicoal - trend
# 第三步：按识别到的周期长度对去趋势序列分组，组内均值就是对应周期位置的季节性/周期性效应
period <- 10 # 替换成你确定的周期长度
cycle_effect <- tapply(detrend_series, rep(1:period, length.out = length(bicoal)), mean, na.rm = T)

• 方法4：调用更灵活的多周期分解函数
  可以使用forecast包中的mstl()函数，该函数支持自定义周期，对低频序列的兼容性比原生stl()更好：

library(forecast)
# 直接指定周期参数即可完成分解
decomp_result <- mstl(time_series, period = 10)
# 提取季节性/周期性效应
cycle_effect <- decomp_result[,"Seasonal1"]

年度数据的“季节性效应”本质是多年商业/生产周期，你需要确认你设定的周期长度有对应的业务逻辑支撑，避免无意义的数值拟合。

内容的提问来源于stack exchange，提问作者rrr