模型介绍

TBATS 模型

TBATS(Trend and Breaks, Adaptive, Truncated, Seasonal) 模型是一种用于时间序列预测的方法，它的名称来源于模型的几个关键组成部分：

Trigonometric seasonality（三角季节性）、Box-Cox transformation（Box-Cox转换）、ARMA errors（ARMA误差）、Trend（趋势）和Seasonal components（季节性成分）。

该模型特别适用于处理具有复杂季节性模式的时间序列数据，例如那些周期不固定或者包含多个季节性周期的数据。

主要特点

TBATS模型的主要特点是它能够通过指数平滑方法来处理复杂的季节性问题。

它可以对时间序列的周期进行多种建模选择，包括是否进行Box-Cox转换、是否考虑趋势、是否使用趋势衰减以及是否对模型残差使用ARIMA(p,q)。

此外，TBATS模型还可以处理非整数和长周期的季节性，这在其他模型中往往难以实现。

ETS 模型

ETS(Exponential Smoothing State Space Model) 模型是一组具有基础状态空间模型的时间序列模型，它是基于指数平滑法，它的名称来源于三个核心成分：

• E (Error): 误差项（加性或乘性）。

• T (Trend): 趋势（无趋势、加性趋势、乘性趋势）。

• S (Seasonal): 季节性（无季节性、加性季节性、乘性季节性）

这三个部分的任意组合构成了ETS模型，ETS可以理解为Error,Trend,Seasonality,也可以理解为Exponen Tial Smoothing模型。ETS模型采用最大似然法进行参数估计。

误差项为相加模型时，最大似然法等价于使SSE最小来进行参数估计。该模型会自动尝试不同的ETS各项组合并通过使AICc最小来选择最佳模型。

STL 模型

STL(Seasonal and Trend decomposition using Loess) 模型是一种通过局部加权回归（LOESS）方法对时间序列进行分解的技术，将数据分为季节性、趋势和残差三部分。它以灵活性和鲁棒性著称，能够处理非线性趋势、非固定季节性和带噪声的数据，适合分析复杂的时间序列。STL允许用户调整窗口宽度以平衡分解的细节程度，适用于各种周期的时间序列（如月度、季度）。模型广泛用于数据分析和预测，特别适合带有长期趋势和周期性波动的场景，但对快速变化或非周期性数据的效果可能有限。

NNETAR模型

NNETAR(Neural Network Autoregressive model) 模型是一种基于神经网络的时间序列预测方法，通过将时间序列建模为自回归问题，利用神经网络捕捉数据中的非线性关系，模型会自动选择滞后阶数和隐藏层神经元数量，适合具有非平稳性或复杂模式的数据。NNETAR能灵活处理季节性和噪声，适用于非线性动态系统的短期预测。

XGBoost模型

XGBoost(Extreme Gradient Boosting) 是一种基于梯度提升框架的高效机器学习模型，以构建和优化决策树为核心，它在大规模、高效、准确性上具有优势。它的原理如下：

梯度提升树（Gradient Boosting Tree）：梯度提升树是一种集成学习算法，通过迭代地训练一系列弱分类器（决策树），每一次迭代都试图拟合前一次迭代的残差。最终将这些弱分类器进行组合，得到一个更强大的模型。
正则化：XGBoost引入了正则化项来控制模型的复杂度，包括L1和L2正则化。
优化策略：XGBoost在损失函数中引入了泰勒展开近似，使用一阶和二阶导数的信息来加速训练过程。此外，它还使用了近似的贪婪算法来选择最佳切分点。

综合模型

以TBATS 模型、ETS 模型（指数平滑模型）、STL 模型（季节分解）和NNETAR模型（神经网络自回归）作为输入，以真实的数据作为输出，训练XGBoost，并进行预测。

研究思路

收集JS地区1990—2020年的每月的乙肝发病率数据，然后将数据分为训练集和测试集，其中将1990年1月至2015年12月的数据作为训练集的数据，将剩余的数据作为测试集的数据。首先，从Excel文件中读取乙肝发病率数据，并将年份转换为日期格式。接着，将发病率数据构建为时间序列对象，并按照时间划分为训练集（1990年1月到2015年12月）和测试集（2016年1月到2020年12月）。

单一模型构建

使用四种不同的时间序列预测模型进行建模：

