Added splice_index()

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: gpindex
 Title: Generalized Price and Quantity Indexes
-Version: 0.6.0.9012
+Version: 0.6.0.9013
 Authors@R: c(
     person("Steve", "Martin", role = c("aut", "cre", "cph"),
     email = "[email protected]",
@@ -28,8 +28,8 @@ URL: https://marberts.github.io/gpindex/, https://github.com/marberts/gpindex
 BugReports: https://github.com/marberts/gpindex/issues
 LazyData: true
 Collate: 'helpers.R' 'means.R' 'weights.R' 'contributions.R' 'price_indexes.R'
-    'geks.R' 'operators.R' 'offset_prices.R' 'outliers.R' 'price_data.R'
-    'gpindex-package.R'
+    'geks.R' 'splice.R' 'operators.R' 'offset_prices.R' 'outliers.R'
+    'price_data.R' 'gpindex-package.R'
 Config/testthat/edition: 3
 VignetteBuilder: knitr
 Roxygen: list(markdown = TRUE)

NAMESPACE

@@ -52,6 +52,7 @@ export(quartile_method)
 export(resistant_fences)
 export(robust_z)
 export(scale_weights)
+export(splice_index)
 export(stuvel_index)
 export(tornqvist_geks)
 export(transmute_weights)

R/contributions.R

@@ -248,7 +248,7 @@ contributions <- function(r) {
   arithmetic_weights <- transmute_weights(r, 1)
 
   function(x, w = NULL) {
-    (x - 1) * scale_weights(arithmetic_weights(x, w))
+    (x - 1) * arithmetic_weights(x, w)
   }
 }
 
@@ -276,7 +276,7 @@ nc <- function(nest_transmute) {
     arithmetic_weights <- nest_transmute(r1, r2, 1, t)
 
     function(x, w1 = NULL, w2 = NULL) {
-      (x - 1) * scale_weights(arithmetic_weights(x, w1, w2))
+      (x - 1) * arithmetic_weights(x, w1, w2)
     }
   }
 }

R/geks.R

@@ -120,26 +120,16 @@ geks_matrix <- function(index, p, q, product, n, nper, window, na.rm) {
 #' cumprod(tornqvist_geks(price, quantity, period, product)[[1]])
 #'
 #' # Calculate the index over a rolling window
+#' 
 #' (tg <- tornqvist_geks(price, quantity, period, product, window = 3))
 #' 
 #' # Use a movement splice to combine the indexes in each window
-#' cumprod(c(tg[[1]], sapply(tg[-1], \(x) x[length(x)])))
+#' 
+#' splice_index(tg, 2)
 #' 
 #' # ... 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)
 #'
 #' #---- Missing data ----
 #'

R/means.R

@@ -214,6 +214,8 @@ generalized_mean <- function(r) {
       }
       if (r == 0) {
         exp(sum(log(x)) / length(x))
+      } else if (r == 1) {
+        sum(x) / length(x)
       } else {
         (sum(x^r) / length(x))^(1 / r)
       }
@@ -228,6 +230,8 @@ generalized_mean <- function(r) {
       }
       if (r == 0) {
         exp(sum(log(x) * w) / sum(w))
+      } else if (r == 1) {
+        sum(x * w) / sum(w)
       } else {
         (sum(x^r * w) / sum(w))^(1 / r)
       }
@@ -419,6 +423,12 @@ 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)
@@ -431,16 +441,16 @@ extended_mean <- function(r, s) {
   }
 
   function(a, b, tol = .Machine$double.eps^0.5) {
-    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)
+    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
+        )
+      }
     } else {
-      res <- ((a^s - b^s) / (a^r - b^r) * r / s)^(1 / (s - r))
+      res <- (rdiff(a, b, s) / rdiff(a, b, r))^(1 / (s - r))
     }
     # set output to a when a == b
     i <- which(abs(a - b) <= tol)

R/splice.R

