R中tsibble做STL时序分解补零致趋势线失真的解决方法

我正在开展时间序列预测相关实践，使用R语言tsibble包，参考Hyndman和Athanasopoulos编写的经典时序教材进行练习，目前使用的是小时级订单时间序列数据，仅包含每日9:00-20:00营业时段的观测值，初始数据结构为符合规范的tsibble对象，前100行数据结构如下：

structure(list(X = 1:100, Time = structure(c(1546351200, 1546354800, 
1546358400, 1546362000, 1546365600, 1546369200, 1546372800, 1546376400, 
1546380000, 1546383600, 1546387200, 1546390800, 1546437600, 1546441200, 
1546444800, 1546448400, 1546452000, 1546455600, 1546459200, 1546462800, 
1546466400, 1546470000, 1546473600, 1546477200, 1546524000, 1546527600, 
1546531200, 1546534800, 1546538400, 1546542000, 1546545600, 1546549200, 
1546552800, 1546556400, 1546560000, 1546563600, 1546610400, 1546614000, 
1546617600, 1546621200, 1546624800, 1546628400, 1546632000, 1546635600, 
1546639200, 1546642800, 1546646400, 1546650000, 1546696800, 1546700400, 
1546704000, 1546707600, 1546711200, 1546714800, 1546718400, 1546722000, 
1546725600, 1546729200, 1546732800, 1546736400, 1546783200, 1546786800, 
1546790400, 1546794000, 1546797600, 1546801200, 1546804800, 1546808400, 
1546812000, 1546815600, 1546819200, 1546822800, 1546869600, 1546873200, 
1546876800, 1546880400, 1546884000, 1546887600, 1546891200, 1546894800, 
1546898400, 1546902000, 1546905600, 1546909200, 1546956000, 1546959600, 
1546963200, 1546966800, 1546970400, 1546974000, 1546977600, 1546981200, 
1546984800, 1546988400, 1546992000, 1546995600, 1547042400, 1547046000, 
1547049600, 1547053200), tzone = "", class = c("POSIXct", "POSIXt"
)), Orders = c(390.9300738, 424.5024938, 459.9507418, 493.879574, 
521.1915476, 543.6420076, 564.1188556, 583.1138192, 599.9870792, 
608.2502946, 589.0774506, 552.9864864, 460.0146478, 513.6096, 
565.54751, 614.3622836, 649.610842, 673.080916, 694.1457822, 
714.687121, 730.0065136, 727.6420116, 704.9715348, 669.8276592, 
596.9598262, 627.6943506, 663.224885, 689.6623702, 705.7821348, 
705.2804398, 702.4425002, 686.257045, 673.0440842, 631.4105592, 
590.2226836, 557.170647, 505.5489378, 514.1362306, 518.0591858, 
519.7312244, 515.538957, 517.0255898, 516.4563428, 519.1586616, 
518.8452174, 494.1823666, 468.0562396, 444.6603772, 465.6368096, 
475.6484144, 481.8889642, 489.3110196, 492.9861102, 495.1878822, 
496.2013992, 500.4736856, 502.9525222, 490.2884824, 465.9459928, 
446.4332428, 468.827488, 475.2297188, 480.3550016, 486.8966308, 
488.641556, 492.2285006, 493.485411, 501.271116, 501.7387056, 
485.8849556, 462.3912654, 444.3381798, 423.7514376, 442.9296904, 
456.1299334, 459.6065968, 466.132121, 468.2358706, 476.4634124, 
481.580409, 484.5936224, 477.3972256, 458.546062, 439.0321916, 
397.7730592, 418.1761574, 426.6949568, 435.0296632, 438.8624322, 
437.5872586, 441.9915442, 445.0284556, 443.4291354, 445.9624284, 
430.1143198, 420.7732792, 485.8293664, 494.2056144, 502.3287016, 
509.0143842)), class = c("tbl_ts", "tbl_df", "tbl", "data.frame"
), row.names = c(NA, -100L), key = structure(list(.rows = structure(list(
    1:100), ptype = integer(0), class = c("vctrs_list_of", "vctrs_vctr", 
"list"))), class = c("tbl_df", "tbl", "data.frame"), row.names = c(NA, 
-1L)), index = structure("Time", ordered = TRUE), index2 = "Time", interval = structure(list(
    year = 0, quarter = 0, month = 0, week = 0, day = 0, hour = 1, 
    minute = 0, second = 0, millisecond = 0, microsecond = 0, 
    nanosecond = 0, unit = 0), .regular = TRUE, class = c("interval", 
"vctrs_rcrd", "vctrs_vctr")))

