diff --git a/DESCRIPTION b/DESCRIPTION index 78cf5a5..4bdcf11 100644 --- a/DESCRIPTION +++ b/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 = "marberts@protonmail.com", @@ -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) diff --git a/NAMESPACE b/NAMESPACE index 532ebfb..27c1618 100644 --- a/NAMESPACE +++ b/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) diff --git a/R/contributions.R b/R/contributions.R index 84584fd..a4bcea1 100644 --- a/R/contributions.R +++ b/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) } } } diff --git a/R/geks.R b/R/geks.R index a48fe46..641a4ab 100644 --- a/R/geks.R +++ b/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 ---- #' diff --git a/R/means.R b/R/means.R index 8a69276..b8bb72d 100644 --- a/R/means.R +++ b/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) diff --git a/R/splice.R b/R/splice.R new file mode 100644 index 0000000..4d037ab --- /dev/null +++ b/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 +} \ No newline at end of file diff --git a/R/weights.R b/R/weights.R index c67c97a..3c065e6 100644 --- a/R/weights.R +++ b/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 diff --git a/man/geks.Rd b/man/geks.Rd index 14c607a..6339584 100644 --- a/man/geks.Rd +++ b/man/geks.Rd @@ -108,26 +108,16 @@ product <- rep(letters[1:2], each = 5) 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 ---- @@ -166,6 +156,7 @@ GEKS index with more options. Other price-indexes: \code{\link{index_weights}()}, -\code{\link{price_indexes}} +\code{\link{price_indexes}}, +\code{\link{splice_index}()} } \concept{price-indexes} diff --git a/man/index_weights.Rd b/man/index_weights.Rd index f48a513..a604040 100644 --- a/man/index_weights.Rd +++ b/man/index_weights.Rd @@ -180,6 +180,7 @@ quantity index. Other price-indexes: \code{\link{geks}()}, -\code{\link{price_indexes}} +\code{\link{price_indexes}}, +\code{\link{splice_index}()} } \concept{price-indexes} diff --git a/man/price_indexes.Rd b/man/price_indexes.Rd index 6e79d91..27df0be 100644 --- a/man/price_indexes.Rd +++ b/man/price_indexes.Rd @@ -372,6 +372,7 @@ multiple groups of products over many time periods. Other price-indexes: \code{\link{geks}()}, -\code{\link{index_weights}()} +\code{\link{index_weights}()}, +\code{\link{splice_index}()} } \concept{price-indexes} diff --git a/man/splice_index.Rd b/man/splice_index.Rd new file mode 100644 index 0000000..fd17699 --- /dev/null +++ b/man/splice_index.Rd @@ -0,0 +1,49 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/splice.R +\name{splice_index} +\alias{splice_index} +\title{Splice an index series} +\usage{ +splice_index(x, periods = NULL, initial = NULL) +} +\arguments{ +\item{x}{A list of numeric vectors giving the period-over-period indexes for +each window.} + +\item{periods}{A vector (usually numeric) used to subscript each element of +\code{x} and give the splice points for each window. The default splices on each +point in the window.} + +\item{initial}{A numeric vector giving an initial period-over-period index +series onto which the elements of \code{x} are spliced. The default uses the +first element of \code{x}.} +} +\value{ +A numeric vector giving the spliced (fixed-base) index series. +} +\description{ +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. +} +\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) + +} +\seealso{ +Other price-indexes: +\code{\link{geks}()}, +\code{\link{index_weights}()}, +\code{\link{price_indexes}} +} +\concept{price-indexes} diff --git a/tests/testthat/test-splice.R b/tests/testthat/test-splice.R new file mode 100644 index 0000000..116fd1a --- /dev/null +++ b/tests/testthat/test-splice.R @@ -0,0 +1,44 @@ +x <- list(1:3, 3:5, 5:7) + +test_that("key splice methods works", { + # Movement splice. + expect_equal(splice_index(x, 3), c(1, 2, 6, 30, 210)) + # Window splice. + expect_equal(splice_index(x, 1), c(1, 2, 6, 60, 420)) + # Mean splice. + expect_equal( + splice_index(x), + c(1, 2, 6, geometric_mean(c(60, 40, 30)), + geometric_mean(c(420, 252, 7 * geometric_mean(c(60, 40, 30))))) + ) +}) + +test_that("result length is correct", { + expect_equal( + splice_index(x[-1], 3, c(1, 1, 1, 1, 2, 3)), + c(1, 1, 1, 1, 2, 6, 30, 210) + ) + expect_equal( + splice_index(x[-1], 1, c(1, 1, 1, 1, 2, 3)), + c(1, 1, 1, 1, 2, 6, 60, 420) + ) +}) + +test_that("NAs return NA", { + expect_equal(splice_index(x, c(1, NA, 4)), c(1, 2, 6, NA, NA)) + x[[2]][2] <- NA + expect_equal(splice_index(x), c(1, 2, 6, NA, NA)) + expect_equal(splice_index(x, 3), c(1, 2, 6, 30, 210)) +}) + +test_that("corner cases work", { + expect_equal(splice_index(list()), numeric(0)) + expect_equal(splice_index(list(), initial = 1:5), cumprod(1:5)) + expect_equal(splice_index(list(1:5)), cumprod(1:5)) + expect_equal(splice_index(list(1, 2, 3, 4, 5)), cumprod(1:5)) +}) + +test_that("errors work", { + expect_error(splice_index(list(1:3, 1:2))) + expect_error(splice_index(list(1:3, 1:3), initial = 1:2)) +}) \ No newline at end of file diff --git a/tests/testthat/test-weights.R b/tests/testthat/test-weights.R index 67266c0..2620d27 100644 --- a/tests/testthat/test-weights.R +++ b/tests/testthat/test-weights.R @@ -5,13 +5,16 @@ w <- runif(15, 0, 2) f <- factor(sample(letters[1:3], 15, TRUE)) test_that("weights transmute correctly", { - expect_equal(transmute_weights(2, 2)(x), rep(1, length(x))) - expect_equal(transmute_weights(0, 0)(xna, w), replace(w, 2, NA)) + expect_equal(transmute_weights(2, 2)(x), rep(1 / 15, length(x))) + expect_equal( + transmute_weights(0, 0)(xna, w), + scale_weights(replace(w, 2, NA)) + ) expect_equal(transmute_weights(2, 1)(c(1, NA)), c(1, NA)) - expect_equal(scale_weights(transmute_weights(-1, 1)(x, w)), - scale_weights(w / x)) + expect_equal(transmute_weights(-1, 1)(x, w), scale_weights(w / x)) expect_equal( - transmute_weights(7, -3)(x, transmute_weights(-3, 7)(x, w)), w + transmute_weights(7, -3)(x, transmute_weights(-3, 7)(x, w)), + scale_weights(w) ) expect_equal( grouped(transmute_weights(1, 2))(x, w, group = f),