基于单变量低频时间序列的最优估计及季度预测方法咨询

假设我仅能获取年度观测数据，需要仅基于单变量时间序列的历史观测值推导当前最优估计值，对应历史数据如下：

我原本的思路是先拟合ARIMA模型，再对ARIMA模型的状态空间表示使用Kalman滤波，推导下一个观测值的最优估计，现有三个问题需要咨询：

1. 我不只想预测下一个年度数据点，还需要推导2021年四个季度的估计值（w、x、y、z），且要求离最后一个真实观测值越远，估计的不确定性越高，该问题是否有可行的解决方法，具体应如何实现？
2. 如何将预测不确定性随观测距离增加而上升的特性纳入估计过程？
3. ARIMA + KF是解决该问题的最适用方案吗？是否有其他更合适的处理方法？

我使用R或Python实现均可，感谢分享相关思路。

数据集的R结构代码如下：

structure(list(date = c("2005-12-31", "2006-12-31", "2007-12-31", 
"2008-12-31", "2009-12-31", "2010-12-31", "2011-12-31", "2012-12-31", 
"2013-12-31", "2014-12-31", "2015-12-31", "2016-12-31", "2017-12-31", 
"2018-12-31", "2019-12-31", "2020-12-31", "2021-03-31", "2021-06-30", 
"2021-09-30", "2021-12-31"), id = c("ABC", "ABC", "ABC", "ABC", 
"ABC", "ABC", "ABC", "ABC", "ABC", "ABC", "ABC", "ABC", "ABC", 
"ABC", "ABC", "ABC", "ABC", "ABC", "ABC", "ABC"), value = c("3150000", 
"5970000", "6640000", "6390000", "7130000", "7270000", "7030000", 
"7360000", "7470000", "7810000", "8690000", "8910000", "2820000", 
"4380000", "2720000", "2480000", "w", "x", "y", "z")), class = "data.frame", row.names = c(NA, 
-20L))

问题1解答 #

该问题完全可解，属于时间拆量（Temporal Disaggregation） 场景，核心实现逻辑如下：

• 第一步预处理：提取2005-2020年的年度实测值，将序列扩展为季度频率，所有非年末的季度值设为缺失值，相当于你只在每4个季度的最后一期有观测值，其余为待补全的缺失值。
• 第二步建模补全：直接基于含缺失值的季度序列拟合状态空间形式的时间序列模型，卡尔曼滤波会自动补全所有缺失的季度值，同时输出每个估计值的标准误。如果要求2021年四个季度的估计值总和等于2021年的年度预测值，可以额外加总和约束，R的tempdisagg包、Python的statsmodels都支持该约束设置。

问题2解答 #

预测不确定性随步长增大而上升是状态空间类时间序列模型的固有特性，不需要额外手动设置。以ARIMA模型为例，h步预测的方差计算公式为$\sigma^2_h = \sigma^2 \sum_{i=0}^{h-1} \psi_i^2$，其中$\psi_i$是模型的脉冲响应系数，天然随h增大而递增。你只需要在输出预测结果时同步输出对应步长的标准误或置信区间即可：

• R中调用predict函数可直接返回每个预测值的标准误；
• Python中statsmodels的get_prediction方法可调用conf_int输出指定置信水平的区间。

问题3解答 #

ARIMA+KF不是唯一方案，也不一定是最优方案，尤其你的数据在2017年出现明显的结构突变，ARIMA对突变的适应性较差，可选的替代方案包括：

• 结构时间序列模型（STS）：可显式定义趋势、季节、随机冲击项，对结构变化的鲁棒性优于普通ARIMA，R的KFAS包、Python的statsmodels都支持实现；
• 贝叶斯结构时间序列模型（BSTS）：小样本下的不确定性估计更稳定，可直接输出参数的后验分布，R的bsts包可快速实现；
• 自助法（Bootstrap）预测：不需要显式构建状态空间模型，通过多次重采样生成预测分布，用分位数表示不确定性，实现门槛更低。

简单实现示例（R） #

library(forecast)
# 构造年度序列
yr_value <- c(3150000,5970000,6640000,6390000,7130000,7270000,7030000,7360000,7470000,7810000,8690000,8910000,2820000,4380000,2720000,2480000)
yr_ts <- ts(yr_value, frequency = 1, start = 2005)
# 扩展为季度序列，非年末设为NA
q_ts <- ts(rep(NA, 16*4 + 4), frequency = 4, start = c(2005,1))
for(i in 1:16) q_ts[4*i] <- yr_ts[i]
# 拟合ARIMA模型
fit <- auto.arima(q_ts)
# 得到2021年四个季度的预测值和95%置信区间
pred <- forecast(fit, h=4)
print(pred)

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