ARIMA模型

使用R实现ARIMA模型， 拟合时间序列， 评估拟合效果。并结合bootstrap进行预测和区间估计

吴嘉栋

日期

基本介绍

ARIMA（AutoRegressive Integrated Moving Average，自动回归积分滑动平均）模型是一种常用于时间序列预测的统计模型。ARIMA模型通过三个主要参数来描述时间序列数据的特征：AR（自回归）、I（差分）和MA（滑动平均）。该模型假设时间序列数据是平稳的或可以通过差分使其平稳，并通过自回归和滑动平均成分来捕捉数据中的依赖关系。ARIMA模型广泛应用于金融、经济、气候等领域的时间序列预测。

为了提高ARIMA模型的预测精度，结合了Bootstrap方法，通过重复抽样生成样本的分布，进行区间估计，从而获得更准确的置信区间。

下载讲解视频

ARIMA

AR（自回归）：通过历史数据的线性组合来预测当前值，即当前时间点的数据依赖于其前几个时间点的值。

I（差分）：用于处理非平稳时间序列，使其变得平稳，从而便于建模。差分操作是通过计算序列当前值与前一个值的差异来实现的。

MA（滑动平均）：考虑到预测误差的影响，MA部分通过历史误差的线性组合来预测当前值。

ARIMA模型适用于平稳时间序列，若数据非平稳，通过差分（I部分）将其转换为平稳序列。该模型的核心思想是结合自回归（AR）和滑动平均（MA）过程来捕捉数据中的依赖性，并通过差分来解决非平稳问题。

更详细内容以及模型公式可见https://blog.csdn.net/weixin_62586597/article/details/120626460

Bootstrap

Bootstrap方法是一种非参数统计方法，主要用于通过重复抽样来估计样本统计量的分布，特别适用于数据分布未知或不满足常规统计假设的情况。在ARIMA模型中，Bootstrap方法常用于以下目的：

估计预测区间：通过生成多组样本的预测值，从而得到预测值的置信区间。

提高模型的稳健性：当ARIMA模型的误差项不服从正态分布时，Bootstrap方法可以帮助获得更可靠的参数估计。

ARIMA+Bootstrap

在时间序列预测中，ARIMA模型能够有效地捕捉数据中的趋势和季节性变化。但ARIMA模型的参数估计往往依赖于假设模型误差服从正态分布，这可能在某些情况下不成立。为了解决这个问题，我们可以使用Bootstrap方法来进行区间估计，通过对历史数据进行重复抽样，模拟出多组时间序列的样本分布，从而估计预测值的置信区间。

具体而言，首先使用ARIMA模型对数据进行拟合，得到残差。然后，从残差中随机抽取样本，生成新的时间序列，并对每个样本进行ARIMA模型拟合，反复进行多次（例如1000次）。通过这种方法，我们得到多个预测值，从而形成预测值的分布。根据这些分布，我们可以计算预测值的置信区间，通常取2.5%和97.5%分位数作为95%置信区间。这一过程能够为ARIMA模型的预测提供更为可靠和精准的区间估计，避免了误差项非正态分布对预测结果的影响，增强了模型的稳健性和精确性。

代码示例

1、数据导入

#加载包（运行前记得先安装）
library(tseries)
library(forecast)
library(ggplot2)
library(dplyr)
data <- data.frame(
  t = c(1:96), 
  cases = c(
1468, 1267, 4922, 12190, 13328, 11939, 7578, 3152, 3307, 2676, 2666, 2331, 1417, 470, 1991, 5156, 10683, 11478, 7413, 4230, 4118, 4619, 7046, 6938, 3827, 3742, 7623, 15232, 19672, 10486, 4848, 1778, 2030, 2552, 5281, 6737, 
3739, 3839, 8506, 10209,8982, 9265, 6752, 3427, 3991, 3347, 3890, 5776, 
3366, 3439, 11286, 21903, 19796, 14231, 8298, 4955, 7596, 8878, 6825, 5214, 
3707, 2967, 5131, 12216, 15203, 10605, 5455, 2985, 2782, 3050, 4152, 5170, 
5383, 3089, 6770, 12533, 16241, 16154, 12090, 7863,6859, 8564, 8886, 5730,
4301, 5057, 5889, 8954, 9947, 9032, 6614,4834, 6172, 7658, 6813, 4970)
) #这里使用处理过的示例数据，通常要导入数据文件

2、数据处理

#将数据被分为两部分：训练和验证
training_data<-subset(data,t >= 1 & t <= 84) #假设是2000-2006
validation_data<-subset(data,t >= 85 & t <= 96) #假设是2007

# 将数据转化为时间序列,使用ts()函数
ts_data <- ts(training_data$cases, frequency =12, start = c(2000, 1))
test_data <- ts(validation_data$cases, frequency = 12, start = c(2007, 1))