@@ -0,0 +1,67 @@
+#' Splice an index series
+#' 
+#' Splice a collection of index series computed over a rolling window into one
+#' 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 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.
+#' @param initial A numeric vector giving an initial period-over-period index
+#' series onto which the elements of `x` are spliced. The default uses the
+#' first element of `x`.
+#' 
+#' @returns
+#' A numeric vector giving the spliced (fixed-base) index series.
+#' 
+#' @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)
+#' 
+#' # Movement splice
+#' splice_index(x, 3)
+#' 
+#' # Window splice
+#' splice_index(x, 1)
+#' 
+#' @family price-indexes
+#' @export
+splice_index <- function(x, periods = NULL, initial = NULL) {
+  x <- as.list(x)
+  if (do.call(different_lengths, x)) {
+    stop("all elements of 'x' must be the same length")
+  }
+
+  if (is.null(initial) && length(x) > 0L) {
+    initial <- x[[1L]]
+    x <- x[-1L]
+  } 
+  initial <- cumprod(initial)
+  if (length(x) == 0L) {
+    return(initial)
+  }
+
+  offset <- length(initial)
+  if (offset < length(x[[1L]])) {
+    stop("'initial' must be at least as long as each element of 'x'")
+  }
+
+  if (is.null(periods)) {
+    periods <- seq_along(x[[1L]])
+  }
+
+  x <- lapply(x, \(z) rev(cumprod(rev(z))))
+  res <- numeric(offset + length(x))
+  res[seq_along(initial)] <- initial
+  for (i in seq_along(x)) {
+    window <- x[[i]]
+    iw <- seq.int(to = i + offset - 1L, length.out = length(window))
+    res[i + offset] <- geometric_mean(window[periods] * res[iw[periods]])
+  }
+  res
+}

R/weights.R

@@ -1,3 +1,22 @@
+#' Simplified extended mean for transmuting weights
+#' @noRd
+extended_mean_ <- function(r, s) {
+  r <- as.numeric(r)
+  s <- as.numeric(s)
+  if (not_finite_scalar(r)) {
+    stop("'r' must be a finite length 1 numeric")
+  }
+  if (not_finite_scalar(s)) {
+    stop("'s' must be a finite length 1 numeric")
+  }
+
+  function(x, m, tol = .Machine$double.eps^0.5) {
+    res <- rdiff(x, m, r, FALSE) / rdiff(x, m, s, FALSE)
+    res[abs(x - m) <= tol] <- m^(r - s)
+    res
+  }
+}
+
 #' Transmute weights
 #'
 #' Transmute weights to turn a generalized mean of order \eqn{r} into a
@@ -166,28 +185,28 @@ transmute_weights <- function(r, s) {
   r <- as.numeric(r)
   s <- as.numeric(s)
   gen_mean <- generalized_mean(r)
-  ext_mean <- extended_mean(r, s)
+  ext_mean <- extended_mean_(r, s)
 
   function(x, w = NULL) {
     if (r == s) {
       if (is.null(w)) {
-        w <- rep.int(1, length(x))
+        return(rep.int(1 / length(x), length(x)))
       }
       if (length(x) != length(w)) {
         stop("'x' and 'w' must be the same length")
       }
       w[is.na(x)] <- NA_real_
-      w
+      scale_weights(w)
     } else {
       m <- gen_mean(x, w, na.rm = TRUE)
       if (is.null(w)) {
-        v <- ext_mean(x, m)^(r - s)
+        v <- ext_mean(x, m)
         attributes(v) <- NULL
       } else {
-        v <- w * ext_mean(x, m)^(r - s)
+        v <- w * ext_mean(x, m)
         attributes(v) <- attributes(w)
       }
-      v
+      scale_weights(v)
     }
   }
 }
@@ -213,12 +232,12 @@ nested_transmute <- function(r1, r2, s, t = c(1, 1)) {
 
   function(x, w1 = NULL, w2 = NULL) {
     if (is.na(t[1L]) && !is.na(t[2L])) {
-      w <- scale_weights(r_weights2(x, w2))
+      w <- r_weights2(x, w2)
     } else if (!is.na(t[1L]) && is.na(t[2L])) {
-      w <- scale_weights(r_weights1(x, w1))
+      w <- r_weights1(x, w1)
     } else {
-      v1 <- scale_weights(r_weights1(x, w1))
-      v2 <- scale_weights(r_weights2(x, w2))
+      v1 <- r_weights1(x, w1)
+      v2 <- r_weights2(x, w2)
       # the calculation is wrong if NAs in w1 or w2 propagate
       if (anyNA(w1)) {
         v1[is.na(v1) & !is.na(v2)] <- 0
@@ -257,12 +276,12 @@ nested_transmute2 <- function(r1, r2, s, t = c(1, 1)) {
     m <- c(mean1(x, w1, na.rm = TRUE), mean2(x, w2, na.rm = TRUE))
     v <- s_weights(m, t)
     if (is.na(v[1L]) && !is.na(v[2L])) {
-      scale_weights(s_weights2(x, w2))
+      s_weights2(x, w2)
     } else if (!is.na(v[1L]) && is.na(v[2L])) {
-      scale_weights(s_weights1(x, w1))
+      s_weights1(x, w1)
     } else {
-      u1 <- scale_weights(s_weights1(x, w1))
-      u2 <- scale_weights(s_weights2(x, w2))
+      u1 <- s_weights1(x, w1)
+      u2 <- s_weights2(x, w2)
       # the calculation is wrong if NAs in w1 or w2 propagate
       if (anyNA(w1)) {
         u1[is.na(u1) & !is.na(u2)] <- 0