问题复现 #

首先绘制原始时序图后，尝试通过STL分解提取趋势、季节性、周期成分，初始建模代码如下：

dcmp <- sales %>%
  model(stl = STL(Orders))

上述代码可正常运行，但后续调用components()提取分解结果时报错：

Error in `transmute()`:
! Problem while computing `cmp = map(.fit, components)`.
Caused by error in `UseMethod()`:
! no applicable method for 'components' applied to an object of class "null_mdl"
Backtrace:
  1. generics::components(dcmp)
  8. fabletools:::map(.fit, components)
  9. base::lapply(.x, .f, ...)
 11. fabletools:::components.mdl_ts(X[[i]], ...)
 12. generics::components(object$fit, ...)

检索报错信息后得知问题源于时间序列存在时间缺口，因此使用fill_gaps()补全所有小时级时间点，将补全后生成的NA值统一替换为0，代码如下：

sales <- sales %>%
  fill_gaps()

sales$Orders[is.na(sales$Orders)] <- 0

处理后可正常完成STL分解并绘图，绘图代码如下：

dcmp <- sales %>%
  model(stl = STL(Orders))

dcmp <- components(dcmp)

dcmp %>%
  as_tsibble() %>%
  autoplot(Orders, colour="gray") +
  geom_line(aes(y=trend), colour = "#D55E00") +
  labs(
    y = "Orders",
    title = "Orders") + 
  labs(caption = "fake data")

但此时得到的趋势线存在严重偏差：橙色趋势线因纳入了非营业时段的0值被大幅拉低，远低于实际订单的真实趋势水平，不符合预期。教材示例中的趋势线贴近真实观测值，不会被无效零值拉低。

核心诉求 #

需要基于tsibble、fable等教材所用tidyverts时序生态包给出解决方案：如何处理仅覆盖每日营业时段、非营业时段无观测的小时级时间序列，才能正确完成STL分解，得到准确的趋势、季节性成分，且不使用手动移动平均替代原生STL算法的结果。

问题根源是将非营业时段的结构性缺失误判为真实0值：STL是基于等间隔连续时间点的局部加权回归算法，强行插入无业务意义的0值必然会拉低趋势估计。正确的处理逻辑是不要强行补全为24小时连续小时序列，而是将序列的时间索引重构为营业时序索引，让观测点在业务逻辑上保持等间隔，同时保留小时、日期等季节性维度供STL识别周期。
具体实现步骤：

1. 加载依赖包，提取原始数据的时间维度特征

   library(tsibble)
   library(fable)
   library(dplyr)
   library(lubridate)
   library(ggplot2)
   
   sales <- sales %>%
     mutate(
       date = as_date(Time),
       hour = hour(Time),
       # 标记日内观测序号：9点为第1个，10点为第2个，以此类推到20点为第12个
       obs_in_day = match(hour, sort(unique(hour))),
       # 标记营业日序号
       day_id = as.integer(date - min(date) + 1)
     )
   

2. 重构tsibble索引，使用业务层面的连续等间隔索引替代自然时间索引，消除无意义的时间缺口

   sales <- sales %>%
     # 按日期、日内观测点排序，保证时序顺序正确
     arrange(date, obs_in_day) %>%
     mutate(
       # 营业时序索引：每个营业小时递增1，全局连续无缺口
       business_idx = row_number()
     ) %>%
     # 重新构建规范tsibble，索引设为营业时序索引
     as_tsibble(index = business_idx, regular = TRUE)
   

3. 调用原生STL算法，手动指定季节性周期长度，无需依赖自然时间间隔自动识别
   • 日内小时周期：每日共12个营业小时，周期长度设为12
   • 周度日周期：每周7天*每日12个营业小时，周期长度设为84
     可根据实际数据的周期规律调整窗口参数：

   dcmp <- sales %>%
     model(
       stl = STL(
         Orders ~ trend(window = 21) +
                 season(period = 12, window = 13) + # 日内小时级季节项
                 season(period = 84, window = 25),  # 周内日度季节项
         robust = TRUE
       )
     ) %>%
     components()
   

4. 可视化分解结果时，使用自然时间Time作为x轴即可，此时得到的趋势、季节项完全基于真实营业观测计算，不会被无效0值干扰：

   dcmp %>%
     autoplot(Orders, colour="gray") +
     geom_line(aes(y=trend, x=Time), colour = "#D55E00", linewidth=1) +
     labs(
       y = "订单量",
       title = "营业