# 拆分季节性、趋势、随机性(时间序列分解，看看各种特性)
decomposition_data <- decompose(ts_data)
plot(decomposition_data)

plot（decomposition_data）可以对数据的流行曲线和特征有一个初步的观察

3、数据检查

#稳定性检验，P值小于0.01说明数据稳定，可以建模
adf.test(ts_data)
#白噪声检验，P值小于0.05，序列不是白噪声，可以建模
Box.test(ts_data,type = "Ljung-Box")

adf.test p-value <0.01

Box.test p-value = 3.13e-11

平稳、非白噪声，可建模

4、模型构建

#使用auto.arima，自动选取最优模型
auto.arima(ts_data, 
           stepwise = TRUE,  # 是否启用逐步搜索算法，默认启用。逐步搜索会自动选择最优的p、d、q值
           trace = TRUE,     # 是否输出每一步的模型选择过程。启用后会显示每一步选择的模型
           stationary = TRUE,  # 指定数据是否为平稳数据。TRUE表示将自动对数据进行平稳性检验并进行差分（如有必要）
           ic = c("aic"),    # 使用AIC准则选择最佳模型，AIC即赤池信息量准则。可以设置为AIC或BIC，这里选择了AIC
           max.p = 5,        # 设置p的最大值。默认最大值是5，表示AR部分的最大阶数为5。你可以根据数据复杂性调整
           max.q = 5,        # 设置q的最大值。默认最大值是5，表示MA部分的最大阶数为5。你可以根据数据复杂性调整
           max.d = 2,        # 设置d的最大值。最大差分次数，通常为1或2
           seasonal = TRUE,  # 是否进行季节性差分
)
#运行上方代码会得到不同公式下，模型的AIC值。AIC值较小表示模型更合适，一般选择最小的AIC所代表的模型
#找到上方代码结果中显示的best model(这个数据的结果显示，pdq、PDQ取值为0、0、2以及2、0、0)
fit <- arima(ts_data, 
             order = c(0, 0, 2),  # ARIMA模型的p=0, d=0, q=2，表示0阶自回归和2阶移动平均项，无差分
             seasonal = list(order = c(2, 0, 0), period = 12)  # 季节性部分的阶数为2，0，0，表示有2阶季节性AR项，period=12表示周期为12（月度数据时，周期为12）
)

除了使用代码自动选择，也可以通过画出数据的ACF、PACF图来进行手动选择 非季节性部分： p（AR阶数）：通过PACF图来确定。若PACF在滞后p处显著衰减且在后续滞后消失，选择p为此滞后值。 d（差分阶数）：若数据非平稳，通常进行差分，差分次数即为d。 q（MA阶数）：通过ACF图来确定。若ACF在滞后q处显著衰减且后续滞后消失，选择q为此滞后值。 季节性部分（可以使用decompose将季节性部分拆解出来）： P（季节性AR阶数）：通过季节性PACF图来确定。若季节性PACF在滞后s（季节周期）时显著衰减，选择P=1或P=2。 D（季节性差分阶数）：若数据非平稳，通常进行差分，差分次数即为D。 Q（季节性MA阶数）：通过季节性ACF图来确定，若季节性滞后期的ACF显著，选择Q=1或Q=2。 更详细的讲解可见：https://blog.csdn.net/weixin_62586597/article/details/120626460

5、模型检查

#正态性检验
shapiro.test(fit$residuals)#残差是否服从正态(0假设残差正态)
qqnorm(fit$residuals, main = "Q-Q Plot of Residuals")#QQ图
#自相关性检验(0假设残差是白噪声)
Box.test(fit$residuals,lag=24,type="Ljung-Box") #非季节性=10,季节性=period*2
tsdisplay(residuals(fit),lag.max=20,main='残差')#残差相关性

shapiro.test p-value = 0.06676

Box.test p-value = 0.8147

残差正态且独立，模型效果良好

6、效果评估

#预测测试集2007年
prediction<-forecast(fit,12)#向后预测12个月，也就是测试集的年份

# 提取预测值和实际值
predicted_values <- prediction$mean  # 对验证集的预测值
actual_values <- test_data  # 验证集的实际值

# 计算残差（实际-预测）
errors <- actual_values - predicted_values

# 计算指标
mean(abs(errors))  # 平均绝对误差（通常情况下越小越好）
mean(abs(errors / actual_values)) * 100

MAE: 1203.096
MAPE: 17.3117 %

mape<20%,精确度不错

7、单次预测

#预测2008年
prediction2<-forecast(fit,24)
plot(prediction2)

将步长延长，就可以获得预测结果。可以看到其实预测部分已经有区间。但值得注意的是，预测的偏保守而不够集中，此时有请重抽样。

8、加入抽样

# 获取训练集的拟合值和残差序列
fitted_values <- fitted(fit)
residuals_values <- residuals(fit)

