将不等时间间隔的xts重采样为等距时间序列

Question

我正在使用xts时间序列进行R语言编程。

我的问题: 一个具有不等间隔时间步的时间序列数据集。

我想要得到的结果: 一个具有等间隔时间步的时间序列，其值对应于重叠时间步的原始值比例（见下面的示例）。

例如：使用原始时间序列如下：

sample_xts <- as.xts(read.zoo(text='
2016-07-01 00:00:20,   0.0
2016-07-01 00:01:20,  60.0
2016-07-01 00:01:50,  30.0
2016-07-01 00:02:30,  40.0
2016-07-01 00:04:20, 110.0
2016-07-01 00:05:30, 140.0
2016-07-01 00:06:00,  97.0
2016-07-01 00:07:12, 144.0
2016-07-01 00:08:09,   0.0
', sep=',', index=1, tz='', format="%Y-%m-%d %H:%M:%S"))
names(sample_xts) <- c('x')

我希望得到一个均匀间隔的时间序列，就像这样:

                         x
2016-07-01 00:00:00,   0.0
2016-07-01 00:01:00,  40.0
2016-07-01 00:02:00,  60.0
2016-07-01 00:03:00,  60.0
2016-07-01 00:04:00,  60.0
2016-07-01 00:05:00, 100.0
2016-07-01 00:06:00, 157.0
2016-07-01 00:07:00, 120.0
2016-07-01 00:08:00,  24.0
2016-07-01 00:09:00,   0.0

注意：

• 一些原始时间步长小于新时间步长，而其他时间步长大于新时间步长。
• x的colSums保持不变（即621）。

这是我用来创建上面示例的草图（可能有助于说明我想做什么）： 我希望的方法不仅限于创建1分钟时间步长系列，而是适用于任何固定时间步长。

我已经查看了许多stackoverflow上的q / a并尝试了许多不同的方法，但没有成功。

非常感谢任何帮助！谢谢。

- Morten Grum

Answer 1

这里是我用zoo编写的一些代码-我没有使用过xts，所以不知道是否可以应用相同的函数。希望有所帮助！

函数

以下函数计算原始数据每个时间间隔与给定时间间隔重叠的分数（注意：在以下所有代码中，变量名ta1和ta2指代给定时间间隔的开始和结束（例如需要作为输出的每个等间隔），而tb1和tb2则指代原始数据的（不等）时间间隔的开始和结束）：

frac.overlap <- function(ta1,ta2,tb1,tb2){
if(tb1 <= ta1 & tb2 >= ta2) {   # Interval 2 starts earlier and ends later than interval 1
    frac <- as.numeric(difftime(ta2,ta1,units="secs"))/as.numeric(difftime(tb2,tb1,units="secs"))
} else if(tb1 >= ta1 & tb2 <= ta2) {    # Interval 2 is fully contained within interval 1
    frac <- 1
} else if(tb1 <= ta1 & tb2 >= ta1) {    # Interval 2 partly overlaps with interval 1 (starts earlier, ends earlier)
    frac <- as.numeric(difftime(tb2,ta1,units="secs"))/as.numeric(difftime(tb2,tb1,units="secs"))
} else if (tb1 <= ta2 & tb2 >= ta2){    # Interval 2 partly overlaps with interval 1 (starts later, ends later)
    frac <- as.numeric(difftime(ta2,tb1,units="secs"))/as.numeric(difftime(tb2,tb1,units="secs"))
        } else {                                # No overlap
            frac <- 0
    }

    return(frac)
}

下一个函数确定原始数据集中与当前考虑的区间 ta1-ta2 重叠的记录：

check.overlap <- function(ta1,ta2,tb1,tb2){
ov <- vector("logical",4)
ov[1] <- (tb1 <= ta1 & tb2 >= ta2)  # Interval 2 starts earlier and ends later than interval 1
ov[2] <- (tb1 >= ta1 & tb2 <= ta2)  # Interval 2 is fully contained within interval 1
ov[3] <- (tb1 <= ta1 & tb2 >= ta1)  # Interval 2 partly overlaps with interval 1 (starts earlier, ends earlier)
ov[4] <- (tb1 <= ta2 & tb2 >= ta2)  # Interval 2 partly overlaps with interval 1 (starts later, ends later)
return(as.logical(sum(ov))) # Gives TRUE if at least one element of ov is TRUE, otherwise FALSE
}

（注意：这在您提供的示例数据中效果很好，但在更大的数据集上，我发现它速度非常慢。由于我编写此代码以重新采样具有固定时间步长的时间序列，因此我通常使用固定间隔来完成此步骤，这样可以显著提高速度。很可能很容易修改代码（请参见下一个函数的代码），以根据原始数据的间隔加快此步骤。）

下一个函数使用前两个函数来计算间隔ta1-ta2的重新采样值：

fracres <- function(tstart,interval,input){
# tstart: POSIX object
# interval: length of interval in seconds
# input: zoo object

ta1 <- tstart
ta2 <- tstart + interval

# First, determine which records of the original data (input) overlap with the current
# interval, to avoid going through the whole object at every iteration
ind <- index(input)
ind1 <- index(lag(input,-1))
recs <- which(sapply(1:length(ind),function(x) check.overlap(ta1,ta2,ind[x],ind1[x])))
#recs <- which(abs(as.numeric(difftime(ind,ta1,units="secs"))) < 601)


# For each record overlapping with the current interval, return the fraction of the input data interval contained in the current interval
if(length(recs) > 0){
    fracs <- sapply(1:length(recs), function(x) frac.overlap(ta1,ta2,ind[recs[x]],ind1[recs[x]]))
    return(sum(coredata(input)[recs]*fracs))

} else {
    return(0)
}
}

（如果原始时间步长和新时间步长之间的最大时间差已知，则可以使用注释掉的行来获取相关记录。）

应用程序

首先，将您的示例数据作为zoo对象读入：

sample_zoo <- read.zoo(text='
2016-07-01 00:00:20,   0.0
2016-07-01 00:01:20,  60.0
2016-07-01 00:01:50,  30.0
2016-07-01 00:02:30,  40.0
2016-07-01 00:04:20, 110.0
2016-07-01 00:05:30, 140.0
2016-07-01 00:06:00,  97.0
2016-07-01 00:07:12, 144.0
2016-07-01 00:08:09,   0.0
', sep=',', index=1, tz='', format="%Y-%m-%d %H:%M:%S")

看起来你的数据集包含瞬时值（“在 01:20 时，x 的值为60”）。由于我编写的代码是针对累加值的，时间戳的含义不同（“从 01:20 开始的记录的值为60”）。为了纠正这个问题，需要对记录进行移位：

sample_zoo <- lag(sample_zoo,1)

然后，我们定义一个序列的POSIXct对象，对应所需的分辨率：

time.out <- seq.POSIXt(from=as.POSIXct("2016-07-01"),to=(as.POSIXct("2016-07-01")+(60*9)),by="1 min")

我们可以接着使用上述描述的fracres函数：

data.out <- sapply(1:length(time.out), function(x) fracres(tstart=time.out[x],interval=60,input=sample_zoo))

索引和数据被合并到一个名为 zoo 的对象中：

zoo.out <- read.zoo(data.frame(time.out,data.out))

最后，时间序列再次向相反的方向移动一步：

zoo.out <- lag(zoo.out,-1)

2016-07-01 00:01:00 2016-07-01 00:02:00 2016-07-01 00:03:00 2016-07-01 00:04:00 2016-07-01 00:05:00 2016-07-01 00:06:00 2016-07-01 00:07:00 2016-07-01 00:08:00 2016-07-01 00:09:00 
             40                  60                  60                  60                 100                 157                 120                  24                   0 

- m.chips

---

感谢@m.chips！终于在我的实时系列中尝试了这个。它完美地工作，但是正如你指出的那样，即使是相当短的系列也会变得“难以承受的缓慢”。似乎执行时间与系列长度不成比例——呈指数或2^N增长。我的系列有30万到100万个观测值，所以受你算法的启发，我决定尝试其他方法。将其发布在问题的答案中。 - Morten Grum

Answer 2

我最终决定使用“while循环”的方式，并创建了以下解决方案。它的效果很好 - 并不是超级快，但执行时间似乎与时间序列的长度成比例。我将其与我在问题中发布的小示例以及一个源时间序列有330,000个观测值和一个约110,000个时间步骤的目标系列进行了测试。

源系列和目标系列都可以拥有不规则的时间步长。生成系列的总和与源相同。

性能: 它速度还算可以，但我相信它可以更快。我猜这是一个明显的Rcpp版本候选项，对于长系列来说应该会快得多。现在对我来说这样做就可以了，如果/当我开始创建Rcpp版本时，我会在这里发布。

如果您有性能改进建议，请发表评论！

谢谢！

sameEndTime <- function(i,j,src_index,dest_index){
  if(src_index[i] == dest_index[j]){
    TRUE
  } else {
    FALSE
  }
}

wholeSourceStepIsWithinDestStep <- function(i,j,src_index,dest_index){
  if(dest_index[j-1] <= src_index[i-1] & src_index[i] <= dest_index[j]){
    TRUE
  } else {
    FALSE
  }
}

wholeDestStepIsWithinSourceStep <- function(i,j,src_index,dest_index){
  if(src_index[i-1] <= dest_index[j-1]  &  dest_index[j] <= src_index[i]){
    TRUE
  } else {
    FALSE
  }
}

onlyEndOfSourceStepIsWithinDestStep <- function(i,j,src_index,dest_index){
  if(src_index[i-1] < dest_index[j-1]  &  src_index[i] < dest_index[j] & src_index[i] > dest_index[j-1]){
    TRUE
  } else {
    FALSE
  }
}

onlyStartOfSourceStepIsWithinDestStep <- function(i,j,src_index,dest_index){
  if(src_index[i-1] < dest_index[j]  &  src_index[i-1] > dest_index[j-1] & src_index[i] > dest_index[j]){
    TRUE
  } else {
    FALSE
  }
}

resampleToDestTimeSteps <- function(src, dest){
  # src and dest are both xts with only one time series each
  # src is the original series and 
  # dest holds the time steps of the final series
  #
  # NB: there is an issue with the very first time step 
  # (gets ignored in this version)
  #
  original_names <- names(src)
  names(src) <- c("value")
  names(dest) <- c("value")
  dest$value <- dest$value*0.0
  dest$value[is.na(dest$value)] <- 0.0

  dest[1]$value = 0.0

  for(k in 2:length(src)){
    src[k]$value <- src[k]$value/as.numeric(difftime(index(src[k]),index(src[k-1]),units="secs"))
  }
  # First value is NA due to lag at this point (we don't want that)
  src$value[1] = 0.0

  i = 2 # source timestep counter
  j = 2 # destination timestep counter

  src_index = index(src)
  dest_index = index(dest)

  src_length = length(src)
  dest_length = length(dest)

  # Make sure we start with an overlap
  if(src_index[2] < dest_index[1]){
    while(src_index[i] < dest_index[1]){
      i = i + 1
    }
  } else if(dest_index[2] < src_index[1]){
    while(dest_index[j] < src_index[1]){
      j = j + 1
    }
  }

  while(i <= src_length & j <= dest_length){
    if( wholeSourceStepIsWithinDestStep(i,j,src_index,dest_index) ){
      dest[j]$value = dest[j]$value + as.numeric(src[i]$value)*as.numeric(difftime(src_index[i],src_index[i-1],units="secs"))
      if(sameEndTime(i,j,src_index,dest_index)){
        j = j+1
      }
      i = i+1
    } else if( wholeDestStepIsWithinSourceStep(i,j,src_index,dest_index) ){
      dest[j]$value = dest[j]$value + as.numeric(src[i]$value)*as.numeric(difftime(dest_index[j],dest_index[j-1],units="secs"))
      if(sameEndTime(i,j,src_index,dest_index)){
        i = i+1
      }
      j = j+1
    } else if( onlyEndOfSourceStepIsWithinDestStep(i,j,src_index,dest_index) ){
      dest[j]$value = dest[j]$value + as.numeric(src[i]$value)*as.numeric(difftime(src_index[i],dest_index[j-1],units="secs"))
      i = i+1
    } else if( onlyStartOfSourceStepIsWithinDestStep(i,j,src_index,dest_index) ){
      diff_time = difftime(dest_index[j],src_index[i-1],units="secs")
      dest[j]$value = dest[j]$value + as.numeric(src[i]$value)*as.numeric(diff_time)
      j = j+1
    } else {
      print("======================================================")
      print(paste0("i=",i,", j=",j))
      print(paste0("src_index[i]   =",src_index[i]))
      print(paste0("dest_index[j]  =",dest_index[j]))
      print(" ")
      print(paste0("src_index[i-1] =",src_index[i-1]))
      print(paste0("dest_index[j-1]=",dest_index[j-1]))
      print("======================================================")
      stop("This should never happen.")
    }
  }
  names(dest) <- original_names
  return(dest)
}

- Morten Grum