Making price indexes

Most price indexes are made with a two-step procedure, where period-over-period elemental indexes are calculated for a collection of elemental aggregates at each point in time, and then aggregated according to a price index aggregation structure. These indexes can then be chained together to form a time series that gives the evolution of prices with respect to a fixed base period. The piar package contains a collections of functions that revolve around this work flow, making it easy to build standard price indexes in R.

The purpose of this vignette is to give several extended examples of how to use the functions in this package to make different types of price indexes. This should serve both as a introduction to the functionality in piar, and a reference for solving specific index-number problems.

Matched-sample index

The first example covers calculating a matched-sample index, where a fixed set of businesses each provide prices for a collection of products over time. The products reported by a businesses can change over time, but the set of businesses is fixed for the duration of the sample. Each businesses has a weight that is established when the sample is drawn, and represents a particular segment of the economy.

The usual approach for calculating a matched-sample index starts by computing the elemental index for each business as an equally-weighted geometric mean of price relatives (i.e., a Jevons index). From there, index values for different segments of the economy are calculated as an arithmetic mean of the elemental indexes, using the businesses-level weights (either a Young or Lowe index, depending how the weights are constructed).

The ms_prices dataset has price data for five businesses over four quarters, and the ms_weights dataset has the weight data. Note that these data have fairly realistic patterns of missing data.

library(piar)

options(stringsAsFactors = FALSE)

head(ms_prices)
#>   period business product price
#> 1 202001       B1       1  1.14
#> 2 202001       B1       2    NA
#> 3 202001       B1       3  6.09
#> 4 202001       B2       4  6.23
#> 5 202001       B2       5  8.61
#> 6 202001       B2       6  6.40

ms_weights
#>   business classification weight
#> 1       B1             11    553
#> 2       B2             11    646
#> 3       B3             11    312
#> 4       B4             12    622
#> 5       B5             12    330

The elemental_index() function makes, well, elemental indexes, using information on price relatives, elemental aggregates (businesses), and time periods (quarters). By default it makes a Jevons index, but any bilateral generalized-mean index is possible. The only wrinkle is that price data here are in levels, and not relatives, but the price_relative() function can make the necessary conversion.

relative <- with(ms_prices, price_relative(price, period, product))

(ms_epr <- with(ms_prices, elemental_index(relative, period, business, na.rm = TRUE)))
#>    202001    202002    202003   202004
#> B1      1 0.8949097 0.3342939      NaN
#> B2      1       NaN       NaN 2.770456
#> B3      1 2.0200036 1.6353355 0.537996
#> B4    NaN       NaN       NaN 4.576286

(Homogeneous elemental aggregates often leads to unit-value elemental indexes that are not based on price relatives. These cases can be dealt with by first aggregating prices for each elemental aggregate, aggregate(price ~ period + product, ms_prices, mean), at each point in time with an arithmetic mean, then forming price relatives to feed into elemental_index().)

As with most functions in R, missing values are contagious by default in piar. Setting na.rm = TRUE in elemental_index() means that missing price relatives are ignored, which is equivalent to imputing these missing relatives with the value of the elemental index for the respective businesses (i.e., parental or overall mean imputation). Other types of imputation are possible, and are the topic of a subsequent example.

The elemental_index() function returns a special index object, and there are a number of methods for working with these objects. Probably the most useful of these methods allows the resulting elemental indexes to be extracted like a matrix, even though it’s not a matrix. (Note that there are only indexes for four businesses, not five, because the fifth business never reports any prices; an elemental index can be made for this business with a small change to the call to elemental_index().)

ms_epr[, "202004"]
#>      202004
#> B1      NaN
#> B2 2.770456
#> B3 0.537996
#> B4 4.576286
ms_epr["B1", ]
#>    202001    202002    202003 202004
#> B1      1 0.8949097 0.3342939    NaN

With the elemental indexes out of the way, it’s time to make a price-index aggregation structure that maps each business to its position in the aggregation hierarchy. The only hiccup is unpacking the digit-wise classification for each businesses that defines the hierarchy. That’s the job of the expand_classification() function.

hierarchy <- with(ms_weights, c(expand_classification(classification), list(business)))

pias <- aggregation_structure(hierarchy, ms_weights$weight)

It is now simple to aggregate the elemental indexes according to this aggregation structure with the aggregate() function. As with the elemental indexes, missing values are ignored by setting na.rm = TRUE, which is equivalent to parentally imputing missing values. Note that, unlike the elemental indexes, missing values are filled in to ensure the index can be chained over time.

