如何从拟合的ARIMA模型模拟100条轨迹并计算平均值？

问题描述 #

我正在解决如下时间序列建模任务：

• 生成随机数据
• 为数据拟合时间序列模型
• 从拟合模型中模拟100条轨迹（每条含100个点）并计算这些轨迹的平均值

已完成步骤 #

1. 模拟数据：使用R的forecast和ggplot2包生成1000个标准正态分布随机数据，转换为频率12的时间序列random_ts
2. 拟合模型：通过auto.arima()得到ARIMA(0,0,0)(0,0,1)[12]模型，参数如下：

Series: random_ts 
ARIMA(0,0,0)(0,0,1)[12] with zero mean 

Coefficients:
         sma1
      -0.0792
s.e.   0.0305

sigma^2 = 0.9771:  log likelihood = -1406.89
AIC=2817.79   AICc=2817.8   BIC=2827.6

遇到的疑问 #

如何从该拟合模型中模拟100条轨迹并计算平均值？是否需要使用arima.sim()函数？原以为可以采用类似预测的方式实现。

尝试代码 #

library(forecast)
library(ggplot2)

set.seed(123)


fit <- auto.arima(random_ts)

simulate_forecast <- function(model, h) {
    sim <- simulate(model, nsim = h)
    return(as.numeric(sim))
}

n_sims <- 100
h <- 100
simulations <- replicate(n_sims, simulate_forecast(fit, h))

avg_prediction <- rowMeans(simulations)

time <- seq_len(n + h)
original_data <- c(as.numeric(random_ts), rep(NA, h))
sim_data <- rbind(matrix(NA, nrow = n, ncol = n_sims), simulations)
avg_data <- c(rep(NA, n), avg_prediction)

plot_data <- data.frame(
    time = rep(time, n_sims + 2),
    value = c(original_data, as.vector(sim_data), avg_data),
    type = factor(rep(c("Original", rep("Simulation", n_sims), "Average"), each = n + h))
)

 ggplot(plot_data, aes(x = time, y = value, group = interaction(type, rep(1:(n_sims+2), each = n + h)), color = type)) +
    geom_line(aes(alpha = type)) +
    scale_color_manual(values = c("Original" = "blue", "Simulation" = "red", "Average" = "black")) +
    scale_alpha_manual(values = c("Original" = 1, "Simulation" = 0.02, "Average" = 1)) +  #
    theme_minimal() +
    labs(title = "Original Data, 100 Simulations, and Average Prediction",
         x = "Time",
         y = "Value",
         color = "Type") +
    theme(legend.position = "bottom")

解决方案 #

你的思路方向是对的，不需要使用arima.sim()，forecast包中的simulate()函数可以直接基于auto.arima()拟合的模型对象生成模拟轨迹，比arima.sim()更便捷，因为它会自动继承模型的所有参数（包括季节性项）。

你的尝试代码存在几个小问题需要修正：

1. 变量n未定义：n应为原始时间序列的长度，即length(random_ts)
2. 模拟轨迹的起始逻辑：simulate()默认从模型拟合的最后一个观测值之后生成新序列，符合需求

修正后的完整代码 #

library(forecast)
library(ggplot2)

# 1. 生成随机数据（补充完整数据生成步骤）
set.seed(123)
random_ts <- ts(rnorm(1000), frequency = 12)

# 2. 拟合模型
fit <- auto.arima(random_ts)

# 3. 模拟100条轨迹，每条100个点
simulate_forecast <- function(model, h) {
  as.numeric(simulate(model, nsim = h))
}

n_sims <- 100
h <- 100
n <- length(random_ts)  # 定义原始数据长度
simulations <- replicate(n_sims, simulate_forecast(fit, h))

# 计算每个时间点的轨迹平均值
avg_prediction <- rowMeans(simulations)

# 整理绘图数据
time <- seq_len(n + h)
original_data <- c(as.numeric(random_ts), rep(NA, h))
sim_data <- rbind(matrix(NA, nrow = n, ncol = n_sims), simulations)
avg_data <- c(rep(NA, n), avg_prediction)

plot_data <- data.frame(
  time = rep(time, n_sims + 2),
  value = c(original_data, as.vector(sim_data), avg_data),
  type = factor(rep(c("Original", rep("Simulation", n_sims), "Average"), each = n + h))
)

# 绘图
ggplot(plot_data, aes(x = time, y = value, group = interaction(type, rep(1:(n_sims+2), each = n + h)), color = type)) +
  geom_line(aes(alpha = type)) +
  scale_color_manual(values = c("Original" = "blue", "Simulation" = "red", "Average" = "black")) +
  scale_alpha_manual(values = c("Original" = 1, "Simulation" = 0.05, "Average" = 1)) +  # 提高模拟轨迹透明度，便于观察
  theme_minimal() +
  labs(title = "原始数据、100条模拟轨迹及平均预测值",
       x = "时间",
       y = "数值",
       color = "类型") +
  theme(legend.position = "bottom")

关键说明 #

• simulate(model, nsim = h)：直接基于拟合好的ARIMA模型生成h个新观测值，自动处理季节性参数（如模型中的sma1项）
• replicate(n_sims, ...)：重复执行模拟函数n_sims次，得到h × n_sims的矩阵，每行对应一个时间点的所有模拟值
• rowMeans(simulations)：计算每个时间点上100条轨迹的平均值，得到平均预测序列

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