将函数应用于时间序列对象

Question

假设我在 R 中有以下 ts 对象。

x <- ts(data = matrix(1:10, 5, 2), start = 1/12, deltat = 1/12)

现在我想对这个时间序列的每个条目应用一个函数f：函数f取决于时间序列的值和对应的时间，例如：

f <- function(z, time){z*time}

我正在寻找一种有效的方法来做到这一点：到目前为止，我只管理了一些解决方案，如下所示：

timex <- seq(from = 1/12, by = 1/12, length = 5)
apply(x, 2, function(y){ apply(cbind(y, timex),1, function(z) f(z[1], z[2]))})

此解决方案提供了正确的结果，但我很确定还有更直接的方法。

我正在寻找一种适用于更“复杂功能”f 的方法：特别是如果 length(x)>1，则无法直接调用 f(x, time(x))。

对我来说，时间序列是否为 ts 格式并不重要，所以如果有针对不同时间序列格式的解决方案，那我完全赞成。对我来说，这里最重要的是性能。 最终结果也不必是时间序列对象。

你们有什么提示吗？ 谢谢。

编辑： 问题似乎是 x 一旦使用 apply 就会失去其 ts 属性。 特别是不能再使用 time 命令了。

apply(x, 2, class)
 Series 1  Series 2 
"integer" "integer"

class(x)
[1] "mts"    "ts"     "matrix"

编辑 2： 我被要求提供一个“更复杂的函数 f： 这是一个计算某些金融产品价值的函数。

f  <- function(s0, t){
  option.times <- 1:20/2
  tau <- option.times[option.times >= t] - t

  d = (log(1.5/s0) - 0.0098*tau)/(0.0098* sqrt(tau))

  p1 <- pnorm(d)
  p2 <- s0 / 1.5* exp((1/2*0.02^2 + 0.0098)*tau) * pnorm(d-0.02*sqrt(tau)) 
  sum((p1-p2)*exp(-0.02*tau))*100

}

我无法提供我正在使用的时间序列：但它是一个包含每周数据的 530x10000 时间序列（时间（x）从 0 开始，并以 deltat = 1/52 进行）。 对于这个函数 f 和我的时间序列，我明白了

> system.time(apply(x, 2, function(y){ apply(cbind(y, timex),1, function(z) f(z[1], z[2]))}))
       User      System verstrichen 
      69.62        0.03       69.75 
> system.time(mapply(f, x, time(x)))
       User      System verstrichen 
      79.01        0.06       79.12 

所以 mapply 比使用 apply 两次稍慢。

Answer 1

您可以使用mapply 来使用多个向量作为输入：

x <- ts(data = matrix(1:10, 5, 2), start = 1/12, deltat = 1/12)
f <- function(z, time){z*time}
mapply(f, x, time(x))
#>  [1] 0.08333333 0.33333333 0.75000000 1.33333333 2.08333333 0.50000000
#>  [7] 1.16666667 2.00000000 3.00000000 4.16666667

编辑： 我可以用真实的功能和真实的数据量重现你的时间：

n <- 530
m <- 10000
x <- ts(data = matrix(rlnorm(n * m), n, m), start = 0, deltat = 1/52)
f  <- function(s0, t){
  option.times <- 1:20/2
  tau <- option.times[option.times >= t] - t

  d = (log(1.5/s0) - 0.0098*tau)/(0.0098* sqrt(tau))

  p1 <- pnorm(d)
  p2 <- s0 / 1.5* exp((1/2*0.02^2 + 0.0098)*tau) * pnorm(d-0.02*sqrt(tau)) 
  sum((p1-p2)*exp(-0.02*tau))*100
}
timex <- time(x)
system.time(r1 <- apply(x, 2, function(y){ apply(cbind(y, timex),1, function(z) f(z[1], z[2]))}))
#>        User      System verstrichen 
#>      67.002       0.059      67.089
system.time(r2 <- matrix(mapply(f, x, time(x)), n, m))
#>        User      System verstrichen 
#>      78.975       0.244      79.250
all(r1 == r2)
#> [1] TRUE

但是，您的函数至少允许部分矢量化，因此足以对所有行进行（显式）循环，而不是对所有矩阵元素进行（隐式）循环：

g  <- function(s0, t){
  option.times <- 1:20/2
  tau <- option.times[option.times >= t] - t

  d = outer(-0.0098*tau, log(1.5/s0), FUN = "+")/(0.0098 * sqrt(tau))

  p1 <- pnorm(d)
  p2 <- outer(exp((1/2*0.02^2 + 0.0098)*tau), s0 / 1.5, FUN = "*") * pnorm(d - 0.02*sqrt(tau)) 
  colSums((p1 - p2)*exp(-0.02*tau))*100
}

r3 <- matrix(0, n ,m)
timex <- time(x)
system.time(for (i in seq_along(timex)) {
  r3[i, ] <- g(x[i, ], timex[i])
})
#>        User      System verstrichen 
#>       4.955       0.136       4.919
all(r3 == r2)
#> [1] TRUE

可能有一些方法可以完全矢量化函数......