基于ARIMA与虚拟变量的预测：疫情虚拟变量引入报错排查

问题背景 #

处理多伦多2014-2021年入室盗窃月度时间序列数据时，想要引入2020年3月疫情起始的虚拟变量，构建带ARIMA误差的回归模型。原auto.arima()输出的ARIMA(1,0,1)模型未考虑疫情影响，仅基于序列均值回归，但创建虚拟变量时出现类型不兼容报错。

报错信息 #

In ifelse(time(BEDATA_GROUPEDtsssarima) >= yearmonth("2020-03"), :
Incompatible methods ("Ops.ts", ">=.vctrs_vctr") for ">="

原代码 #

# Create a binary time series that indicates the start of the pandemic
library(fpp3)
library(forecast)
library(zoo)

# Check if timeseries
class(BEDATA_GROUPED)

#Convert timeseries
BEDATA_GROUPEDtsssarima <- ts(BEDATA_GROUPED[,2], frequency = 12, start = c(2014, 1))
class(BEDATA_GROUPEDtsssarima)

#Plot
forecast::autoplot(BEDATA_GROUPEDtsssarima)

# Assume that the pandemic began in March 2020
pandemic_dummy <- ifelse(time(BEDATA_GROUPEDtsssarima) >= yearmonth("2020-03"), 1, 0)

# Use auto.arima() to fit an ARIMA model with the dummy variable as an exogenous variable
beddatamodel <- auto.arima(BEDATA_GROUPEDtsssarima, xreg = pandemic_dummy, ic="aic", trace = TRUE)

# Create a binary time series that indicates the start of the pandemic
# In this example, we will assume that the pandemic began in March 2020
pandemic_dummy <- ifelse(time(BEDATA_GROUPEDtsssarima) >= yearmonth("2020-03"), 1, 0)

# Use auto.arima() to fit an ARIMA model with the dummy variable as an exogenous variable
beddatamodel <- auto.arima(BEDATA_GROUPEDtsssarima, xreg = pandemic_dummy, ic="aic", trace = TRUE)

# Create a binary time series for the forecast period that includes the pandemic dummy variable
forecast_period <- time(BEDATA_GROUPEDtsssarima)["2022/01/01/":"2023/12/31/"]
pandemic_dummy_forecast <- ifelse(forecast_period >= yearmonth("2020-03"), 1, 0)

# Use the forecast()
forecast(pandemic_dummy_forecast)

数据集 #

structure(list(occurrence_yrmn = c("2014-January", "2014-February", 
"2014-March", "2014-April", "2014-May", "2014-June", "2014-July", 
"2014-August", "2014-September", "2014-October", "2014-November", 
"2014-December", "2015-January", "2015-February", "2015-March", 
"2015-April", "2015-May", "2015-June", "2015-July", "2015-August", 
"2015-September", "2015-October", "2015-November", "2015-December", 
"2016-January", "2016-February", "2016-March", "2016-April", 
"2016-May", "2016-June", "2016-July", "2016-August", "2016-September", 
"2016-October", "2016-November", "2016-December", "2017-January", 
"2017-February", "2017-March", "2017-April", "2017-May", "2017-June", 
"2017-July", "2017-August", "2017-September", "2017-October", 
"2017-November", "2017-December", "2018-January", "2018-February", 
"2018-March", "2018-April", "2018-May", "2018-June", "2018-July", 
"2018-August", "2018-September", "2018-October", "2018-November", 
"2018-December", "2019-January", "2019-February", "2019-March", 
"2019-April", "2019-May", "2019-June", "2019-July", "2019-August", 
"2019-September", "2019-October", "2019-November", "2019-December", 
"2020-January", "2020-February", "2020-March", "2020-April", 
"2020-May", "2020-June", "2020-July", "2020-August", "2020-September", 
"2020-October", "2020-November", "2020-December", "2021-January", 
"2021-February", "2021-March", "2021-April", "2021-May", "2021-June", 
"2021-July", "2021-August", "2021-September", "2021-October", 
"2021-November", "2021-December"), MCI = c(586, 482, 567, 626, 
625, 610, 576, 634, 636, 663, 657, 556, 513, 415, 510, 542, 549, 
618, 623, 666, 641, 632, 593, 617, 541, 523, 504, 536, 498, 552, 
522, 519, 496, 541, 602, 570, 571, 492, 560, 525, 507, 523, 593, 
623, 578, 657, 683, 588, 664, 582, 619, 512, 630, 644, 563, 654, 
635, 732, 639, 748, 719, 567, 607, 746, 739, 686, 805, 762, 696, 
777, 755, 675, 704, 617, 732, 609, 464, 487, 565, 609, 513, 533, 
505, 578, 526, 418, 428, 421, 502, 452, 509, 492, 478, 469, 457, 
457)), class = c("tbl_df", "tbl", "data.frame"), row.names = c(NA, 
-96L))

报错原因与解决方法 #

报错核心是类型不匹配：time(BEDATA_GROUPEDtsssarima)返回ts类型时间值，yearmonth("2020-03")返回vctrs_vctr类型对象，两者无法直接比较。以下是三种可行的解决方式：

方法1：统一为数值型时间格式 #

将ts时间转换为YYYY.MM格式的数值（如2020年3月对应2020.25），再与目标月份的数值比较：

# 转换目标月份为数值格式
pandemic_start <- 2020 + (3-1)/12
pandemic_dummy <- ifelse(time(BEDATA_GROUPEDtsssarima) >= pandemic_start, 1, 0)

方法2：用zoo包统一时间类型 #

借助zoo包的as.yearmon将两种时间转换为同一种类型：

# 将ts时间转换为yearmon类型
ts_time <- as.yearmon(time(BEDATA_GROUPEDtsssarima))
# 目标时间也用as.yearmon创建
pandemic_start <- as.yearmon("2020-03")
pandemic_dummy <- ifelse(ts_time >= pandemic_start, 1, 0)

方法3：按数据行位置生成虚拟变量 #

已知数据共96行，2020年3月是第75行（2014-2019共72行，2020年1-2月为73、74行），直接按位置赋值：

pandemic_dummy <- rep(0, length(BEDATA_GROUPEDtsssarima))
pandemic_dummy[75:96] <- 1

修正后的完整建模代码 #

library(fpp3)
library(forecast)
library(zoo)

# 加载数据集（假设已导入BEDATA_GROUPED）
# BEDATA_GROUPED <- structure(...)

# 转换为ts对象
BEDATA_GROUPEDtsssarima <- ts(BEDATA_GROUPED[,2], frequency = 12, start = c(2014, 1))

# 用方法2创建虚拟变量
ts_time <- as.yearmon(time(BEDATA_GROUPEDtsssarima))
pandemic_start <- as.yearmon("2020-03")
pandemic_dummy <- ifelse(ts_time >= pandemic_start, 1, 0)

# 拟合带外生变量的ARIMA模型
beddatamodel <- auto.arima(BEDATA_GROUPEDtsssarima, xreg = pandemic_dummy, ic="aic", trace = TRUE)

# 生成预测期虚拟变量：2022-2023全处于疫情期
forecast_period <- seq(as.yearmon("2022-01"), as.yearmon("2023-12"), by = 1/12)
pandemic_dummy_forecast <- rep(1, length(forecast_period))

# 生成并查看预测
bed_forecast <- forecast(beddatamodel, xreg = pandemic_dummy_forecast)
print(bed_forecast)
forecast::autoplot(bed_forecast)

额外说明 #

• 原代码中重复创建虚拟变量和拟合模型，建议删除冗余代码。
• 预测期时间不要依赖原ts对象的time（原数据仅到2021年12月），直接用seq生成目标区间的月度时间。
• 拟合带外生变量的ARIMA模型后，模型会自动纳入疫情的影响，可对比原模型结果观察疫情对入室盗窃数量的显著作用。

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