Optimized transmute_weights()

NEWS.md

@@ -8,6 +8,11 @@ averaged over the rolling window.
 - Fixed a bug where `transmute_weights()` and `factor_weights()` could return
 a result with a different length than `w`.
 
+- Added a new function `splice_index()` for splicing indexes calculated over
+a rolling window (this was previously sketched in an example).
+
+- `transmute_weights()` is now faster.
+
 ## Version 0.6.0
 
 - Bumped minimum version of R to at least 4.0.

R/means.R

@@ -423,12 +423,6 @@ harmonic_mean <- generalized_mean(-1)
 #' 1 / integrate(function(t) 1 / (2 * (1 - t) + 3 * t), 0, 1)$value
 #'
 #' @family means
-rdiff <- function(a, b, r, const = TRUE) {
-  if (r == 0) return(log(a / b))
-  res <- if (r == 1) a - b else a^r - b^r
-  if (const) res / r else res
-}
-
 #' @export
 extended_mean <- function(r, s) {
   r <- as.numeric(r)
@@ -439,18 +433,18 @@ extended_mean <- function(r, s) {
   if (not_finite_scalar(s)) {
     stop("'s' must be a finite length 1 numeric")
   }
-
+  
   function(a, b, tol = .Machine$double.eps^0.5) {
-    if (r == s) {
-      if (r == 0) {
-        res <- sqrt(a * b)
-      } else {
-        res <- exp(
-          (a^r * log(a) - b^r * log(b)) / rdiff(a, b, r, FALSE) - 1 / r
-        )
-      }
+    if (r == 0 && s == 0) {
+      res <- sqrt(a * b)
+    } else if (r == 0) {
+      res <- ((a^s - b^s) / log(a / b) / s)^(1 / s)
+    } else if (s == 0) {
+      res <- ((a^r - b^r) / log(a / b) / r)^(1 / r)
+    } else if (r == s) {
+      res <- exp((a^r * log(a) - b^r * log(b)) / (a^r - b^r) - 1 / r)
     } else {
-      res <- (rdiff(a, b, s) / rdiff(a, b, r))^(1 / (s - r))
+      res <- ((a^s - b^s) / (a^r - b^r) * r / s)^(1 / (s - r))
     }
     # set output to a when a == b
     i <- which(abs(a - b) <= tol)

R/splice.R

@@ -4,8 +4,8 @@
 #' index series. Splicing on multiple points combines the results with a
 #' geometric mean.
 #' 
-#' @param x A list of numeric vectors giving the period-over-period indexes for
-#' each window.
+#' @param x A list of equal-length numeric vectors giving the period-over-period
+#' indexes for each window.
 #' @param periods A vector (usually numeric) used to subscript each element of 
 #' `x` and give the splice points for each window. The default splices on each
 #' point in the window.
@@ -21,12 +21,15 @@
 #' x <- list(c(1.1, 0.9, 1.2), c(0.8, 1.3, 1.4), c(1.3, 1.3, 0.8))
 #' 
 #' # Mean splice
+#' 
 #' splice_index(x)
 #' 
 #' # Movement splice
+#' 
 #' splice_index(x, 3)
 #' 
 #' # Window splice
+#' 
 #' splice_index(x, 1)
 #' 
 #' @family price-indexes

R/weights.R

@@ -1,5 +1,15 @@
 #' Simplified extended mean for transmuting weights
 #' @noRd
+rdiff <- function(a, b, r) {
+  if (r == 0) {
+    log(a / b)
+  } else if (r == 1) {
+    a - b
+  } else {
+    a^r - b^r
+  }
+}
+
 extended_mean_ <- function(r, s) {
   r <- as.numeric(r)
   s <- as.numeric(s)
@@ -11,7 +21,7 @@ extended_mean_ <- function(r, s) {
   }
 
   function(x, m, tol = .Machine$double.eps^0.5) {
-    res <- rdiff(x, m, r, FALSE) / rdiff(x, m, s, FALSE)
+    res <- rdiff(x, m, r) / rdiff(x, m, s)
     res[abs(x - m) <= tol] <- m^(r - s)
     res
   }

_pkgdown.yml

@@ -20,6 +20,7 @@ reference:
   - price_indexes
   - index_weights
   - geks
+  - splice_index
 
 - title: Percent-change contributions
   contents:

tests/Examples/gpindex-Ex.Rout.save

@@ -1,6 +1,6 @@
 
-R version 4.3.2 (2023-10-31) -- "Eye Holes"
-Copyright (C) 2023 The R Foundation for Statistical Computing
+R version 4.3.3 (2024-02-29) -- "Angel Food Cake"
+Copyright (C) 2024 The R Foundation for Statistical Computing
 Platform: x86_64-pc-linux-gnu (64-bit)
 
 R is free software and comes with ABSOLUTELY NO WARRANTY.
@@ -483,6 +483,7 @@ Warning in back_period(period) :
 1.413257 1.835676 2.284565 2.789856 
 > 
 > # Calculate the index over a rolling window
+> 
 > (tg <- tornqvist_geks(price, quantity, period, product, window = 3))
 [[1]]
        2        3 
@@ -498,25 +499,13 @@ Warning in back_period(period) :
 
 > 
 > # Use a movement splice to combine the indexes in each window
-> cumprod(c(tg[[1]], sapply(tg[-1], \(x) x[length(x)])))
-       2        3        4        5 
-1.391443 1.801142 2.230521 2.689833 
+> 
+> splice_index(tg, 2)
+[1] 1.391443 1.801142 2.230521 2.689833
 > 
 > # ... or use a mean splice
-> mean_splice <- function(x, init) {
-+   offset <- length(init)
-+   x <- lapply(x, \(z) rev(cumprod(rev(z))))
-+   res <- numeric(offset + length(x))
-+   res[seq_along(init)] <- init
-+   for (i in seq_along(x)) {
-+     res[i + offset] <- geometric_mean(
-+       x[[i]] * res[seq(to = i + offset - 1, length.out = length(x[[i]]))]
-+     )
-+   }
-+   res
-+ }
 > 
-> mean_splice(tg[-1], cumprod(tg[[1]]))
+> splice_index(tg)
 [1] 1.391443 1.801142 2.228836 2.688842
 > 
 > #---- Missing data ----
@@ -1412,6 +1401,39 @@ t5   1.6720 1.0712    1.2477             0.9801 1.1850  1.2540  1.2678 0.9770
 > 
 > 
 > cleanEx()
+> nameEx("splice_index")
+> ### * splice_index
+> 
+> flush(stderr()); flush(stdout())
+> 
+> ### Name: splice_index
+> ### Title: Splice an index series
+> ### Aliases: splice_index
+> 
+> ### ** Examples
+> 
+> # Make an index series over a rolling window
+> x <- list(c(1.1, 0.9, 1.2), c(0.8, 1.3, 1.4), c(1.3, 1.3, 0.8))
+> 
+> # Mean splice
+> 
+> splice_index(x)
+[1] 1.100000 0.990000 1.188000 1.686819 1.306793
+> 
+> # Movement splice
+> 
+> splice_index(x, 3)
+[1] 1.10000 0.99000 1.18800 1.66320 1.33056
+> 
+> # Window splice
+> 
+> splice_index(x, 1)
+[1] 1.10000 0.99000 1.18800 1.60160 1.33848
+> 
+> 
+> 
+> 
+> cleanEx()
 > nameEx("transmute_weights")
 > ### * transmute_weights
 > 
@@ -1487,7 +1509,7 @@ t5   1.6720 1.0712    1.2477             0.9801 1.1850  1.2540  1.2678 0.9770
 +     x, nested_transmute2(0, c(1, -1), 2)(x, w1, w2)
 +   )
 + )
-[1] "Mean relative difference: 0.04366643"
+[1] "Mean relative difference: 17.81836"
 > 
 > #---- Monotonicity ----
 > 
@@ -1533,7 +1555,7 @@ t5   1.6720 1.0712    1.2477             0.9801 1.1850  1.2540  1.2678 0.9770
 > cleanEx()
 > options(digits = 7L)
 > base::cat("Time elapsed: ", proc.time() - base::get("ptime", pos = 'CheckExEnv'),"\n")
-Time elapsed:  0.35 0.044 0.445 0 0 
+Time elapsed:  0.095 0.023 0.118 0 0 
 > grDevices::dev.off()
 null device 
           1 

tests/testthat/test-means.R

@@ -57,6 +57,11 @@ test_that("generalized means satifies key properties", {
 test_that("generalized means work with transmuted weights", {
   expect_equal(generalized_mean(-4.4)(x),
                generalized_mean(0)(x, transmute_weights(-4.4, 0)(x)))
+
+  y <- c(x, geometric_mean(x))
+  expect_equal(generalized_mean(0)(y),
+               generalized_mean(2)(y, transmute_weights(0, 2)(y)))
+
   expect_equal(generalized_mean(3.8)(x, w),
                generalized_mean(-1)(x, transmute_weights(3.8, -1)(x, w)))
   expect_equal(

tests/testthat/test-splice.R

@@ -24,6 +24,24 @@ test_that("result length is correct", {
   )
 })
 
+test_that("splicing is invariant", {
+  # Movement
+  expect_equal(splice_index(x, 3), splice_index(x[-1], 3, x[[1]]))
+  expect_equal(splice_index(x, 3), splice_index(x[-(1:2)], 3, c(1, 2, 3, 5)))
+  # Window
+  expect_equal(splice_index(x, 1), splice_index(x[-1], 1, x[[1]]))
+  expect_equal(splice_index(x, 1), splice_index(x[-(1:2)], 1, c(1, 2, 3, 10)))
+  # Mean
+  expect_equal(splice_index(x), splice_index(x[-1], initial = x[[1]]))
+  expect_equal(
+    splice_index(x),
+    splice_index(
+      x[-(1:2)],
+      initial = c(1, 2, 3, geometric_mean(c(10, 40 / 6, 5)))
+    )
+  )
+})
+
 test_that("NAs return NA", {
   expect_equal(splice_index(x, c(1, NA, 4)), c(1, 2, 6, NA, NA))
   x[[2]][2] <- NA