使用R实现ETS模型, 拟合时间序列, 评估拟合效果。并结合bootstrap进行预测和区间估计
指数平滑模型(Exponential Smoothing Models)是一类广泛使用的时间序列预测方法, 其核心思想是通过对过去数据赋予指数递减的权重,来生成对未来的预测。
指数平滑模型根据数据的特性(如是否存在趋势和季节性)可以分为不同的类型: 简单指数平滑(Simple Exponential Smoothing):适用于没有明显趋势和季节性的平稳时间序列。 霍尔特线性趋势模型(Holt’s Linear Trend Model):扩展了简单指数平滑,适用于具有线性趋势但无季节性的时间序列。 霍尔特-温特斯季节性模型(Holt-Winters Seasonal Model):适用于同时具有趋势和季节性的时间序列。包括加法型和乘法型两种季节性模型。
其中ETS模型是指数平滑方法的一个现代化、统一和扩展的框架,由 Hyndman 等人 提出。它不仅包含了传统的指数平滑方法(如霍尔特-温特斯模型),还通过状态空间表示,提供了一个系统化的方式来表示和选择不同的误差(Error)、趋势(Trend)和季节性(Seasonality)成分的组合。 E(Error):误差成分,可以是加法(A)或乘法(M)。 T(Trend):趋势成分,可以是无(N)、加法(A)或乘法(M)。 S(Seasonality):季节性成分,可以是无(N)、加法(A)或乘法(M)
通过不同的 E、T、S 组合,ETS 模型能够涵盖所有传统的指数平滑方法。使得不同类型的指数平滑方法可以在同一模型框架下进行比较和选择。
在 R 语言的 forecast 包中,使用 ets() 函数时,则模型会自动选择最优的 E、T、S 组合,基于信息准则(如 AIC、BIC)进行选择。 模型处理比较灵活,能够自动识别并选择最适合数据特性的误差、趋势和季节性成分组合。
Bootstrap方法通过对原始样本进行重复抽样来估计分布,进而进行推断。这种方法特别适用于无法直接获取分布的情形。 Bootstrap是一种非参数方法,因为它不依赖于对数据的分布假设。过计算Bootstrap样本的统计量的分布,可以构建原始统计量的置信区间。
#加载包
library(ggplot2)
library(forecast)
library(tidyr)
library(foreach)
library(doParallel)
# 读取文件
data <- read.csv("sxdata.csv")
# 将字符串转换为日期类型
data$t <- as.Date(data$t) # 假设't'列包含日期字符串
#定义时间范围
start_date <- as.Date("2000-01-01")
end_date <- as.Date("2014-12-01")
validation_start <- as.Date("2015-01-01")
validation_end <- as.Date("2017-06-01")
# 筛选训练集和验证集
training_data <- subset(data, t >= start_date & t <= end_date)
validation_data <- subset(data, t >= validation_start & t <= validation_end)
# 将数据转换为时间序列
ts_data <- ts(training_data$value, frequency = 12, start = c(2000, 1))
test_data <- ts(validation_data$value, frequency = 12, start = c(2015, 1))
# 查看趋势
plot.ts(ts_data, main = " Time Series", ylab = "Value", xlab = "Time")
# 使用 STL 进行季节性和趋势分解
stl_decomposition <- stl(ts_data, s.window = "periodic")
plot(stl_decomposition)
# ETS 模型拟合
ets_model <- ets(ts_data)
ets_fitted <- fitted(ets_model)
# 使用 ETS 模型进行预测
forecast_horizon <- 30 # 预测步长
ets_forecast <- forecast(ets_model, h = forecast_horizon)
# 提取预测值(均值)
ets_forecast_values <- ets_forecast$mean
# 绘制预测结果
autoplot(ets_forecast) +
autolayer(fitted(ets_forecast), series = "Fitted") +
labs(title = "ETS Model Forecast", x = "Time", y = "Value") +
theme_minimal()
##bootstrap
# 定义Bootstrap参数
n_bootstrap <- 1000 # 生成1000个Bootstrap样本
# 设置并行环境
num_cores <- parallel::detectCores() - 1 # 根据CPU核心数调整
cl <- makeCluster(num_cores)
registerDoParallel(cl)
# 提取ETS模型的残差和拟合值
ets_residuals <- residuals(ets_model)
ets_fitted_values <- ets_fitted
# 开始并行Bootstrap过程
bootstrap_predictions <- foreach(i = 1:n_bootstrap, .combine = rbind, .packages = "forecast") %dopar% {
# 从残差中进行有放回抽样
bootstrap_residuals <- sample(ets_residuals, size = length(ets_residuals), replace = TRUE)
# 生成新的时间序列数据
new_ts_data <- ets_fitted_values + bootstrap_residuals
# 重新拟合ETS模型
new_ets_model <- ets(new_ts_data)
# 进行预测
new_forecast <- forecast(new_ets_model, h = forecast_horizon)
# 返回预测均值
as.numeric(new_forecast$mean)
}
# 关闭并行集群
stopCluster(cl)
# 计算预测的95%置信区间(使用百分位数)
lower_bound <- apply(bootstrap_predictions, 2, quantile, probs = 0.025)
upper_bound <- apply(bootstrap_predictions, 2, quantile, probs = 0.975)
# 计算预测的中位数(避免极端值)
median_prediction <- apply(bootstrap_predictions, 2, median)
# 计算预测误差
errors2 <- validation_data$value[1:forecast_horizon] - median_prediction
# 计算评价指标
MAE2 <- mean(abs(errors2)) # 平均绝对误差
MAPE2 <- mean(abs(errors2 / validation_data$value[1:forecast_horizon])) * 100 # 平均绝对百分比误差
# 输出评价结果
print(paste("MAE2:", round(MAE2, 3)))
[1] "MAE2: 3.883"
print(paste("MAPE2:", round(MAPE2, 5), "%"))
[1] "MAPE2: 3.22982 %"
# 准备可视化数据
# 创建预测数据框
prediction_data <- data.frame(
t = time(ets_forecast$mean),
lower = lower_bound,
upper = upper_bound,
mean = median_prediction
)
# 创建验证数据框用于绘图
validation_plot_data <- data.frame(
t = time(test_data),
cases = as.numeric(test_data)
)
# 绘制实际值、预测均值及置信区间
ggplot() +
# 绘制实际值
geom_line(data = validation_plot_data, aes(x = t, y = cases, color = "Actual Data"), size = 1) +
# 绘制预测中位数
geom_line(data = prediction_data, aes(x = t, y = mean, color = "Predicted Median"), linetype = "dashed", size = 1) +
# 绘制95%置信区间
geom_ribbon(data = prediction_data, aes(x = t, ymin = lower, ymax = upper, fill = "95% CI"), alpha = 0.3) +
# 添加标签和主题
labs(title = "ETS Model Forecast with Bootstrap Confidence Intervals",
x = "Time",
y = "Value",
color = "Legend", # 图例标题
fill = "Legend") + # 图例标题
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5)) + # 居中标题
scale_color_manual(values = c("Actual Data" = "blue", "Predicted Median" = "magenta")) + # 设置线条颜色
scale_fill_manual(values = c("95% CI" = "green")) # 设置填充颜色
# 定义评价函数
evaluate_model <- function(actual, fitted_values, predicted, actual_test) {
# 计算训练集 R²
ss_res <- sum((actual - fitted_values)^2) # 残差平方和
ss_tot <- sum((actual - mean(actual))^2) # 总平方和
r_squared <- 1 - (ss_res / ss_tot)
# 计算训练集 MSE 和 MAE
mse <- mean((actual - fitted_values)^2)
mae <- mean(abs(actual - fitted_values))
# 验证集 R²
ss_res_test <- sum((actual_test - predicted)^2)
ss_tot_test <- sum((actual_test - mean(actual_test))^2)
r_squared_test <- 1 - (ss_res_test / ss_tot_test)
# 验证集 MSE 和 MAE
mse_test <- mean((actual_test - predicted)^2)
mae_test <- mean(abs(actual_test - predicted))
return(data.frame(
Metric = c("R-squared", "MSE", "MAE"),
Train = c(round(r_squared, 4), round(mse, 4), round(mae, 4)),
Test = c(round(r_squared_test, 4), round(mse_test, 4), round(mae_test, 4))
))
}
# 模型评价
ets_results <- evaluate_model(ts_data, ets_fitted, median_prediction, test_data)
print(ets_model) # 打印ETS模型参数
ETS(A,A,A)
Call:
ets(y = ts_data)
Smoothing parameters:
alpha = 0.0183
beta = 1e-04
gamma = 1e-04
Initial states:
l = 22.6894
b = 0.4932
s = -2.411 -10.0215 -17.4241 -20.4455 -15.5605 -8.7317
0.4715 9.4264 18.3544 19.2057 16.7531 10.3832
sigma: 4.9472
AIC AICc BIC
1527.549 1531.327 1581.830
print(ets_results) # 打印模型评价指标
Metric Train Test
1 R-squared 0.9737 0.8689
2 MSE 22.2989 26.7532
3 MAE 3.7115 3.8828
checkresiduals(ets_model)
Ljung-Box test
data: Residuals from ETS(A,A,A)
Q* = 22.675, df = 24, p-value = 0.5391
Model df: 0. Total lags used: 24
根据图1:观察原始数据趋势可知时间序列从2000年到2015年整体呈现上升趋势,表明数据中存在一个明确的长期增长趋势。 并且有明显的周期性波动(季节性)
因此,采用ETS建模是适合的,尤其是需要选择包含趋势(T)和季节性(S)成分的ETS模型, 例如 ETS(A,A,A)(加性趋势和加性季节性)或 ETS(M,A,M)(乘性趋势和乘性季节性)
STL分解之前的时间序列图一致,显示出数据具有明显的趋势性和周期性(季节性波动)
根据拟合结果图3可以看到ETS模型较好地捕捉了时间序列中的趋势和季节性特征,并且验证数据集的实际值在预测置信区间内,说明模型在验证集上的表现是可靠的。
加入重抽样后结果如图4
训练集上: R² = 0.9737,表明模型在训练集上解释了97%以上的时间序列变化,拟合效果较好。
MSE = 22.30,MAE = 3.71,误差较低,说明模型对训练数据的拟合精度较高。
验证集上: R² = 0.8696,模型在验证集上的表现稍逊于训练集,但达到了86%以上的解释能力,表现较为优秀。
MSE = 26.60,MAE = 3.87,与训练集误差接近,表明模型没有明显的过拟合。
验证集的MAPE为3.22%,这表明预测误差占实际值的比例很低,模型预测精度较高。
残差检验结果和图表中显示: 残差接近白噪声,无明显的模式或自相关性,表明模型已经很好地捕捉了数据中的主要趋势和季节性成分。
ETS(A,A,A)模型中 alpha = 0.0183: 平滑水平参数接近0表明模型对最近观测值的响应非常缓慢,主要依赖于历史数据。
beta = 1e-04: 趋势平滑参数接近0,意味着趋势变化非常缓慢。
gamma = 1e-04: 季节性平滑参数接近0,表明季节性成分更新非常缓慢。
平滑参数较低可能是因为随机数据非常稳定且变化不大,导致模型不能很好的响应新数据,可以手动调参,但不适当的参数设置可能导致模型性能下降。调整后需重新评估模型的表现。