# 定义重抽样次数和预测长度
n_bootstrap <- 1000  # 生成1000个Bootstrap样本
forecast_horizon <- 12  # 预测步长12个月

# 建立存储每个样本预测值的矩阵
bootstrap_predictions <- matrix(NA, nrow = n_bootstrap, ncol = forecast_horizon)

set.seed(123)  # 复现性，种子

#循环1000次
for (i in 1:n_bootstrap) {
  # 从残差中进行有放回抽样
  bootstrap_residuals <- sample(residuals_values, size = length(residuals_values), replace = TRUE)
  
  # 将抽样的残差与拟合值相加，得到新的原始序列
  new_ts_data <- fitted_values + bootstrap_residuals
  
  # 重新拟合ARIMA模型
  new_fit <- arima(new_ts_data, order = c(0, 0, 2), seasonal = list(order = c(2, 0, 0), period = 12))
  
  # 对新的序列进行预测
  bootstrap_prediction <- forecast(new_fit, h = forecast_horizon)
  
  # 存储每次的预测值
  bootstrap_predictions[i, ] <- bootstrap_prediction$mean 
}

#现在已经获得了1000个预测路径
# 计算预测的95%置信区间（使用百分位数）
lower_bound <- apply(bootstrap_predictions, 2, function(x) quantile(x, 0.025))
upper_bound <- apply(bootstrap_predictions, 2, function(x) quantile(x, 0.975))

# 计算预测的中位数（避免极端值）
median_prediction <- apply(bootstrap_predictions, 2, median)

#看看重抽样的结果效果如何
errors2 <- validation_data$cases - median_prediction

# 计算指标
mean(abs(errors2))  # 平均绝对误差
mean(abs(errors2 / validation_data$cases)) * 100  # # 输出结果

MAE2: 1063.391

MAPE2: 17.21197 %

9、来个简图

# 把重抽样结果转为数据框
prediction_data <- data.frame(
  t = 85:96,  # 预测的时间段
  lower = lower_bound, 
  upper = upper_bound, 
  mean = median_prediction
)

ggplot() +
  # 绘制实际值
  geom_line(data = validation_data, aes(x = t, y = cases, color = "Actual Data"), linewidth = 1) +
  # 绘制预测均值
  geom_line(data = prediction_data, aes(x = t, y = mean, color = "Predicted Mean"), linetype = "dashed", linewidth = 1) +
  # 绘制95%置信区间
  geom_ribbon(data = prediction_data, aes(x = t, ymin = lower, ymax = upper, fill = "95% CI"), alpha = 0.3) +
  labs(title = "Report vs Simulate", 
       x = "Time (t)", 
       y = "Cases", 
       color = " ",  # 图例标题
       fill = " ") +  # 图例标题
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5)) +  # 居中标题
  scale_color_manual(values = c("Actual Data" = "blue", "Predicted Mean" = "red")) +  # 设置线条颜色
  scale_fill_manual(values = c("95% CI" = "lightblue"))  # 设置填充颜色

分析讨论

  1. 在数据平稳、合适的情况下，ARIMA作为经典的时间序列模型，预测的效果并不会差（实际情况中，很多数据使用arima是不够的。这个模型不够复杂，抓特征的能力并没有那么强，特别是非线性部分）。R包自带的预测区间，是基于单次的残差序列来运算的，随着预测步长的前进，误差会越来越大，就导致预测区间的宽度比较宽，而且会变大，这意味着预测的非常保守。而使用重抽样的方法，通过生成多个预测路径，就能够得到更为精准的预测区间。

  2. 然而，ARIMA+Bootstrap方法仍然有其局限性。首先，ARIMA模型自身无法捕捉非线性特征和长期依赖，尽管Bootstrap通过重抽样弥补了ARIMA的短期预测能力，但如果残差存在系统性误差，重抽样的效果也会受到限制。通俗来讲，arima本身预测的很差，那bootstrap再怎么弥补也没太大意义。此外，Bootstrap方法本身的计算复杂度较高，尤其在大数据集上，其抽样和计算过程需要较长时间，会影响了效率和实时性。

  3. 如何改进？首先，针对真实数据，在拟合的时候选择最合适的模型。尽可能的捕捉数据本身的复杂趋势，在初步拟合效果不错的情况下，重抽样可以达到锦上添花的效果；其次，可以优化Bootstrap方法，采用动态抽样策略以适应时间序列的季节性波动和变化。同时，为了提升计算效率，可以借助GPU加速等技术提升在大规模数据上的实时预测能力；最后，可以使用残差建模的方法，运用深度学习在捕捉复杂关系的优势，进一步拟合残差，使用组合模型来提升预测精度。这些改进将帮助模型在处理复杂、非线性、高噪声数据时，提供更精准、可靠和高效的预测。

吴嘉栋
Still KD