TBATS模型：适合处理具有复杂季节性和趋势特性的时间序列。
ETS模型：基于指数平滑方法，适合处理具有趋势和季节性但较简单的时间序列。
STL分解模型：通过对时间序列分解为趋势、季节性和残差部分，分别进行建模。
NNETAR模型：使用神经网络自回归建模，能够捕捉时间序列中的非线性关系。

组合模型

考虑到单一模型的预测能力可能具有局限性，进一步构建了一个基于 XGBoost的组合模型。其核心思路是将单一模型的预测值作为特征，与时间序列的时间特征结合，生成新的训练数据集。使用XGBoost进行建模时，通过超参数调优和交叉验证技术，以确定最佳的迭代次数，从而优化模型性能。最终通过XGBoost模型对测试集进行预测，实现对单一模型预测结果的集成。

结果可视化与模型比较

将测试集的实际值与各单一模型的预测值及组合模型的预测值进行对比，并通过折线图显示不同模型的预测表现。计算训练集和测试集上的均方根误差（RMSE），量化各模型的预测误差，从而评估不同模型的性能。

建模过程

# 加载必要的库
library(forecast)
library(ggplot2)
library(caret)
library(xgboost) 
library(stargazer)
library(cowplot)

# 数据读取
df <- readxl::read_excel('HBV--JS.xlsx')

# 转换年份为日期类型
df$年份 <- as.Date(df$年份)
ts_data <- ts(df$发病率, start = c(1990, 1), frequency = 12)  

# --- 划分训练集和测试集 ---
train_data <- window(ts_data, end = c(2015, 12))
test_data <- window(ts_data, start = c(2016, 1))

— 构建单一模型 —

# TBATS 模型
tbats_model <- tbats(train_data)
tbats_forecast <- forecast(tbats_model, h = length(test_data))

# ETS 模型（指数平滑模型）
ets_model <- ets(train_data)
ets_forecast <- forecast(ets_model, h = length(test_data))

# STL 模型（季节分解）
stlm_model <- stlm(train_data, s.window = "periodic")
stlm_forecast <- forecast(stlm_model, h = length(test_data))

# NNETAR 模型（神经网络自回归）
set.seed(1234)
nnetar_model <- nnetar(train_data)
nnetar_forecast <- forecast(nnetar_model, h = length(test_data))

— 构建组合模型 —

# 增加时间特征
ensemble_data_train <- data.frame(
  tbats = as.numeric(fitted(tbats_model)),
  ets = as.numeric(fitted(ets_model)),
  stlm = as.numeric(fitted(stlm_model)),
  nnetar = as.numeric(fitted(nnetar_model)),
  month = as.numeric(format(seq(as.Date("1990-01-01"), by = "month", length.out = length(train_data)), "%m")),
  year = as.numeric(format(seq(as.Date("1990-01-01"), by = "month", length.out = length(train_data)), "%Y")),
  actual = as.numeric(train_data)
)

# 准备特征和目标变量
train_matrix <- xgb.DMatrix(data = as.matrix(ensemble_data_train[, -7]), label = ensemble_data_train$actual)

# 设置XGBoost参数
params <- list(
  objective = "reg:squarederror",
  eta = 0.01,  # 降低学习率
  max_depth = 20,  # 增加深度
  subsample = 0.9,  # 增加随机采样
  colsample_bytree = 0.9,  # 增加特征采样
  lambda = 1,  # L2 正则化
  alpha = 0.5  # L1 正则化
)

# 使用交叉验证选择最佳迭代次数
set.seed(123)
cv <- xgb.cv(
  params = params,
  data = train_matrix,
  nrounds = 3000,
  nfold = 5,
  verbose = 0,
  early_stopping_rounds = 10
)

# 使用最佳迭代次数训练模型
best_nrounds <- cv$best_iteration
xgb_model <- xgb.train(
  params = params,
  data = train_matrix,
  nrounds = best_nrounds,
  verbose = 0
)

