如何使用不均匀和缺失数据创建按子集分组的滚动预测？

问题描述

我有一个数据集（160 万行，感兴趣的 4 列），由定向国家/地区 dyad-year 组织。每个非交换二元组（year1-stateA-stateB 并不总是等于 year1-stateB-stateA）都有一个输出值“var1”。

简化数据示例

library(forecast)
library(dplyr)

df=data.frame(year=c(1994,1995,1996,1997,1998,1964,1965,1967,1968,1969,1988,1987,1989),stateA=c(1,1,138,20,20),stateB=c(2,2,87,55,55),var1=c(0.101,0.132,0.136,0.148,-0.287,-0.112,0.088,0.101,0.121,0.387,NA,0.377,0.388)
)

> df
   year stateA stateB   var1
1  1994      1      2  0.101
2  1995      1      2  0.132
3  1996      1      2  0.136
4  1997      1      2  0.136
5  1998      1      2  0.148
6  1964    138     87 -0.287
7  1965    138     87 -0.112
8  1967    138     87  0.088
9  1968    138     87  0.101
10 1969    138     87  0.121
11 1988     20     55  0.387
12 1987     20     55     NA
13 1988     20     55  0.377
14 1989     20     55  0.388

我想要做的是将每组国家/地区二元组分解为时间序列，并使用 holt 模型使用过去 5 年的数据创建下一年的预测预测。

预期结果：我希望添加一个新变量，其中包含基于前几年的 yearX+1 的预测值到 yearX 的行。

并发症：并非每个国家/地区对每年都存在，并且在某些年份，尽管数据集中存在国家/地区对联，但没有数据。

到目前为止我所做的：

首先，请原谅我最近才开始在 R 中使用时间序列。

首先，我使用 dplr 按年份组织数据（因此它将按正确的时间序列顺序排列）然后按 stateA、stateB 分组

 rolldata <- df %>%
  dplyr::arrange(year) %>% 
  dplyr::group_by(stateA,stateB) %>% [...]

我之前所做的是 5 年滚动平均值，它不符合我的分析需求，所以看起来像这样：

rolldata <- df %>%
  dplyr::arrange(year) %>% 
  dplyr::group_by(stateA,stateB) %>% 
  dplyr::mutate(
    point_5a = zoo::rollmean(var1,k = 5,fill = NA,align='right'))

这里的问题是我需要为每一行创建一个时间序列对象以传递给 holt() 以输出预测值 (fvar)。

dat_ts <- ts(df$var1,start = c(STARTYEAR,1),end = c(ROWYEAR,frequency = 1)
holt_model <- holt(dat_ts,h = 5)
fvar[i] <-holt_model$x[1]

我希望我已经以易于理解的方式阐述了这个问题。非常感谢您的帮助，我随时准备澄清和回答任何可能对您有帮助的问题。

附言效率不是必须的，只有结果。

编辑：我认为我之前没有说清楚，但我的主要目标是为每一行而不是整个子集生成一个预测对象。在我的国家 1 和国家 2 的示例数据中：将根据 1994 年的时间序列对 1994 年进行预测；将在 1994-1995 年的基础上对 1995 年作出预测；基于 1994-1996 年的 1996 年预测。然后对 (138,87) 也是如此，每一行都有自己的预测。

解决方法

paste 将列一起用于使用 by 进行拆分/应用。 holt 警告因缺失而采取的行动。你可以插入它们，我不确定，但你想用它们做什么。

rolldata <- df[order(df$year),]

library(forecast)
res <- by(rolldata,Reduce(paste,rolldata[c("stateA","stateB")]),function(x) {
  STARTYEAR <- x$year[1]; ROWYEAR <- x$year[nrow(x)]
  x0 <- x$var1
  x_ts <- ts(x0,start = c(STARTYEAR,1),end = c(ROWYEAR,frequency = 1)
  holt_mod <- holt(x_ts,h=5)
  holt_mod$x[1]
})
# Warning message:
#   In ets(x,"AAN",alpha = alpha,beta = beta,phi = phi,damped = damped,:
#            Missing values encountered. Using longest contiguous portion of time series

as.list(res)
# $`1 2`
# [1] 0.101
# 
# $`138 87`
# [1] -0.287
# 
# $`20 55`
# [1] 0.387

,

您可以在数据帧本身中存储类 forecast 的对象并从中提取相关值。

library(dplyr)
library(purrr)
library(forecast)

df %>%
  arrange(year) %>% 
  group_by(stateA,stateB) %>%
  summarise(ts = list(ts(var1,start = c(min(year),end = c(max(year),frequency = 1)),holt = map(ts,holt,h = 5),all_value = map(holt,~.x$x),first_value = map_dbl(all_value,first)) %>%
    ungroup

#  stateA stateB ts       holt       all_value first_value
#   <dbl>  <dbl> <list>   <list>     <list>          <dbl>
#1      1      2 <ts [5]> <forecast> <ts [5]>        0.101
#2     20     55 <ts [3]> <forecast> <ts [2]>        0.387
#3    138     87 <ts [6]> <forecast> <ts [6]>       -0.287