(ms_index <- aggregate(ms_epr, pias, na.rm = TRUE))
#>    202001    202002    202003   202004
#> 1       1 1.3007239 1.0630743 2.734761
#> 11      1 1.3007239 1.0630743 1.574515
#> 12      1 1.3007239 1.0630743 4.576286
#> B1      1 0.8949097 0.3342939 1.574515
#> B2      1 1.3007239 1.0630743 2.770456
#> B3      1 2.0200036 1.6353355 0.537996
#> B4      1 1.3007239 1.0630743 4.576286
#> B5      1 1.3007239 1.0630743 4.576286

Although simple, this example covers the core functionality of piar. The remaining examples in the vignette build on this one by adding complexities that often arise in practice.

Chaining

The elemental_index() function makes period-over-period elemental indexes by default, which can then be aggregated to make a period-over-period index. Chaining an index is the process of taking the cumulative product of each of these period-over-period indexes to make a time series that compares prices to a fixed base period.

The chain() function can be used to chain the values in an index object.

(ms_index_chained <- chain(ms_index))
#>    202001    202002    202003    202004
#> 1       1 1.3007239 1.3827662 3.7815355
#> 11      1 1.3007239 1.3827662 2.1771866
#> 12      1 1.3007239 1.3827662 6.3279338
#> B1      1 0.8949097 0.2991629 0.4710366
#> B2      1 1.3007239 1.3827662 3.8308934
#> B3      1 2.0200036 3.3033836 1.7772072
#> B4      1 1.3007239 1.3827662 6.3279338
#> B5      1 1.3007239 1.3827662 6.3279338

This gives almost the same result as directly manipulating the index as a matrix, except that the former returns an index object (not a matrix).

t(apply(as.matrix(ms_index), 1, cumprod))
#>    202001    202002    202003    202004
#> 1       1 1.3007239 1.3827662 3.7815355
#> 11      1 1.3007239 1.3827662 2.1771866
#> 12      1 1.3007239 1.3827662 6.3279338
#> B1      1 0.8949097 0.2991629 0.4710366
#> B2      1 1.3007239 1.3827662 3.8308934
#> B3      1 2.0200036 3.3033836 1.7772072
#> B4      1 1.3007239 1.3827662 6.3279338
#> B5      1 1.3007239 1.3827662 6.3279338

Chained indexes often need be to rebased, and this can be done with the rebase() function. For example, rebasing the index so that 202004 is the base period just requires dividing the chained index by the slice for 202004.

rebase(ms_index_chained, ms_index_chained[, "202004"])
#>       202001    202002    202003 202004
#> 1  0.2644428 0.3439671 0.3656626      1
#> 11 0.4593084 0.5974334 0.6351161      1
#> 12 0.1580295 0.2055527 0.2185178      1
#> B1 2.1229774 1.8998731 0.6351161      1
#> B2 0.2610357 0.3395354 0.3609514      1
#> B3 0.5626806 1.1366169 1.8587499      1
#> B4 0.1580295 0.2055527 0.2185178      1
#> B5 0.1580295 0.2055527 0.2185178      1

In some cases the base period is the average of several periods; setting the base period to the second half of 2020 just requires dividing by the row-wise mean.

rebase(ms_index_chained, rowMeans(as.matrix(ms_index_chained)[, c("202003", "202004")]))
#>       202001    202002    202003    202004
#> 1  0.3872740 0.5037366 0.5355095 1.4644905
#> 11 0.5618052 0.7307535 0.7768452 1.2231548
#> 12 0.2593798 0.3373815 0.3586616 1.6413384
#> B1 2.5967299 2.3238388 0.7768452 1.2231548
#> B2 0.3836077 0.4989677 0.5304398 1.4695602
#> B3 0.3936550 0.7951845 1.3003935 0.6996065
#> B4 0.2593798 0.3373815 0.3586616 1.6413384
#> B5 0.2593798 0.3373815 0.3586616 1.6413384

Multi-dimensional aggregation structures

Price indexes are often aggregated over multiple dimensions. Matched sample indexes that use sequential Poisson sampling are a good example, as there are usually take-all and take-some strata in addition to, say, an industry classification.

(ms_weights <- transform(ms_weights, stratum = c("TS", "TA", "TS", "TS", "TS")))
#>   business classification weight stratum
#> 1       B1             11    553      TS
#> 2       B2             11    646      TA
#> 3       B3             11    312      TS
#> 4       B4             12    622      TS
#> 5       B5             12    330      TS

