R中predict()无法识别sarima对象的问题求助

问题场景 #

我正在针对汽油需求做时间序列分析，采用TSDL包中的销售数据，步骤如下：

1. 加载并可视化数据：

library(tsdl)
options(warn=-1)
tsdl_gasdemand <- subset(tsdl, "Sales")
X_t <- tsdl_gasdemand[[1]]
plot(X_t, main = "Monthly gasoline demand Ontario gallon millions 1960 – 1975", ylab = "")

2. 完成Box-Cox变换与差分处理后，通过ACF/PACF图选定SARIMA(3,0,0)x(1,1,0)₁₂模型，用astsa包的sarima()拟合：

model4 <- sarima(xdata=train.bc, p=3, d=0, q=0, P=1, D=1, Q=0, S=12, details = F)

拟合结果包含constant项：

Coefficients:
          ar1     ar2     ar3     sar1  constant
      -0.2432  0.0356  0.5622  -0.4113    1.6712
s.e.   0.0692  0.0735  0.0734   0.0812    0.0970

但调用predict()时触发错误：

pred.demand <- predict(model4$fit, n.ahead = 12)
# Error in eval(expr, p) : object 'constant' not found

尝试用base R的arima()拟合同款模型：

SARIMA <- arima(train.bc, order=c(3,0,0), method="ML",
      seasonal = list(order = c(1,1,0), period = 12))

但残差PACF在滞后5阶处存在显著峰值，模型效果不佳。另外astsa::sarima.for()可生成预测图，但无法查看计算过程。

解决方案 #

方案1：修复predict()对constant项的识别 #

astsa::sarima()返回的模型对象中，constant项的存储方式导致base R的predict.Arima()无法识别，推荐用forecast包的forecast()函数替代，它对这类模型的兼容性更好：

library(forecast)
# 生成预测，h指定预测步数
pred.demand <- forecast(model4$fit, h = 12)
# 查看预测结果
print(pred.demand)
# 可视化
plot(pred.demand)

如果坚持用predict()，可手动修正模型对象的调用参数，让constant变量被识别：

# 提取constant系数值
const_coef <- coef(model4$fit)["constant"]
# 修改模型对象的call，添加constant参数
model4$fit$call$constant <- const_coef
# 重新执行预测
pred.demand <- predict(model4$fit, n.ahead = 12)

方案2：为base R arima()添加constant项 #

base R的arima()默认在存在季节性差分（D=1）时不包含常数项，可通过以下两种方式手动添加：

方式1：使用include.mean参数

当模型包含季节性差分D=1时，include.mean=TRUE会为差分后的序列添加均值，对应原模型的constant项：

SARIMA <- arima(train.bc, order = c(3,0,0), method = "ML",
                seasonal = list(order = c(1,1,0), period = 12),
                include.mean = TRUE)

方式2：通过xreg传入常数外生变量

手动构造全1的外生变量，将constant作为回归项加入模型：

# 训练集的常数外生变量
train_xreg <- rep(1, length(train.bc))
# 拟合模型
SARIMA <- arima(train.bc, order = c(3,0,0), method = "ML",
                seasonal = list(order = c(1,1,0), period = 12),
                xreg = train_xreg)
# 预测时需传入对应步数的常数外生变量
pred_xreg <- rep(1, 12)
pred.demand <- predict(SARIMA, n.ahead = 12, xreg = pred_xreg)

若残差仍存在滞后5阶峰值，可尝试调整模型参数（比如加入MA(5)项，即order=c(3,0,1)），或用auto.arima()自动筛选最优模型：

library(forecast)
auto_sarima <- auto.arima(train.bc, seasonal = TRUE, trace = TRUE)
# 生成预测
pred.demand <- forecast(auto_sarima, h = 12)

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