# 准备测试集数据
ensemble_data_test <- data.frame(
  tbats = as.numeric(tbats_forecast$mean),
  ets = as.numeric(ets_forecast$mean),
  stlm = as.numeric(stlm_forecast$mean),
  nnetar = as.numeric(nnetar_forecast$mean),
  month = as.numeric(format(seq(as.Date("2016-01-01"), by = "month", length.out = length(test_data)), "%m")),
  year = as.numeric(format(seq(as.Date("2016-01-01"), by = "month", length.out = length(test_data)), "%Y"))
)
test_matrix <- xgb.DMatrix(data = as.matrix(ensemble_data_test))

# 组合模型预测
xgb_forecast <- predict(xgb_model, newdata = test_matrix)

# 使用 XGBoost 模型拟合训练集
train_combined <- data.frame(
  tbats = as.numeric(fitted(tbats_model)),
  ets = as.numeric(fitted(ets_model)),
  stlm = as.numeric(fitted(stlm_model)),
  nnetar = as.numeric(fitted(nnetar_model)),
  month = as.numeric(format(seq(as.Date("1990-01-01"), by = "month", length.out = length(train_data)), "%m")),
  year = as.numeric(format(seq(as.Date("1990-01-01"), by = "month", length.out = length(train_data)), "%Y"))
)

train_combined$Combined <- predict(xgb_model, newdata = train_matrix)

# 构建训练集数据框
train_forecast_df <- data.frame(
  时间 = seq(as.Date("1990-01-01"), by = "month", length.out = length(train_data)),
  TBATS = as.numeric(fitted(tbats_model)),
  ETS = as.numeric(fitted(ets_model)),
  STL = as.numeric(fitted(stlm_model)),
  NNETAR = as.numeric(fitted(nnetar_model)),
  Combined = train_combined$Combined,
  Actual = as.numeric(train_data) 
)
# 绘制训练集图
ggplot(train_forecast_df, aes(x = 时间)) +
  geom_line(aes(y = TBATS, color = "TBATS"), linewidth = 1) +
  geom_line(aes(y = ETS, color = "ETS"), linewidth = 1, linetype = "dotted") +
  geom_line(aes(y = STL, color = "STL"), linewidth = 1, linetype = "dashed") +
  geom_line(aes(y = NNETAR, color = "NNETAR"), linewidth = 1, linetype = "dotdash") +
  geom_line(aes(y = Combined, color = "Combined"), linewidth = 1.1, linetype = "solid") +
  geom_line(aes(y = Actual, color = "Actual"), linewidth = 1.1, linetype = "solid") +
  labs(
    title = "训练集预测结果对比",
    subtitle = "TBATS、ETS、STL、NNETAR、组合模型与实际值对比",
    y = "发病率",
    x = "时间",
    color = "模型类型"
  ) +
  theme_bw(base_size = 14) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5),
    legend.position = "bottom",
    legend.title = element_text(face = "bold"),
    legend.text = element_text(size = 12),
    legend.spacing.x = unit(0.5, 'cm'),
    legend.background = element_rect(fill = "white", color = "black"),
    panel.grid.major = element_line(color = "grey", linetype = "dashed"),
    panel.grid.minor = element_blank()
  ) +
  scale_color_manual(
    values = c(
      "TBATS" = "blue",
      "ETS" = "red",
      "STL" = "green",
      "NNETAR" = "purple",
      "Combined" = "black",
      "Actual" = "orange"
    )
  ) +
  guides(
    color = guide_legend(
      override.aes = list(linewidth = 2)
    )
  )

# 准备测试集预测结果数据框
forecast_df <- data.frame(
  时间 = seq(as.Date("2016-01-01"), by = "month", length.out = length(test_data)),
  TBATS = as.numeric(tbats_forecast$mean),
  ETS = as.numeric(ets_forecast$mean),
  STL = as.numeric(stlm_forecast$mean),
  NNETAR = as.numeric(nnetar_forecast$mean),
  Combined = xgb_forecast, 
  Actual = as.numeric(test_data)  
)