The easiest way to deal with multiple digit-wise classifications is to turn them into one classification. In this example the “stratum” dimension comes before the “classification” dimension for the purposes of parental imputation.

(classification_sps <- with(ms_weights, paste0(classification, stratum)))
#> [1] "11TS" "11TA" "11TS" "12TS" "12TS"

This classification can be expanded with the expand_classification() function as before, just with an extra instruction to say that the last “digit” in the classification is two characters wide, not one.

(classification_sps <- expand_classification(classification_sps, width = c(1, 1, 2)))
#> [[1]]
#> [1] "1" "1" "1" "1" "1"
#> 
#> [[2]]
#> [1] "11" "11" "11" "12" "12"
#> 
#> [[3]]
#> [1] "11TS" "11TA" "11TS" "12TS" "12TS"
pias_sps <- with(
  ms_weights, 
  aggregation_structure(c(classification_sps, list(business)), weight)
)

The elemental indexes can now be aggregated according to this new aggregation structure.

aggregate(ms_epr, pias_sps, na.rm = TRUE)
#>      202001    202002    202003   202004
#> 1         1 1.3007239 1.0630743 2.684412
#> 11        1 1.3007239 1.0630743 1.492443
#> 12        1 1.3007239 1.0630743 4.576286
#> 11TA      1 1.3007239 1.0630743 2.770456
#> 11TS      1 1.3007239 1.0630743 0.537996
#> 12TS      1 1.3007239 1.0630743 4.576286
#> B1        1 0.8949097 0.3342939 0.537996
#> B2        1 1.3007239 1.0630743 2.770456
#> B3        1 2.0200036 1.6353355 0.537996
#> B4        1 1.3007239 1.0630743 4.576286
#> B5        1 1.3007239 1.0630743 4.576286

Non-parental imputation during aggregation

Parental imputation is the usual way to impute missing index values during aggregation, and it is simple to do with aggregate(). In some cases, however, a business-level index may get imputed with the value for, say, another business, rather than for an entire group of businesses. The simplest way to do this sort of imputation is to alter the elemental indexes prior to aggregation. It is also possible to augment the aggregation structure with an imputation layer, but this is more complex.

Suppose that missing index values for business B2 should be imputed as 1, rather than the value for group 11. This replacement can be done as if the index was a matrix.

ms_epr2 <- ms_epr
ms_epr2["B2", 2:3] <- 1
ms_epr2
#>    202001    202002    202003   202004
#> B1      1 0.8949097 0.3342939      NaN
#> B2      1 1.0000000 1.0000000 2.770456
#> B3      1 2.0200036 1.6353355 0.537996
#> B4    NaN       NaN       NaN 4.576286

The index can now be aggregated as usual.

aggregate(ms_epr2, pias, na.rm = TRUE)
#>    202001    202002    202003   202004
#> 1       1 1.1721550 1.0400686 2.626560
#> 11      1 1.1721550 1.0400686 1.398142
#> 12      1 1.1721550 1.0400686 4.576286
#> B1      1 0.8949097 0.3342939 1.398142
#> B2      1 1.0000000 1.0000000 2.770456
#> B3      1 2.0200036 1.6353355 0.537996
#> B4      1 1.1721550 1.0400686 4.576286
#> B5      1 1.1721550 1.0400686 4.576286

Alternate index-number formulas

By default, the elemental_index() function calculates a Jevons index. Although this is the standard index-number formula for making elemental indexes, many other types of index-numbers are possible. The Carli index (equally-weighted arithmetic mean of price relatives) is the main competitor to the Jevons, and requires specifying the order of the index r when calling elemental_index(). An order of 1 corresponds to an arithmetic mean.

with(ms_prices, elemental_index(relative, period, business, na.rm = TRUE, r = 1))
#>    202001     202002    202003   202004
#> B1      1  0.8949097 0.3342939      NaN
#> B2      1        NaN       NaN 5.155942
#> B3      1 23.7480455 2.4900997 0.607197
#> B4    NaN        NaN       NaN 9.368610

The Coggeshall index (equally-weighted harmonic mean of price relatives) is another competitor to the Jevons, but is seldom used in practice. Despite it being more exotic, it is just as easy to make by specifying an order r of -1.

with(ms_prices, elemental_index(relative, period, business, na.rm = TRUE, r = -1))
#>    202001    202002    202003    202004
#> B1      1 0.8949097 0.3342939       NaN
#> B2      1       NaN       NaN 1.7205750
#> B3      1 0.6591433 0.8185743 0.4746769
#> B4    NaN       NaN       NaN 2.2353790

