Skip to content

Commit

Permalink
Added splicing and walsh index example
Browse files Browse the repository at this point in the history
  • Loading branch information
marberts committed Jan 26, 2024
1 parent c9eb65e commit 0a002ce
Show file tree
Hide file tree
Showing 6 changed files with 137 additions and 78 deletions.
2 changes: 1 addition & 1 deletion DESCRIPTION
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
Package: gpindex
Title: Generalized Price and Quantity Indexes
Version: 0.6.0.9006
Version: 0.6.0.9008
Authors@R: c(
person("Steve", "Martin", role = c("aut", "cre", "cph"),
email = "[email protected]",
Expand Down
37 changes: 29 additions & 8 deletions R/geks.R
Original file line number Diff line number Diff line change
Expand Up @@ -109,16 +109,36 @@ geks_matrix <- function(index, p, q, product, n, nper, window, na.rm) {
#' *Journal of Econometrics*, 161(1): 24--35.
#'
#' @examples
#' price <- 1:6
#' quantity <- 6:1
#' period <- rep(1:3, 2)
#' product <- rep(letters[1:2], each = 3)
#' price <- 1:10
#' quantity <- 10:1
#' period <- rep(1:5, 2)
#' product <- rep(letters[1:2], each = 5)
#'
#' tornqvist_geks(price, quantity, period, product)
#' cumprod(tornqvist_geks(price, quantity, period, product)[[1]])
#'
#' tornqvist_geks(price, quantity, period, product, window = 2)
#' # 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)])))
#'
#' # ... 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_len(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]]))
#'
#' # Missing data
#' #---- Missing data ----
#'
#' quantity[2] <- NA
#'
Expand All @@ -131,7 +151,8 @@ geks_matrix <- function(index, p, q, product, n, nper, window, na.rm) {
#' fg <- geks(balanced(fisher_index))
#' fg(price, quantity, period, product, na.rm = TRUE)
#'
#' # Make a Jevons GEKS index
#' #---- Make a Jevons GEKS index ----
#'
#' jevons_geks <- geks(\(p1, p0, ..., na.rm) jevons_index(p1, p0, na.rm))
#' jevons_geks(price, quantity, period, product)
#'
Expand Down
20 changes: 8 additions & 12 deletions R/means.R
Original file line number Diff line number Diff line change
Expand Up @@ -676,18 +676,9 @@ contraharmonic_mean <- lehmer_mean(2)
#' # A function to make the superlative quadratic mean price index by
#' # Diewert (1976) as a product of generalized means
#'
#' quadratic_mean_index <- function(x, w0, w1, r) {
#' x <- sqrt(x)
#' generalized_mean(r)(x, w0) * generalized_mean(-r)(x, w1)
#' }
#'
#' quadratic_mean_index(x, w1, w2, 2)
#'
#' # Same as the nested generalized mean (with the order halved)
#'
#' quadratic_mean_index2 <- function(r) nested_mean(0, c(r / 2, -r / 2))
#' quadratic_mean_index <- function(r) nested_mean(0, c(r / 2, -r / 2))
#'
#' quadratic_mean_index2(2)(x, w1, w2)
#' quadratic_mean_index(2)(x, w1, w2)
#'
#' # The arithmetic AG mean index by Lent and Dorfman (2009)
#'
Expand All @@ -705,7 +696,7 @@ contraharmonic_mean <- lehmer_mean(2)
#' q1 <- quantity6[[2]]
#' q0 <- quantity6[[1]]
#'
#' walsh <- quadratic_mean_index2(1)
#' walsh <- quadratic_mean_index(1)
#'
#' sum(p1 * q1) / sum(p0 * q0) / walsh(q1 / q0, p0 * q0, p1 * q1)
#'
Expand All @@ -715,6 +706,11 @@ contraharmonic_mean <- lehmer_mean(2)
#' # quadratic mean price index of order 1
#'
#' walsh(p1 / p0, p0 * q0, p1 * q1)
#'
#' # That requires using the average value share as weights
#'
#' walsh_weights <- sqrt(scale_weights(p0 * q0) * scale_weights(p1 * q1))
#' walsh(p1 / p0, walsh_weights, walsh_weights)
#'
#' #---- Missing values ----
#'
Expand Down
39 changes: 30 additions & 9 deletions man/geks.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

20 changes: 8 additions & 12 deletions man/nested_mean.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

97 changes: 61 additions & 36 deletions tests/Examples/gpindex-Ex.Rout.save
Original file line number Diff line number Diff line change
Expand Up @@ -473,54 +473,83 @@ Warning in back_period(period) :
>
> ### ** Examples
>
> price <- 1:6
> quantity <- 6:1
> period <- rep(1:3, 2)
> product <- rep(letters[1:2], each = 3)
> price <- 1:10
> quantity <- 10:1
> period <- rep(1:5, 2)
> product <- rep(letters[1:2], each = 5)
>
> tornqvist_geks(price, quantity, period, product)
[[1]]
2 3
1.530869 1.376227

> cumprod(tornqvist_geks(price, quantity, period, product)[[1]])
2 3 4 5
1.413257 1.835676 2.284565 2.789856
>
> tornqvist_geks(price, quantity, period, product, window = 2)
> # Calculate the index over a rolling window
> (tg <- tornqvist_geks(price, quantity, period, product, window = 3))
[[1]]
2
1.520408
2 3
1.391443 1.294442

[[2]]
3
1.366822
3 4
1.292486 1.238393

[[3]]
4 5
1.238417 1.205921

>
> # Missing data
> # 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
>
> # ... 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_len(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]]))
Warning in seq_len(init) : first element used of 'length.out' argument
Warning in res[seq_len(init)] <- init :
number of items to replace is not a multiple of replacement length
[1] 1.391443 0.000000 0.000000 0.000000
>
> #---- Missing data ----
>
> quantity[2] <- NA
>
> # Use all non-missing data
>
> fisher_geks(price, quantity, period, product, na.rm = TRUE)
[[1]]
2 3
1.455853 1.370478
2 3 4 5
1.438137 1.234230 1.234212 1.216746

>
> # Remove records with any missing data
>
> fg <- geks(balanced(fisher_index))
> fg(price, quantity, period, product, na.rm = TRUE)
[[1]]
2 3
1.403078 1.346954
2 3 4 5
1.501481 1.148250 1.219688 1.199513

>
> # Make a Jevons GEKS index
> #---- Make a Jevons GEKS index ----
>
> jevons_geks <- geks(\(p1, p0, ..., na.rm) jevons_index(p1, p0, na.rm))
> jevons_geks(price, quantity, period, product)
[[1]]
2 3
1.581139 1.341641
2 3 4 5
1.527525 1.309307 1.224745 1.178511

>
>
Expand Down Expand Up @@ -962,19 +991,9 @@ Warning in back_period(period) :
> # A function to make the superlative quadratic mean price index by
> # Diewert (1976) as a product of generalized means
>
> quadratic_mean_index <- function(x, w0, w1, r) {
+ x <- sqrt(x)
+ generalized_mean(r)(x, w0) * generalized_mean(-r)(x, w1)
+ }
>
> quadratic_mean_index(x, w1, w2, 2)
[1] 1.912366
>
> # Same as the nested generalized mean (with the order halved)
>
> quadratic_mean_index2 <- function(r) nested_mean(0, c(r / 2, -r / 2))
> quadratic_mean_index <- function(r) nested_mean(0, c(r / 2, -r / 2))
>
> quadratic_mean_index2(2)(x, w1, w2)
> quadratic_mean_index(2)(x, w1, w2)
[1] 1.912366
>
> # The arithmetic AG mean index by Lent and Dorfman (2009)
Expand All @@ -994,7 +1013,7 @@ Warning in back_period(period) :
> q1 <- quantity6[[2]]
> q0 <- quantity6[[1]]
>
> walsh <- quadratic_mean_index2(1)
> walsh <- quadratic_mean_index(1)
>
> sum(p1 * q1) / sum(p0 * q0) / walsh(q1 / q0, p0 * q0, p1 * q1)
[1] 1.401718
Expand All @@ -1008,6 +1027,12 @@ Warning in back_period(period) :
> walsh(p1 / p0, p0 * q0, p1 * q1)
[1] 1.401534
>
> # That requires using the average value share as weights
>
> walsh_weights <- sqrt(scale_weights(p0 * q0) * scale_weights(p1 * q1))
> walsh(p1 / p0, walsh_weights, walsh_weights)
[1] 1.401718
>
> #---- Missing values ----
>
> x[1] <- NA
Expand Down Expand Up @@ -1511,7 +1536,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.335 0.025 0.361 0 0
Time elapsed: 0.35 0.044 0.445 0 0
> grDevices::dev.off()
null device
1
Expand Down

0 comments on commit 0a002ce

Please sign in to comment.