# 绘制测试集图
ggplot(forecast_df, aes(x = 时间)) +
  geom_line(aes(y = TBATS, color = "TBATS"), linewidth = 1) +
  geom_line(aes(y = ETS, color = "ETS"), linewidth = 1, linetype = "dotted") +
  geom_line(aes(y = STL, color = "STL"), linewidth = 1, linetype = "dashed") +
  geom_line(aes(y = NNETAR, color = "NNETAR"), linewidth = 1, linetype = "dotdash") +
  geom_line(aes(y = Combined, color = "Combined"), linewidth = 1.5, linetype = "solid") +
  geom_line(aes(y = Actual, color = "Actual"), linewidth = 1.2, linetype = "solid") +
  labs(
    title = "测试集预测结果对比",
    subtitle = "TBATS、ETS、STL、NNETAR、组合模型与实际值对比",
    y = "发病率",
    x = "时间",
    color = "模型类型"
  ) +
  theme_bw(base_size = 14) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5),
    legend.position = "bottom",
    legend.title = element_text(face = "bold"),
    legend.text = element_text(size = 12),
    legend.spacing.x = unit(0.5, 'cm'),
    legend.background = element_rect(fill = "white", color = "black"),
    panel.grid.major = element_line(color = "grey", linetype = "dashed"),
    panel.grid.minor = element_blank()
  ) +
  scale_color_manual(
    values = c(
      "TBATS" = "blue",
      "ETS" = "red",
      "STL" = "green",
      "NNETAR" = "purple",
      "Combined" = "black",
      "Actual" = "orange"
    )
  ) +
  guides(
    color = guide_legend(
      override.aes = list(linewidth = 2)
    )
  )

— 模型评估 —

# 计算训练集误差
tbats_rmse_train <- sqrt(mean((fitted(tbats_model) - train_data)^2))
ets_rmse_train <- sqrt(mean((fitted(ets_model) - train_data)^2))
stlm_rmse_train <- sqrt(mean((fitted(stlm_model) - train_data)^2))
nnetar_rmse_train <- sqrt(mean((fitted(nnetar_model) - train_data)^2,na.rm = T))
xgb_rmse_train <- sqrt(mean((predict(xgb_model, newdata = train_matrix) - train_data)^2))

# 计算测试集误差
tbats_rmse_test <- sqrt(mean((tbats_forecast$mean - test_data)^2))
ets_rmse_test <- sqrt(mean((ets_forecast$mean - test_data)^2))
stlm_rmse_test <- sqrt(mean((stlm_forecast$mean - test_data)^2))
nnetar_rmse_test <- sqrt(mean((nnetar_forecast$mean - test_data)^2))
xgb_rmse_test <- sqrt(mean((predict(xgb_model, newdata = test_matrix) - test_data)^2))

# 创建数据框
error_comparison <- data.frame(
  Model = c("TBATS", "ETS", "STL", "NNETAR", "Combined"),
  RMSE_Train = c(tbats_rmse_train, ets_rmse_train, stlm_rmse_train, nnetar_rmse_train, xgb_rmse_train),
  RMSE_Test = c(tbats_rmse_test, ets_rmse_test, stlm_rmse_test, nnetar_rmse_test, xgb_rmse_test)
)

# 使用 stargazer 输出三线表
stargazer(
  error_comparison, 
  type = "text",   
  summary = FALSE,  
  rownames = FALSE, 
  title = "Comparison of RMSE for Different Models", 
  header = FALSE    
)

Comparison of RMSE for Different Models
=============================
Model    RMSE_Train RMSE_Test
-----------------------------
TBATS      0.239      0.436  
ETS        0.203      0.247  
STL        0.223      0.505  
NNETAR     0.227      0.194  
Combined   0.119      0.190  
-----------------------------

结果分析

从误差对比结果可以看出，各模型在训练集和测试集上的表现存在显著差异。其中，Combined模型 在训练集（RMSE = 0.119）和测试集（RMSE =0.190）上均表现最优，误差最低，说明其在平衡拟合能力和泛化能力方面表现较好，可能是由于该模型综合了各种时间序列模型中的提取的原始数据的有效信息，然后进行学习，以致于其拟合和预测效果优异。相比之下，TBATS 和 STL 模型在测试集上的误差较高（RMSE 分别为 0.436 和0.505），表明这两个模型可能存在过拟合问题，以致于对测试集数据的适应性不足。而NNETAR 模型在测试集上的误差相对较低（RMSE = 0.194），仅次于 Combined模型，显示其在捕捉数据非线性特征方面具有较强的能力，但是预测能力相对于Combined模型 仍然较弱，因此最终选择 Combined 模型 作为最终模型对未来乙肝发病率进行预测。

多个时间序列与神经网络综合

模型介绍

研究思路

建模过程

结果分析

方康

postgraduate