The type of mean used to aggregate elemental indexes can be controlled in the same way in the call to aggregate(). The default makes an arithmetic index, but any type of generalized-mean index is possible.

Many superlative indexes can be made by supplying unequal and (usually) time-varying weights when making the elemental indexes. These weights often come from information on quantities.

ms_prices2 <- transform(ms_prices, quantity = 10 - price)

The Tornqvist index is a popular superlative index-number formula, using average period-over-period value shares as the weights in a geometric mean. The only tricky part is making the weights from data on prices and quantities.

library(gpindex)

tw <- grouped(index_weights("Tornqvist"))

ms_prices2[c("back_price", "back_quantity")] <- 
  ms_prices2[back_period(ms_prices2$period, ms_prices2$product), c("price", "quantity")]

ms_prices2 <- na.omit(ms_prices2) # can't have NAs for Tornqvist weights

ms_prices2$weight <- with(
  ms_prices2,
  tw(price, back_price, quantity, back_quantity, group = interaction(period, business))
)

As elemental_index() makes a geometric index by default, all that is needed to make a Tornqvist index is to provide the weights.

with(ms_prices2, elemental_index(price / back_price, period, business, weight))
#>    202001    202002    202003   202004
#> B1      1 0.8949097 0.3342939      NaN
#> B2      1       NaN       NaN 2.165152
#> B3      1 0.9520982 1.5913929 0.542372
#> B4    NaN       NaN       NaN 5.904237

Quote contributions

It’s often convenient to decompose an index into the (additive) contribution of each price relative, also known as the quote contributions. This can be done with the same work flow as in the previous examples, specifying contrib = TRUE when calling elemental_index().

ms_epr <- with(ms_prices, elemental_index(relative, period, business, contrib = TRUE, na.rm = TRUE))

As with index values, quote contributions for a given level of the index can be extracted as a matrix.

contrib(ms_epr)
#>   202001     202002     202003 202004
#> 1      0  0.0000000  0.0000000      0
#> 2     NA         NA -0.6657061      0
#> 3      0 -0.1050903         NA     NA

Aggregating the elemental indexes automatically aggregates quote contributions, so no extra steps are needed after the elemental indexes are made.

Updating a basket

The functions in piar are all designed to work within a “basket”, which is a fancy way of saying within a given aggregation structure. Over time, however, aggregation structures change as the weights used to aggregate an index get updated, and new samples of businesses are drawn. The general approach to keep the time series going is to “chain” the index across baskets.

It is easier to see how to chain an index over time with a simple example that just splits the ms_prices data in two.

ms_prices1 <- subset(ms_prices, period <= "202003")
ms_prices2 <- subset(ms_prices, period >= "202003")

The index for the first basket can be calculate as usual.

ms_epr1 <- with(
  ms_prices1, 
  elemental_index(price_relative(price, period, product), period, business, na.rm = TRUE)
)

(ms_index1 <- aggregate(ms_epr1, pias, na.rm = TRUE))
#>    202001    202002    202003
#> 1       1 1.3007239 1.0630743
#> 11      1 1.3007239 1.0630743
#> 12      1 1.3007239 1.0630743
#> B1      1 0.8949097 0.3342939
#> B2      1 1.3007239 1.0630743
#> B3      1 2.0200036 1.6353355
#> B4      1 1.3007239 1.0630743
#> B5      1 1.3007239 1.0630743

Nothing special needs to be done to make the elemental indexes for the new basket, but it’s easier to remove the index values of 1 for quarter 3 2020.

ms_epr2 <- with(
  subset(transform(ms_prices2, rel = price_relative(price, period, product)), period > "202003"),
  elemental_index(rel, period, business, na.rm = TRUE)
)

Aggregating these elemental indexes, however, requires an aggregation structure. The results of the first example can be reproduced by simply “price updating” the original weights, then building the aggregation structure as usual.

(ms_index2 <- aggregate(ms_epr2, update(pias, ms_index1), na.rm = TRUE))
#>      202004
#> 1  2.734761
#> 11 1.574515
#> 12 4.576286
#> B1 1.574515
#> B2 2.770456
#> B3 0.537996
#> B4 4.576286
#> B5 4.576286

This produces two sets of period-over-period indexes that can be stacked together and then chained.

chain(stack(ms_index1, ms_index2))
#>    202001    202002    202003    202004
#> 1       1 1.3007239 1.3827662 3.7815355
#> 11      1 1.3007239 1.3827662 2.1771866
#> 12      1 1.3007239 1.3827662 6.3279338
#> B1      1 0.8949097 0.2991629 0.4710366
#> B2      1 1.3007239 1.3827662 3.8308934
#> B3      1 2.0200036 3.3033836 1.7772072
#> B4      1 1.3007239 1.3827662 6.3279338
#> B5      1 1.3007239 1.3827662 6.3279338

Carry-forward imputation

The previous examples used parental imputation to both impute missing price relatives when calculating the elemental indexes, and to impute missing elemental indexes during aggregation. Another common imputation strategy when making elemental indexes is to carry forward the previous price to impute for missing prices, and parentally impute missing elemental indexes during aggregation. As the elemental_index() function accepts price relatives as its input, other types of imputations can be done prior to passing price relatives to this function.

(ms_epr2 <- with(
  ms_prices, 
  elemental_index(price_relative(carry_forward(price, period, product), period, product), 
                  period, business, na.rm = TRUE)
  )
)
#>    202001    202002    202003   202004
#> B1      1 0.8949097 0.5781816 1.000000
#> B2      1 1.0000000 0.1777227 2.770456
#> B3      1 2.0200036 1.6353355 0.537996
#> B4    NaN       NaN       NaN 4.576286

Aggregation follows the same steps as in the previous examples, with missing values set to be ignored in order to parentally impute missing elemental indexes.

(ms_index <- aggregate(ms_epr2, pias, na.rm = TRUE))
#>    202001    202002    202003    202004
#> 1       1 1.1721550 0.8082981 2.2653614
#> 11      1 1.1721550 0.8082981 0.8093718
#> 12      1 1.1721550 0.8082981 4.5762862
#> B1      1 0.8949097 0.5781816 1.0000000
#> B2      1 1.0000000 0.1777227 2.7704563
#> B3      1 2.0200036 1.6353355 0.5379960
#> B4      1 1.1721550 0.8082981 4.5762862
#> B5      1 1.1721550 0.8082981 4.5762862

Calculating an index from multiple sources

All of the examples so far have built an index from a single source of price data. In many cases the elemental indexes are built from multiple sources of data, either because no single source of data has the necessary coverage, or different index-number formulas are employed for different elemental aggregates.

It is straightforward to merge index objects together, provided they’re for the same time periods. To keep the example simple, suppose that ms_prices is split in two.

ms_prices1 <- subset(ms_prices, business %in% c("B1", "B2", "B3"))
ms_prices2 <- subset(ms_prices, business == "B4")

Elemental indexes can be made for both groups separately with the usual recipe. Note that there is no data for business B4 in the first two periods, so the time periods need to be made explicit.

ms_epr1 <- with(
  ms_prices1, 
  elemental_index(price_relative(price, period, product), period, business, na.rm = TRUE)
)
ms_epr1
#>    202001    202002    202003   202004
#> B1      1 0.8949097 0.3342939      NaN
#> B2      1       NaN       NaN 2.770456
#> B3      1 2.0200036 1.6353355 0.537996
ms_epr2 <- with(
  transform(ms_prices2, period = factor(period, levels = time(ms_epr1))), 
  elemental_index(price_relative(price, period, product), period, business, na.rm = TRUE)
)
ms_epr2
#>    202001 202002 202003   202004
#> B4    NaN    NaN    NaN 4.576286

Once the elemental indexes are made, they can be merged together and then aggregated.

aggregate(merge(ms_epr1, ms_epr2), pias, na.rm = TRUE)
#>    202001    202002    202003   202004
#> 1       1 1.3007239 1.0630743 2.734761
#> 11      1 1.3007239 1.0630743 1.574515
#> 12      1 1.3007239 1.0630743 4.576286
#> B1      1 0.8949097 0.3342939 1.574515
#> B2      1 1.3007239 1.0630743 2.770456
#> B3      1 2.0200036 1.6353355 0.537996
#> B4      1 1.3007239 1.0630743 4.576286
#> B5      1 1.3007239 1.0630743 4.576286

A slightly more complex case is when some of the input data are already a price index. For example, suppose the index values for businesses B4 and B5 come from some outside process, and are taken as inputs.

ms_prices2 <- subset(
  as.data.frame(aggregate(ms_epr, pias, na.rm = TRUE)),
  level %in% c("B4", "B5")
)
ms_prices2
#>    period level    value
#> 7  202001    B4 1.000000
#> 8  202001    B5 1.000000
#> 15 202002    B4 1.300724
#> 16 202002    B5 1.300724
#> 23 202003    B4 1.063074
#> 24 202003    B5 1.063074
#> 31 202004    B4 4.576286
#> 32 202004    B5 4.576286

All that is required is to pass the pre-existing indexes to elemental_index() to cast it into the correct form. This won’t affect their values, but will allow them to be merged with the other elemental indexes, and aggregated.

ms_epr2 <- with(ms_prices2, elemental_index(value, period, level))
aggregate(merge(ms_epr1, ms_epr2), pias, na.rm = TRUE)
#>    202001    202002    202003   202004
#> 1       1 1.3007239 1.0630743 2.734761
#> 11      1 1.3007239 1.0630743 1.574515
#> 12      1 1.3007239 1.0630743 4.576286
#> B1      1 0.8949097 0.3342939 1.574515
#> B2      1 1.3007239 1.0630743 2.770456
#> B3      1 2.0200036 1.6353355 0.537996
#> B4      1 1.3007239 1.0630743 4.576286
#> B5      1 1.3007239 1.0630743 4.576286

Paasche index

All of the examples so far have used a single set of weights to aggregate an index. Although this is by far the most common case, there are situations where the weights to aggregate change every period. The Paasche index is the notable example, as the weights for aggregation are the current-period revenue shares in each period.

weights <- data.frame(period = rep(c("202001", "202002", "202003", "202004"), each = 5),
                      classification = ms_weights$classification,
                      weight = 1:20)
head(weights)
#>   period classification weight
#> 1 202001             11      1
#> 2 202001             11      2
#> 3 202001             11      3
#> 4 202001             12      4
#> 5 202001             12      5
#> 6 202002             11      6

The only new tools needed to deal with time-varying weights are the stack() and unstack() functions. stack() appends a later index series onto an earlier one for the same levels, whereas unstack() pulls apart an index series for many periods into a collection of one-period indexes.

The first step to making a Paasche index is to unstack the elemental indexes into a list of elemental indexes for each period. (Trying to make the elemental indexes period-by-period can be dangerous when there are missing values.)

(ms_epr <- unstack(ms_epr))
#> [[1]]
#>    202001
#> B1      1
#> B2      1
#> B3      1
#> B4    NaN
#> 
#> [[2]]
#>       202002
#> B1 0.8949097
#> B2       NaN
#> B3 2.0200036
#> B4       NaN
#> 
#> [[3]]
#>       202003
#> B1 0.3342939
#> B2       NaN
#> B3 1.6353355
#> B4       NaN
#> 
#> [[4]]
#>      202004
#> B1      NaN
#> B2 2.770456
#> B3 0.537996
#> B4 4.576286

The second step is to make a sequence of aggregation structures for each set of weights.

pias <- with(weights, 
             Map(aggregation_structure, 
                 list(hierarchy), 
                 split(weight, period))
)

Making the Paasche index for each period is now just a case of mapping the aggregate() function to each elemental index and aggregation structure, and then reducing the result with the stack() function.

(paasche <- Reduce(stack, Map(aggregate, ms_epr, pias, na.rm = TRUE, r = -1)))
#>    202001    202002    202003    202004
#> 1       1 1.3127080 0.5874490 1.3591916
#> 11      1 1.3127080 0.5874490 0.8839797
#> 12      1 1.3127080 0.5874490 4.5762862
#> B1      1 0.8949097 0.3342939 0.8839797
#> B2      1 1.3127080 0.5874490 2.7704563
#> B3      1 2.0200036 1.6353355 0.5379960
#> B4      1 1.3127080 0.5874490 4.5762862
#> B5      1 1.3127080 0.5874490 4.5762862

With a Paasche index in hand, it is now trivial to make a Fisher index by first making the chained Laspeyres index, and then doing a simple matrix operation.

laspeyres <- Reduce(stack, Map(aggregate, ms_epr, pias[c(1, 1, 2, 3)], na.rm = TRUE))
sqrt(as.matrix(laspeyres) * as.matrix(paasche))
#>    202001    202002    202003   202004
#> 1       1 1.5107763 0.7956890 1.996688
#> 11      1 1.5107763 0.7956890 1.192826
#> 12      1 1.5107763 0.7956890 4.576286
#> B1      1 0.8949097 0.3342939 1.192826
#> B2      1 1.5107763 0.7956890 2.770456
#> B3      1 2.0200036 1.6353355 0.537996
#> B4      1 1.5107763 0.7956890 4.576286
#> B5      1 1.5107763 0.7956890 4.576286