diff --git a/DESCRIPTION b/DESCRIPTION
index 7097a7d..5ebf82b 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: piar
 Title: Price Index Aggregation
-Version: 0.8.3.9005
+Version: 0.8.3.9006
 Authors@R: c(
     person("Steve", "Martin", role = c("aut", "cre", "cph"),
            email = "marberts@protonmail.com",
diff --git a/NAMESPACE b/NAMESPACE
index 22a0154..d60d47b 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -33,10 +33,6 @@ S3method(as_index,default)
 S3method(as_index,direct_piar_index)
 S3method(as_index,matrix)
 S3method(as_index,mts)
-S3method(carry_backward,data.frame)
-S3method(carry_backward,default)
-S3method(carry_forward,data.frame)
-S3method(carry_forward,default)
 S3method(chain,chainable_piar_index)
 S3method(chain,default)
 S3method(chain,direct_piar_index)
@@ -46,6 +42,10 @@ S3method(elementary_index,default)
 S3method(elementary_index,numeric)
 S3method(end,piar_index)
 S3method(head,piar_index)
+S3method(impute_prices,data.frame)
+S3method(impute_prices,default)
+S3method(impute_prices,matrix)
+S3method(impute_prices,numeric)
 S3method(is.na,piar_index)
 S3method(levels,piar_aggregation_structure)
 S3method(levels,piar_index)
@@ -61,8 +61,6 @@ S3method(print,piar_index)
 S3method(rebase,chainable_piar_index)
 S3method(rebase,default)
 S3method(rebase,direct_piar_index)
-S3method(shadow_price,data.frame)
-S3method(shadow_price,default)
 S3method(split,piar_index)
 S3method(stack,chainable_piar_index)
 S3method(stack,direct_piar_index)
@@ -99,6 +97,7 @@ export(contrib2DF)
 export(elemental_index)
 export(elementary_index)
 export(expand_classification)
+export(impute_prices)
 export(interact_classifications)
 export(is_aggregation_structure)
 export(is_chainable_index)
diff --git a/NEWS.md b/NEWS.md
index bd2e80b..39a23b5 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -34,11 +34,21 @@ aggregation structure that spans multiple hierarchical classifications.
 - Subscripting an index object with a matrix is now more flexible and can
 return/replace a list of index objects.
 
+- Added a new function `impute_prices()` to consolidate the different
+price-imputation functions.
+
 ## Bug fixes
 
 - Setting `contrib = FALSE` in `aggregate(index)` always returns an index with
 no contributions.
 
+- `set_weights()` no longer returns its result invisibly.
+
+## Deprecations
+
+- `carry_forward()`, `carry_backward()`, and `shadow_price()` are deprecated.
+Use `impute_prices()` instead.
+
 # piar 0.8.3
 
 ## Improvements
diff --git a/R/aggregate.piar_index.R b/R/aggregate.piar_index.R
index c48ad58..4246570 100644
--- a/R/aggregate.piar_index.R
+++ b/R/aggregate.piar_index.R
@@ -225,8 +225,8 @@ aggregate_index <- function(
       stop("'pias' and 'pias2' must represent the same aggregation structure")
     }
     if (
-      any(missing_weights(pias$weights) != missing_weights(pias2$weights)) &&
-        contrib
+      contrib &&
+        any(missing_weights(pias$weights) != missing_weights(pias2$weights))
     ) {
       stop(
         "any NA or zero weights must appear in both 'pias' and 'pias2' when",
@@ -261,7 +261,7 @@ aggregate_index <- function(
     lev <- unlist(drop_last(pias$levels), use.names = FALSE)
   }
 
-  piar_index(res$index, res$contrib, lev, x$time, chainable)
+  piar_index(res$index, res$contrib, lev, x$time, chainable = chainable)
 }
 
 aggregate_ <- function(
diff --git a/R/aggregation_structure-class.R b/R/aggregation_structure-class.R
index c5d6c3d..1234032 100644
--- a/R/aggregation_structure-class.R
+++ b/R/aggregation_structure-class.R
@@ -32,7 +32,7 @@ validate_pias_levels <- function(x) {
   if (anyDuplicated(lev)) {
     stop("cannot make an aggregation structure with duplicated levels")
   }
-  invisible(x)
+  x
 }
 
 validate_pias_structure <- function(x) {
@@ -49,7 +49,7 @@ validate_pias_structure <- function(x) {
       "hierachy"
     )
   }
-  invisible(x)
+  x
 }
 
 validate_pias_weights <- function(x) {
@@ -62,14 +62,13 @@ validate_pias_weights <- function(x) {
   if (any(x$weights < 0, na.rm = TRUE)) {
     stop("cannot make an aggregation structure with negative weights")
   }
-  invisible(x)
+  x
 }
 
 validate_piar_aggregation_structure <- function(x) {
   validate_pias_levels(x)
   validate_pias_structure(x)
   validate_pias_weights(x)
-  x
 }
 
 #---- Undocumented methods ----
diff --git a/R/as_index.R b/R/as_index.R
index 4a5b858..ddde2ec 100644
--- a/R/as_index.R
+++ b/R/as_index.R
@@ -79,9 +79,13 @@ as_index.matrix <- function(x, ..., chainable = TRUE, contrib = FALSE) {
   index <- as.numeric(x)
   dim(index) <- c(length(levels), length(periods))
 
-  contributions <- if (contrib) index2contrib(index, levels, periods)
-
-  piar_index(index, contributions, levels, periods, chainable)
+  piar_index(
+    index,
+    if (contrib) index2contrib(index, levels, periods),
+    levels,
+    periods,
+    chainable = chainable
+  )
 }
 
 #' @rdname as_index
diff --git a/R/contrib.R b/R/contrib.R
index 569afe2..c28e638 100644
--- a/R/contrib.R
+++ b/R/contrib.R
@@ -63,8 +63,12 @@
 #' @export contrib
 #' @family index methods
 contrib <- function(x, level = NULL, period = NULL, pad = 0) {
-  level <- match_levels(as.character(level %||% x$levels[1L]), x)
-  period <- match_time(as.character(period %||% x$time), x, several = TRUE)
+  level <- if (!is.null(level)) match_levels(as.character(level), x) else 1L
+  period <- if (!is.null(period)) {
+    match_time(as.character(period), x, several = TRUE)
+  } else {
+    seq_along(x$time)
+  }
   pad <- as.numeric(pad)
   if (length(pad) != 1L) {
     stop("'pad' must be a length 1 numeric value")
@@ -99,12 +103,16 @@ contrib <- function(x, level = NULL, period = NULL, pad = 0) {
 #' @rdname contrib
 #' @export
 contrib2DF <- function(x, levels = NULL, period = NULL) {
-  level <- match_levels(
-    as.character(levels %||% x$levels[1L]),
-    x,
-    several = TRUE
-  )
-  period <- match_time(as.character(period %||% x$time), x, several = TRUE)
+  level <- if (!is.null(levels)) {
+    match_levels(as.character(levels), x, several = TRUE)
+  } else {
+    1L
+  }
+  period <- if (!is.null(period)) {
+    match_time(as.character(period), x, several = TRUE)
+  } else {
+    seq_along(x$time)
+  }
   if (is.null(x$contrib)) {
     return(
       data.frame(
@@ -141,8 +149,12 @@ contrib2DF <- function(x, levels = NULL, period = NULL) {
 #' @rdname contrib
 #' @export
 `contrib<-` <- function(x, level = NULL, period = NULL, value) {
-  level <- match_levels(as.character(level %||% x$levels[1L]), x)
-  period <- match_time(as.character(period %||% x$time), x, several = TRUE)
+  level <- if (!is.null(level)) match_levels(as.character(level), x) else 1L
+  period <- if (!is.null(period)) {
+    match_time(as.character(period), x, several = TRUE)
+  } else {
+    seq_along(x$time)
+  }
 
   if (is.null(x$contrib)) {
     x$contrib <- contrib_skeleton(x$levels, x$time)
@@ -156,7 +168,12 @@ contrib2DF <- function(x, levels = NULL, period = NULL) {
       "number of items to replace is not a multiple of replacement length"
     )
   }
-  value[] <- as.numeric(value)
+  value <- matrix(
+    as.numeric(value),
+    nrow = nrow(value),
+    ncol = ncol(value),
+    dimnames = dimnames(value)
+  )
 
   products <- if (nrow(value) > 0L) {
     if (is.null(rownames(value))) {
@@ -173,7 +190,7 @@ contrib2DF <- function(x, levels = NULL, period = NULL) {
     names(con) <- products
     x$contrib[level, t] <- list(con)
   }
-  validate_piar_index(x)
+  validate_contrib(x)
 }
 
 #' @rdname contrib
@@ -184,5 +201,5 @@ set_contrib <- `contrib<-`
 #' @export
 set_contrib_from_index <- function(x) {
   x$contrib <- index2contrib(x$index, x$levels, x$time)
-  x
+  validate_contrib(x)
 }
diff --git a/R/elementary_index.R b/R/elementary_index.R
index e44c408..4ac7475 100644
--- a/R/elementary_index.R
+++ b/R/elementary_index.R
@@ -82,8 +82,9 @@
 #'   a Paasche index). Other values are possible; see
 #'   [gpindex::generalized_mean()] for details.
 #' @param ... Further arguments passed to or used by methods.
-#' @param formula A two-sided formula with price relatives on the left-hand
-#'   side, and time periods and elementary aggregates (in that order) on the
+#' @param formula A two-sided formula, or something that can be coerced into
+#'   one, with price relatives on the left-hand
+#'   side and time periods and elementary aggregates (in that order) on the
 #'   right-hand side.
 #'
 #' @returns
@@ -196,9 +197,6 @@ elementary_index.numeric <- function(
   }
   period <- as.factor(period %||% gl(1, length(x)))
   ea <- as.factor(ea %||% gl(1, length(x)))
-  ea_by_period <- period:ea
-  time <- levels(period)
-  levels <- levels(ea)
 
   if (different_length(x, period, ea, weights)) {
     stop("input vectors must be the same length")
@@ -206,6 +204,9 @@ elementary_index.numeric <- function(
   if (any(x <= 0, na.rm = TRUE)) {
     stop("all elements of 'x' must be strictly positive")
   }
+  ea_by_period <- period:ea
+  time <- levels(period)
+  levels <- levels(ea)
 
   if (contrib) {
     if (!is.null(product)) {
@@ -219,11 +220,7 @@ elementary_index.numeric <- function(
   }
 
   x <- split(x, ea_by_period)
-  if (is.null(weights)) {
-    weights <- list(NULL)
-  } else {
-    weights <- split(weights, ea_by_period)
-  }
+  weights <- if (is.null(weights)) list(NULL) else split(weights, ea_by_period)
 
   index <- mapply(
     gpindex::generalized_mean(r),
@@ -243,9 +240,9 @@ elementary_index.numeric <- function(
       USE.NAMES = FALSE
     )
     dim(contributions) <- c(nlevels(ea), nlevels(period))
-    piar_index(index, contributions, levels, time, chainable)
+    piar_index(index, contributions, levels, time, chainable = chainable)
   } else {
-    piar_index(index, NULL, levels, time, chainable)
+    piar_index(index, NULL, levels, time, chainable = chainable)
   }
 }
 
diff --git a/R/extract.piar_index.R b/R/extract.piar_index.R
index d578e1d..c6c5dff 100644
--- a/R/extract.piar_index.R
+++ b/R/extract.piar_index.R
@@ -69,7 +69,7 @@ extract_index <- function(x, levels, periods) {
   }
   x$levels <- x$levels[levels]
   x$time <- x$time[periods]
-  validate_piar_index(x)
+  validate_index_structure(x)
 }
 
 #' @rdname sub-.piar_index
@@ -86,18 +86,18 @@ extract_index <- function(x, levels, periods) {
     if (is_index(value)) {
       value <- list(value)
     }
-    if (is.list(value)) {
-      x <- replace_matrix_list(x, i, value)
+    x <- if (is.list(value)) {
+      replace_matrix_list(x, i, value)
     } else {
-      x <- replace_matrix_numeric(x, i, value)
+      replace_matrix_numeric(x, i, value)
     }
   } else {
     levels <- subscript_index(x$levels, i)
     periods <- subscript_index(x$time, j)
-    if (is_index(value)) {
-      x <- replace_index(x, levels, periods, value)
+    x <- if (is_index(value)) {
+      replace_index(x, levels, periods, value)
     } else {
-      x <- replace_numeric(x, levels, periods, value)
+      replace_numeric(x, levels, periods, value)
     }
   }
   # Replacement value should be validated; e.g., x[1] <- -1.
@@ -151,9 +151,8 @@ replace_matrix_list <- function(x, i, value) {
     stop("'value' must be a list of indexes with one level and time period")
   }
   # Make `value` the same length as replacement to avoid two warnings.
-  value <- rep_len(value, nrow(i))
-  index <- vapply(value, \(x) x$index, numeric(1L))
-  contributions <- unlist(lapply(value, \(x) x$contrib), recursive = FALSE)
+  index <- rep_len(vapply(value, \(x) x$index[[1]], numeric(1L)), nrow(i))
+  contributions <- rep_len(lapply(value, \(x) x$contrib[[1]]), nrow(i))
   has_contrib <- any(vapply(contributions, Negate(is.null), logical(1L)))
   # It's possible that only some elements of `value` have contributions.
   if (has_contrib || !is.null(x$contrib)) {
diff --git a/R/head.piar_index.R b/R/head.piar_index.R
index c31a1f8..d0e6e70 100644
--- a/R/head.piar_index.R
+++ b/R/head.piar_index.R
@@ -31,7 +31,7 @@ head.piar_index <- function(x, n = 6L, ...) {
   if (!is.na(n[2L])) {
     x$time <- head(x$time, n[2L])
   }
-  validate_piar_index(x)
+  validate_index_structure(x)
 }
 
 #' @rdname head.piar_index
@@ -47,5 +47,5 @@ tail.piar_index <- function(x, n = 6L, ...) {
   if (!is.na(n[2L])) {
     x$time <- tail(x$time, n[2L])
   }
-  validate_piar_index(x)
+  validate_index_structure(x)
 }
diff --git a/R/helpers.R b/R/helpers.R
index 267e527..fd7d33a 100644
--- a/R/helpers.R
+++ b/R/helpers.R
@@ -14,7 +14,7 @@ missing_weights <- function(x) {
 }
 
 missing_names <- function(x) {
-  anyNA(x) || any(x == "")
+  anyNA(x) || !all(nzchar(x))
 }
 
 last <- function(x) {
@@ -104,6 +104,7 @@ different_length <- function(...) {
 }
 
 formula_vars <- function(formula, x, n = 2L) {
+  formula <- stats::as.formula(formula)
   if (length(formula) != 3L) {
     stop("'formula' must have a left-hand and right-hand side")
   }
@@ -195,19 +196,10 @@ index_skeleton <- function(levels, time) {
   matrix(NA_real_, length(levels), length(time))
 }
 
-empty_contrib <- function(x) {
-  res <- rep.int(list(numeric(0L)), length(x))
-  list(res)
-}
-
 contrib_skeleton <- function(levels, time) {
   matrix(list(numeric(0L)), length(levels), length(time))
 }
 
-has_contrib <- function(x) {
-  any(lengths(x$contrib) > 0L)
-}
-
 # Backport Reduce and %||%
 # TODO: Remove once min R version gets bumped.
 if (getRversion() < "4.4.0") {
diff --git a/R/impute_prices.R b/R/impute_prices.R
index 43a11af..1969612 100644
--- a/R/impute_prices.R
+++ b/R/impute_prices.R
@@ -1,6 +1,7 @@
 #' Impute missing prices
 #'
-#' Impute missing prices using the carry forward or shadow price method.
+#' Impute missing prices using the carry forward/backward or the
+#' self-correcting overall mean method.
 #'
 #' The carry forward method replaces a missing price for a product by the price
 #' for the same product in the previous period. It tends to push an index value
@@ -8,27 +9,32 @@
 #' (2020). The carry backwards method does the opposite, but this is rarely
 #' used in practice.
 #'
-#' The shadow price method recursively imputes a missing price by the value of
+#' The self-correcting overall mean method recursively imputes a missing price
+#' by the value of
 #' the price for the same product in the previous period multiplied by the
 #' value of the period-over-period elementary index for the elementary aggregate
 #' to which that product belongs. This requires computing and aggregating an
 #' index (according to `pias`, unless `pias` is not supplied) for
-#' each `period`, and so these imputations can take a while. The index
+#' each `period`. The index
 #' values used to do the imputations are not returned because the index needs
-#' to be recalculated to get correct percent-change contributions.
-#'
-#' Shadow price imputation is referred to as self-correcting overall mean
-#' imputation in chapter 6 of the CPI manual (2020). It is identical to simply
+#' to be recalculated to get correct percent-change contributions. It is
+#' identical to simply
 #' excluding missing price relatives in the index calculation, except in the
 #' period that a missing product returns. For this reason care is needed when
 #' using this method. It is sensitive to the assumption that a product does not
 #' change over time, and in some cases it is safer to simply omit the missing
 #' price relatives instead of imputing the missing prices.
 #'
+#' Imputation works slightly differently depending on whether data are in a long
+#' or wide format. When `x` is a two-column of matrix of current and back prices
+#' (in that order), then imputation is done separately on the current price
+#' at a point in time and the back price at the next point in time. When `x` is
+#' a numeric vector then these two prices are necessarily the same.
+#'
 #' @name impute_prices
-#' @aliases impute_prices
-#' @param x Either a numeric vector (or something that can be coerced into one)
-#'   or data frame of prices.
+#' @param x Either a numeric vector (or something that can be coerced into one),
+#'   a data frame of prices, or a two-column matrix of current prices
+#'   and back prices (in that order).
 #' @param period A factor, or something that can be coerced into one, giving
 #'   the time period associated with each price in `x`. The ordering of time
 #'   periods follows of the levels of `period`, to agree with
@@ -36,7 +42,9 @@
 #' @param product A factor, or something that can be coerced into one, giving
 #'   the product associated with each price in `x`.
 #' @param ea A factor, or something that can be coerced into one, giving the
-#'   elementary aggregate associated with each price in `x`.
+#'   elementary aggregate associated with each price in `x`. This is evaluated
+#'   in `x` for the data frame method. The default pools all data into one
+#'   elementary aggregate.
 #' @param pias A price index aggregation structure, or something that can be
 #'   coerced into one, as made with [aggregation_structure()]. The default
 #'   imputes from elementary indexes only (i.e., not recursively).
@@ -44,23 +52,23 @@
 #'   product weights), or something that can be coerced into one. The default is
 #'   to give each price equal weight. This is evaluated in `x` for the data
 #'   frame method.
-#' @param r1 Order of the generalized-mean price index used to calculate the
-#'   elementary price indexes: 0 for a geometric index (the default), 1 for an
-#'   arithmetic index, or -1 for a harmonic index. Other values are possible;
-#'   see [gpindex::generalized_mean()] for details.
-#' @param r2 Order of the generalized-mean price index used to aggregate the
-#'   elementary price indexes: 0 for a geometric index, 1 for an arithmetic
-#'   index (the default), or -1 for a harmonic index. Other values are possible;
+#' @param r A pair of numeric values. The first gives the order of the
+#'   generalized-mean price index used to calculate the
+#'   elementary price indexes, defaulting to a geometric index. The second
+#'   gives the order of the generalized-mean price index used to aggregate the
+#'   elementary price indexes, defaulting to an arithmetic index. Other values
+#'   are possible;
 #'   see [gpindex::generalized_mean()] for details.
-#' @param formula A two-sided formula with prices on the left-hand
-#'   side. For `carry_forward()` and `carry_backward()`, the right-hand side
-#'   should have time periods and products (in that order); for
-#'   `shadow_price()`, the right-hand side should have time period, products,
-#'   and elementary aggregates (in that order).
+#' @param formula A two-sided formula, or something that can be coerced into
+#'   one, with prices on the left-hand
+#'   side and time periods and products on the right-hand side (in that order).
+#' @param method Name of the imputation method, one of `"overall-mean"`,
+#'   `"carry-forward"`, or `"carry-backward"`.
 #' @param ... Further arguments passed to or used by methods.
 #'
 #' @returns
-#' A numeric vector of prices with missing values replaced (where possible).
+#' A numeric vector or matrix of prices with missing values replaced
+#' (where possible).
 #'
 #' @seealso
 #' [price_relative()] for making price relatives for the
@@ -79,123 +87,156 @@
 #'   ea = rep(letters[1:2], 4)
 #' )
 #'
-#' carry_forward(prices, price ~ period + product)
+#' impute_prices(prices, price ~ period + product, method = "carry-forward")
 #'
-#' shadow_price(prices, price ~ period + product + ea)
+#' impute_prices(
+#'   prices,
+#'   price ~ period + product,
+#'   ea = ea,
+#'   method = "overall-mean"
+#' )
+#'
+#' # Can also be done with current price-back price formulation.
+#' prices$back_price <- with(
+#'   prices,
+#'   price[gpindex::back_period(period, product)]
+#' )
+#'
+#' impute_prices(
+#'   prices,
+#'   cbind(price, back_price) ~ period + product,
+#'   ea = ea,
+#'   method = "overall-mean"
+#' )
 #'
 #' @export
-shadow_price <- function(x, ...) {
-  UseMethod("shadow_price")
+impute_prices <- function(x, ...) {
+  UseMethod("impute_prices")
 }
 
 #' @rdname impute_prices
 #' @export
-shadow_price.default <- function(
+impute_prices.default <- function(x, ...) {
+  impute_prices(as.numeric(x), ...)
+}
+
+#' @rdname impute_prices
+#' @export
+impute_prices.matrix <- function(
   x,
   period,
   product,
-  ea,
   ...,
-  pias = NULL,
+  ea = NULL,
   weights = NULL,
-  r1 = 0,
-  r2 = 1
+  pias = NULL,
+  r = c(0, 1),
+  method = c("overall-mean", "carry-forward")
 ) {
   # This is mostly a combination of gpindex::back_period() and aggregate()
   # it just does it period-by-period and keeps track of prices to impute.
   chkDots(...)
-  x <- as.numeric(x)
+  method <- match.arg(method)
   period <- as.factor(period)
   product <- as.factor(product)
   attributes(product) <- NULL
-  ea <- as.factor(ea)
+  if (!is.null(ea)) {
+    ea <- as.factor(ea)
+  }
   if (!is.null(weights)) {
     weights <- as.numeric(weights)
   }
 
-  if (different_length(x, period, product, ea, weights)) {
+  if (different_length(x[, 1L], period, product, ea, weights)) {
     stop("input vectors must be the same length")
   }
   if (nlevels(period) == 0L) {
-    return(rep.int(NA_integer_, length(period)))
+    return(matrix(NA_real_, nrow = length(period), ncol = 2))
   }
 
-  res <- split(x, period)
+  res <- split.data.frame(x, period)
   product <- split(product, period)
   if (duplicate_products(product)) {
     warning("there are duplicated period-product pairs")
   }
-  ea <- split(ea, period)
-  if (is.null(weights)) {
-    w <- rep.int(list(NULL), nlevels(period))
-  } else {
-    w <- split(weights, period)
+  if (!is.null(ea)) {
+    ea <- split(ea, period)
+  }
+  if (!is.null(weights)) {
+    weights <- split(weights, period)
   }
   if (!is.null(pias)) {
     pias <- as_aggregation_structure(pias)
   }
-  for (t in seq_along(res)[-1L]) {
-    # Calculate relatives.
-    matches <- match(product[[t]], product[[t - 1L]], incomparables = NA)
-    back_price <- res[[t - 1L]][matches]
-    price <- res[[t]]
-    # Calculate indexes.
-    epr <- elementary_index(
-      price / back_price,
-      period = gl(1, length(price)),
-      ea = ea[[t]],
-      weights = w[[t]],
-      na.rm = TRUE,
-      r = r1
-    )
-    if (!is.null(pias)) {
-      epr <- aggregate(epr, pias, na.rm = TRUE, r = r2)
-      pias <- update(pias, epr, r = r2)
+  for (t in seq_along(res)) {
+    impute <- which(is.na(res[[t]][, 1L]))
+    if (method == "overall-mean") {
+      index <- elementary_index(
+        res[[t]][, 1L] / res[[t]][, 2L],
+        ea = ea[[t]],
+        weights = weights[[t]],
+        na.rm = TRUE,
+        r = r[1L]
+      )
+      if (!is.null(pias)) {
+        index <- aggregate(index, pias, na.rm = TRUE, r = r[2L])
+        pias <- update(pias, index, r = r[2L])
+      }
+      eas <- if (!is.null(ea)) {
+        match(as.character(ea[[t]][impute]), index$levels)
+      } else {
+        1L
+      }
+      res[[t]][impute, 1L] <- as.numeric(index)[eas] * res[[t]][impute, 2L]
+    } else {
+      res[[t]][impute, 1L] <- res[[t]][impute, 2L]
+    }
+    if (t < length(res)) {
+      impute2 <- which(is.na(res[[t + 1L]][, 2L]))
+      matches <- match(product[[t + 1L]][impute2], product[[t]][impute])
+      res[[t + 1L]][impute2, 2L] <- res[[t]][impute, 1L][matches]
     }
-    # Add shadow prices to 'x'.
-    impute <- is.na(price)
-    eas <- match(as.character(ea[[t]][impute]), epr$levels)
-    res[[t]][impute] <- epr$index[, 1L][eas] * back_price[impute]
   }
-  unsplit(res, period)
-}
-
-#' @rdname impute_prices
-#' @export
-shadow_price.data.frame <- function(x, formula, ..., weights = NULL) {
-  vars <- formula_vars(formula, x, 3L)
-  weights <- eval(substitute(weights), x, parent.frame())
-
-  shadow_price(
-    vars[[1L]],
-    period = vars[[2L]],
-    product = vars[[3L]],
-    ea = vars[[4L]],
-    weights = weights,
-    ...
-  )
-}
-
-#' @rdname impute_prices
-#' @export
-carry_forward <- function(x, ...) {
-  UseMethod("carry_forward")
+  # unsplit() doesn't work nicely with matrices.
+  split(x, period) <- res
+  x
 }
 
 #' @rdname impute_prices
 #' @export
-carry_forward.default <- function(x, period, product, ...) {
+impute_prices.numeric <- function(
+  x,
+  period,
+  product,
+  ...,
+  ea = NULL,
+  weights = NULL,
+  pias = NULL,
+  r = c(0, 1),
+  method = c("overall-mean", "carry-forward", "carry-backward")
+) {
+  # This is mostly a combination of gpindex::back_period() and aggregate()
+  # it just does it period-by-period and keeps track of prices to impute.
   chkDots(...)
-  x <- as.numeric(x)
+  method <- match.arg(method)
   period <- as.factor(period)
+  if (method == "carry-backward") {
+    period <- factor(period, rev(levels(period)))
+  }
   product <- as.factor(product)
   attributes(product) <- NULL
+  if (!is.null(ea)) {
+    ea <- as.factor(ea)
+  }
+  if (!is.null(weights)) {
+    weights <- as.numeric(weights)
+  }
 
-  if (different_length(x, period, product)) {
+  if (different_length(x, period, product, ea, weights)) {
     stop("input vectors must be the same length")
   }
   if (nlevels(period) == 0L) {
-    return(rep.int(NA_integer_, length(period)))
+    return(rep.int(NA_real_, length(period)))
   }
 
   res <- split(x, period)
@@ -203,47 +244,90 @@ carry_forward.default <- function(x, period, product, ...) {
   if (duplicate_products(product)) {
     warning("there are duplicated period-product pairs")
   }
+  if (!is.null(ea)) {
+    ea <- split(ea, period)
+  }
+  if (!is.null(weights)) {
+    weights <- split(weights, period)
+  }
+  if (!is.null(pias)) {
+    pias <- as_aggregation_structure(pias)
+  }
   for (t in seq_along(res)[-1L]) {
-    impute <- is.na(res[[t]])
-    matches <- match(
-      product[[t]][impute],
-      product[[t - 1L]],
-      incomparables = NA
-    )
-    res[[t]][impute] <- res[[t - 1L]][matches]
+    impute <- which(is.na(res[[t]]))
+    if (method == "overall-mean") {
+      matches <- match(product[[t]], product[[t - 1L]], incomparables = NA)
+      back_price <- res[[t - 1L]][matches]
+
+      index <- elementary_index(
+        res[[t]] / back_price,
+        ea = ea[[t]],
+        weights = weights[[t]],
+        na.rm = TRUE,
+        r = r[1L]
+      )
+      if (!is.null(pias)) {
+        index <- aggregate(index, pias, na.rm = TRUE, r = r[2L])
+        pias <- update(pias, index, r = r[2L])
+      }
+      eas <- if (!is.null(ea)) {
+        match(as.character(ea[[t]][impute]), index$levels)
+      } else {
+        1L
+      }
+      res[[t]][impute] <- as.numeric(index)[eas] * back_price[impute]
+    } else {
+      matches <- match(
+        product[[t]][impute],
+        product[[t - 1L]],
+        incomparables = NA
+      )
+      res[[t]][impute] <- res[[t - 1L]][matches]
+    }
   }
   unsplit(res, period)
 }
 
 #' @rdname impute_prices
 #' @export
-carry_forward.data.frame <- function(x, formula, ...) {
-  chkDots(...)
-  vars <- formula_vars(formula, x)
+impute_prices.data.frame <- function(
+  x,
+  formula,
+  ...,
+  ea = NULL,
+  weights = NULL
+) {
+  vars <- formula_vars(formula, x, 2L)
+  ea <- eval(substitute(ea), x, parent.frame())
+  weights <- eval(substitute(weights), x, parent.frame())
 
-  carry_forward(vars[[1L]], period = vars[[2L]], product = vars[[3L]])
+  impute_prices(
+    vars[[1L]],
+    period = vars[[2L]],
+    product = vars[[3L]],
+    ea = ea,
+    weights = weights,
+    ...
+  )
 }
 
 #' @rdname impute_prices
 #' @export
-carry_backward <- function(x, ...) {
-  UseMethod("carry_backward")
+carry_forward <- function(x, ...) {
+  warning("'carry_forward() is deprecated; use 'impute_prices()' instead")
+  impute_prices(x, method = "carry-forward", ...)
 }
 
 #' @rdname impute_prices
 #' @export
-carry_backward.default <- function(x, period, product, ...) {
-  chkDots(...)
-  period <- as.factor(period)
-  levels <- rev(levels(period))
-  carry_forward(x, period = factor(period, levels), product = product)
+carry_backward <- function(x, ...) {
+  warning("'carry_backward() is deprecated; use 'impute_prices()' instead")
+  impute_prices(x, method = "carry-backward", ...)
 }
 
 #' @rdname impute_prices
 #' @export
-carry_backward.data.frame <- function(x, formula, ...) {
-  chkDots(...)
-  vars <- formula_vars(formula, x)
-
-  carry_backward(vars[[1L]], period = vars[[2L]], product = vars[[3L]])
+shadow_price <- function(x, ...) {
+  warning("'shadow_price() is deprecated; use 'impute_prices()' instead")
+  impute_prices(x, method = "overall-mean", ...)
 }
diff --git a/R/index-class.R b/R/index-class.R
index d7224f6..e6d7398 100644
--- a/R/index-class.R
+++ b/R/index-class.R
@@ -39,8 +39,10 @@ new_piar_index <- function(index, contrib, levels, time, chainable) {
 }
 
 piar_index <- function(index, contrib, levels, time, chainable) {
-  index <- unname(index)
-  contrib <- unname(contrib)
+  index <- structure(as.numeric(index), dim = dim(index))
+  if (!is.null(contrib)) {
+    contrib <- structure(as.list(contrib), dim = dim(contrib))
+  }
   levels <- as.character(levels)
   time <- as.character(time)
   validate_piar_index(
@@ -59,7 +61,10 @@ validate_levels <- function(x) {
   if (anyDuplicated(x$levels)) {
     stop("cannot make an index with duplicate levels")
   }
-  invisible(x)
+  if (nrow(x$index) != length(x$levels)) {
+    stop("number of levels does not agree with number of index values")
+  }
+  x
 }
 
 validate_time <- function(x) {
@@ -72,20 +77,17 @@ validate_time <- function(x) {
   if (anyDuplicated(x$time)) {
     stop("cannot make an index with duplicate time periods")
   }
-  invisible(x)
-}
-
-validate_index_values <- function(x) {
   if (ncol(x$index) != length(x$time)) {
     stop("number of time periods does not agree with number of index values")
   }
-  if (nrow(x$index) != length(x$levels)) {
-    stop("number of levels does not agree with number of index values")
-  }
+  x
+}
+
+validate_index_values <- function(x) {
   if (any(x$index <= 0, na.rm = TRUE)) {
     stop("cannot make an index with non-positive values")
   }
-  invisible(x)
+  x
 }
 
 validate_contrib <- function(x) {
@@ -97,15 +99,18 @@ validate_contrib <- function(x) {
       stop("number of levels does not agree with number of contributions")
     }
   }
-  invisible(x)
+  x
 }
 
-validate_piar_index <- function(x) {
+validate_index_structure <- function(x) {
   validate_levels(x)
   validate_time(x)
-  validate_index_values(x)
   validate_contrib(x)
-  x
+}
+
+validate_piar_index <- function(x) {
+  validate_index_values(x)
+  validate_index_structure(x)
 }
 
 #---- Undocumented methods ----
diff --git a/R/levels.piar_index.R b/R/levels.piar_index.R
index d4b7ee3..33ed67e 100644
--- a/R/levels.piar_index.R
+++ b/R/levels.piar_index.R
@@ -22,7 +22,7 @@ levels.piar_index <- function(x) {
 #' @export
 `levels<-.piar_index` <- function(x, value) {
   x$levels <- as.character(value)
-  validate_piar_index(x)
+  validate_levels(x)
 }
 
 #' @rdname levels.piar_index
diff --git a/R/mean.piar_index.R b/R/mean.piar_index.R
index 407aa31..336b084 100644
--- a/R/mean.piar_index.R
+++ b/R/mean.piar_index.R
@@ -134,8 +134,8 @@ mean_index <- function(
   }
 
   window <- as.integer(window %||% ntime(x))
-  if (length(window) > 1L || window < 1L) {
-    stop("'window' must be a positive length 1 integer")
+  if (window < 1L) {
+    stop("'window' must be a positive integer")
   }
   if (window > ntime(x)) {
     stop("'x' must have at least 'window' time periods")
@@ -180,6 +180,6 @@ mean_index <- function(
     do.call(cbind, contrib),
     x$levels,
     periods,
-    chainable
+    chainable = chainable
   )
 }
diff --git a/R/merge.piar_index.R b/R/merge.piar_index.R
index 19fd8d7..122f679 100644
--- a/R/merge.piar_index.R
+++ b/R/merge.piar_index.R
@@ -32,7 +32,7 @@
 merge.chainable_piar_index <- function(x, y, ...) {
   y <- as_index(y, chainable = TRUE)
   res <- NextMethod("merge")
-  new_piar_index(res$index, res$contrib, res$levels, res$time, TRUE)
+  new_piar_index(res$index, res$contrib, res$levels, res$time, chainable = TRUE)
 }
 
 #' @rdname merge.piar_index
@@ -40,7 +40,13 @@ merge.chainable_piar_index <- function(x, y, ...) {
 merge.direct_piar_index <- function(x, y, ...) {
   y <- as_index(y, chainable = FALSE)
   res <- NextMethod("merge")
-  new_piar_index(res$index, res$contrib, res$levels, res$time, FALSE)
+  new_piar_index(
+    res$index,
+    res$contrib,
+    res$levels,
+    res$time,
+    chainable = FALSE
+  )
 }
 
 #' @export
diff --git a/R/price_relative.R b/R/price_relative.R
index 8176ea8..8099de3 100644
--- a/R/price_relative.R
+++ b/R/price_relative.R
@@ -12,8 +12,9 @@
 #' @param product A factor, or something that can be coerced into one, that
 #'   gives the corresponding product identifier for each element in `x`.
 #' @param ... Further arguments passed to or used by methods.
-#' @param formula A two-sided formula with prices on the left-hand
-#'   side, and time periods and products (in that order) on the
+#' @param formula A two-sided formula, or something that can be coerced into
+#'   one, with prices on the left-hand
+#'   side and time periods and products (in that order) on the
 #'   right-hand side.
 #'
 #' @returns
diff --git a/R/split_classification.R b/R/split_classification.R
index 4523359..15d19db 100644
--- a/R/split_classification.R
+++ b/R/split_classification.R
@@ -12,7 +12,7 @@
 #' @param sep A character used to delineate levels in `x` in the result. The
 #'   default separates levels by `"."`.
 #' @param pad A string used to pad the shorter labels for an unbalanced
-#'   classification. The default pads with NA.
+#'   classification. The default pads with `NA`.
 #'
 #' @returns
 #' A list with a entry for each level in `x` giving the "digits" that
diff --git a/R/stack.piar_index.R b/R/stack.piar_index.R
index 2d95165..7c2ce9a 100644
--- a/R/stack.piar_index.R
+++ b/R/stack.piar_index.R
@@ -46,7 +46,7 @@
 stack.chainable_piar_index <- function(x, y, ...) {
   y <- as_index(y, chainable = TRUE)
   res <- NextMethod("stack")
-  new_piar_index(res$index, res$contrib, res$levels, res$time, TRUE)
+  new_piar_index(res$index, res$contrib, res$levels, res$time, chainable = TRUE)
 }
 
 #' @rdname stack.piar_index
@@ -54,7 +54,13 @@ stack.chainable_piar_index <- function(x, y, ...) {
 stack.direct_piar_index <- function(x, y, ...) {
   y <- as_index(y, chainable = FALSE)
   res <- NextMethod("stack")
-  new_piar_index(res$index, res$contrib, res$levels, res$time, FALSE)
+  new_piar_index(
+    res$index,
+    res$contrib,
+    res$levels,
+    res$time,
+    chainable = FALSE
+  )
 }
 
 #' @export
@@ -98,6 +104,7 @@ unstack.direct_piar_index <- function(x, ...) {
 
 #' @export
 unstack.piar_index <- function(x, ..., chainable) {
+  chkDots(...)
   res <- vector("list", ntime(x))
   names(res) <- x$time
   for (t in seq_along(res)) {
@@ -106,7 +113,7 @@ unstack.piar_index <- function(x, ..., chainable) {
       x$contrib[, t, drop = FALSE],
       x$levels,
       x$time[t],
-      chainable
+      chainable = chainable
     )
   }
   res
diff --git a/R/time.piar_index.R b/R/time.piar_index.R
index 46bf0f6..4187829 100644
--- a/R/time.piar_index.R
+++ b/R/time.piar_index.R
@@ -34,7 +34,7 @@ time.piar_index <- function(x, ...) {
 #' @export
 `time<-.piar_index` <- function(x, value) {
   x$time <- as.character(value)
-  validate_piar_index(x)
+  validate_time(x)
 }
 
 #' @rdname time.piar_index
@@ -46,7 +46,7 @@ set_time <- `time<-`
 #' @export
 start.piar_index <- function(x, ...) {
   chkDots(...)
-  x$time[1L]
+  x$time[[1L]]
 }
 
 #' @rdname time.piar_index
diff --git a/R/update.aggregation_structure.R b/R/update.aggregation_structure.R
index 69d1f87..b78c1c5 100644
--- a/R/update.aggregation_structure.R
+++ b/R/update.aggregation_structure.R
@@ -59,7 +59,11 @@ update.piar_aggregation_structure <- function(
   chkDots(...)
   price_update <- gpindex::factor_weights(r)
   index <- chain(as_index(index))
-  period <- match_time(as.character(period %||% last(index$time)), index)
+  period <- if (!is.null(period)) {
+    match_time(as.character(period), index)
+  } else {
+    ntime(index)
+  }
   eas <- match_eas(object, index)
   if (anyNA(eas)) {
     warning("not all weights in 'object' have a corresponding index value")
diff --git a/R/window.piar_index.R b/R/window.piar_index.R
index d6c3f2f..11be575 100644
--- a/R/window.piar_index.R
+++ b/R/window.piar_index.R
@@ -44,8 +44,8 @@ window.piar_index <- function(x, start = NULL, end = NULL, ...) {
 #' Get the indexes for a window of time periods
 #' @noRd
 index_window <- function(x, start, end) {
-  start <- match_time(as.character(start %||% x$time[1L]), x)
-  end <- match_time(as.character(end %||% last(x$time)), x)
+  start <- if (!is.null(start)) match_time(as.character(start), x) else 1L
+  end <- if (!is.null(end)) match_time(as.character(end), x) else ntime(x)
 
   if (start > end) {
     stop("'start' must refer to a time period before 'end'")
diff --git a/man/elementary_index.Rd b/man/elementary_index.Rd
index 9327c7f..340c106 100644
--- a/man/elementary_index.Rd
+++ b/man/elementary_index.Rd
@@ -73,8 +73,9 @@ averaging indexes over subperiods), or -1 for a harmonic index (usually for
 a Paasche index). Other values are possible; see
 \code{\link[gpindex:generalized_mean]{gpindex::generalized_mean()}} for details.}
 
-\item{formula}{A two-sided formula with price relatives on the left-hand
-side, and time periods and elementary aggregates (in that order) on the
+\item{formula}{A two-sided formula, or something that can be coerced into
+one, with price relatives on the left-hand
+side and time periods and elementary aggregates (in that order) on the
 right-hand side.}
 }
 \value{
diff --git a/man/impute_prices.Rd b/man/impute_prices.Rd
index f20647f..6253cba 100644
--- a/man/impute_prices.Rd
+++ b/man/impute_prices.Rd
@@ -2,48 +2,55 @@
 % Please edit documentation in R/impute_prices.R
 \name{impute_prices}
 \alias{impute_prices}
-\alias{shadow_price}
-\alias{shadow_price.default}
-\alias{shadow_price.data.frame}
+\alias{impute_prices.default}
+\alias{impute_prices.matrix}
+\alias{impute_prices.numeric}
+\alias{impute_prices.data.frame}
 \alias{carry_forward}
-\alias{carry_forward.default}
-\alias{carry_forward.data.frame}
 \alias{carry_backward}
-\alias{carry_backward.default}
-\alias{carry_backward.data.frame}
+\alias{shadow_price}
 \title{Impute missing prices}
 \usage{
-shadow_price(x, ...)
+impute_prices(x, ...)
+
+\method{impute_prices}{default}(x, ...)
 
-\method{shadow_price}{default}(
+\method{impute_prices}{matrix}(
   x,
   period,
   product,
-  ea,
   ...,
+  ea = NULL,
+  weights = NULL,
   pias = NULL,
+  r = c(0, 1),
+  method = c("overall-mean", "carry-forward")
+)
+
+\method{impute_prices}{numeric}(
+  x,
+  period,
+  product,
+  ...,
+  ea = NULL,
   weights = NULL,
-  r1 = 0,
-  r2 = 1
+  pias = NULL,
+  r = c(0, 1),
+  method = c("overall-mean", "carry-forward", "carry-backward")
 )
 
-\method{shadow_price}{data.frame}(x, formula, ..., weights = NULL)
+\method{impute_prices}{data.frame}(x, formula, ..., ea = NULL, weights = NULL)
 
 carry_forward(x, ...)
 
-\method{carry_forward}{default}(x, period, product, ...)
-
-\method{carry_forward}{data.frame}(x, formula, ...)
-
 carry_backward(x, ...)
 
-\method{carry_backward}{default}(x, period, product, ...)
-
-\method{carry_backward}{data.frame}(x, formula, ...)
+shadow_price(x, ...)
 }
 \arguments{
-\item{x}{Either a numeric vector (or something that can be coerced into one)
-or data frame of prices.}
+\item{x}{Either a numeric vector (or something that can be coerced into one),
+a data frame of prices, or a two-column matrix of current prices
+and back prices (in that order).}
 
 \item{...}{Further arguments passed to or used by methods.}
 
@@ -56,38 +63,41 @@ periods follows of the levels of \code{period}, to agree with
 the product associated with each price in \code{x}.}
 
 \item{ea}{A factor, or something that can be coerced into one, giving the
-elementary aggregate associated with each price in \code{x}.}
-
-\item{pias}{A price index aggregation structure, or something that can be
-coerced into one, as made with \code{\link[=aggregation_structure]{aggregation_structure()}}. The default
-imputes from elementary indexes only (i.e., not recursively).}
+elementary aggregate associated with each price in \code{x}. This is evaluated
+in \code{x} for the data frame method. The default pools all data into one
+elementary aggregate.}
 
 \item{weights}{A numeric vector of weights for the prices in \code{x} (i.e.,
 product weights), or something that can be coerced into one. The default is
 to give each price equal weight. This is evaluated in \code{x} for the data
 frame method.}
 
-\item{r1}{Order of the generalized-mean price index used to calculate the
-elementary price indexes: 0 for a geometric index (the default), 1 for an
-arithmetic index, or -1 for a harmonic index. Other values are possible;
-see \code{\link[gpindex:generalized_mean]{gpindex::generalized_mean()}} for details.}
+\item{pias}{A price index aggregation structure, or something that can be
+coerced into one, as made with \code{\link[=aggregation_structure]{aggregation_structure()}}. The default
+imputes from elementary indexes only (i.e., not recursively).}
 
-\item{r2}{Order of the generalized-mean price index used to aggregate the
-elementary price indexes: 0 for a geometric index, 1 for an arithmetic
-index (the default), or -1 for a harmonic index. Other values are possible;
+\item{r}{A pair of numeric values. The first gives the order of the
+generalized-mean price index used to calculate the
+elementary price indexes, defaulting to a geometric index. The second
+gives the order of the generalized-mean price index used to aggregate the
+elementary price indexes, defaulting to an arithmetic index. Other values
+are possible;
 see \code{\link[gpindex:generalized_mean]{gpindex::generalized_mean()}} for details.}
 
-\item{formula}{A two-sided formula with prices on the left-hand
-side. For \code{carry_forward()} and \code{carry_backward()}, the right-hand side
-should have time periods and products (in that order); for
-\code{shadow_price()}, the right-hand side should have time period, products,
-and elementary aggregates (in that order).}
+\item{method}{Name of the imputation method, one of \code{"overall-mean"},
+\code{"carry-forward"}, or \code{"carry-backward"}.}
+
+\item{formula}{A two-sided formula, or something that can be coerced into
+one, with prices on the left-hand
+side and time periods and products on the right-hand side (in that order).}
 }
 \value{
-A numeric vector of prices with missing values replaced (where possible).
+A numeric vector or matrix of prices with missing values replaced
+(where possible).
 }
 \description{
-Impute missing prices using the carry forward or shadow price method.
+Impute missing prices using the carry forward/backward or the
+self-correcting overall mean method.
 }
 \details{
 The carry forward method replaces a missing price for a product by the price
@@ -96,22 +106,27 @@ towards 1, and is usually avoided; see paragraph 6.61 in the CPI manual
 (2020). The carry backwards method does the opposite, but this is rarely
 used in practice.
 
-The shadow price method recursively imputes a missing price by the value of
+The self-correcting overall mean method recursively imputes a missing price
+by the value of
 the price for the same product in the previous period multiplied by the
 value of the period-over-period elementary index for the elementary aggregate
 to which that product belongs. This requires computing and aggregating an
 index (according to \code{pias}, unless \code{pias} is not supplied) for
-each \code{period}, and so these imputations can take a while. The index
+each \code{period}. The index
 values used to do the imputations are not returned because the index needs
-to be recalculated to get correct percent-change contributions.
-
-Shadow price imputation is referred to as self-correcting overall mean
-imputation in chapter 6 of the CPI manual (2020). It is identical to simply
+to be recalculated to get correct percent-change contributions. It is
+identical to simply
 excluding missing price relatives in the index calculation, except in the
 period that a missing product returns. For this reason care is needed when
 using this method. It is sensitive to the assumption that a product does not
 change over time, and in some cases it is safer to simply omit the missing
 price relatives instead of imputing the missing prices.
+
+Imputation works slightly differently depending on whether data are in a long
+or wide format. When \code{x} is a two-column of matrix of current and back prices
+(in that order), then imputation is done separately on the current price
+at a point in time and the back price at the next point in time. When \code{x} is
+a numeric vector then these two prices are necessarily the same.
 }
 \examples{
 prices <- data.frame(
@@ -121,9 +136,27 @@ prices <- data.frame(
   ea = rep(letters[1:2], 4)
 )
 
-carry_forward(prices, price ~ period + product)
+impute_prices(prices, price ~ period + product, method = "carry-forward")
+
+impute_prices(
+  prices,
+  price ~ period + product,
+  ea = ea,
+  method = "overall-mean"
+)
+
+# Can also be done with current price-back price formulation.
+prices$back_price <- with(
+  prices,
+  price[gpindex::back_period(period, product)]
+)
 
-shadow_price(prices, price ~ period + product + ea)
+impute_prices(
+  prices,
+  cbind(price, back_price) ~ period + product,
+  ea = ea,
+  method = "overall-mean"
+)
 
 }
 \references{
diff --git a/man/price_relative.Rd b/man/price_relative.Rd
index 76d2480..6a21b3e 100644
--- a/man/price_relative.Rd
+++ b/man/price_relative.Rd
@@ -26,8 +26,9 @@ ordering of time periods follows the levels of \code{period} to agree with
 \item{product}{A factor, or something that can be coerced into one, that
 gives the corresponding product identifier for each element in \code{x}.}
 
-\item{formula}{A two-sided formula with prices on the left-hand
-side, and time periods and products (in that order) on the
+\item{formula}{A two-sided formula, or something that can be coerced into
+one, with prices on the left-hand
+side and time periods and products (in that order) on the
 right-hand side.}
 }
 \value{
diff --git a/man/split_classification.Rd b/man/split_classification.Rd
index d0c26cb..4939bb9 100644
--- a/man/split_classification.Rd
+++ b/man/split_classification.Rd
@@ -20,7 +20,7 @@ See \code{\link[=strsplit]{strsplit()}}.}
 default separates levels by \code{"."}.}
 
 \item{pad}{A string used to pad the shorter labels for an unbalanced
-classification. The default pads with NA.}
+classification. The default pads with \code{NA}.}
 }
 \value{
 A list with a entry for each level in \code{x} giving the "digits" that
diff --git a/tests/Examples/piar-Ex.Rout.save b/tests/Examples/piar-Ex.Rout.save
index a4354e3..e60d655 100644
--- a/tests/Examples/piar-Ex.Rout.save
+++ b/tests/Examples/piar-Ex.Rout.save
@@ -710,10 +710,9 @@ levels 1 2 3
 > 
 > ### Name: impute_prices
 > ### Title: Impute missing prices
-> ### Aliases: impute_prices shadow_price shadow_price.default
-> ###   shadow_price.data.frame carry_forward carry_forward.default
-> ###   carry_forward.data.frame carry_backward carry_backward.default
-> ###   carry_backward.data.frame
+> ### Aliases: impute_prices impute_prices.default impute_prices.matrix
+> ###   impute_prices.numeric impute_prices.data.frame carry_forward
+> ###   carry_backward shadow_price
 > 
 > ### ** Examples
 > 
@@ -724,11 +723,38 @@ levels 1 2 3
 +   ea = rep(letters[1:2], 4)
 + )
 > 
-> carry_forward(prices, price ~ period + product)
+> impute_prices(prices, price ~ period + product, method = "carry-forward")
 [1] 1 2 3 4 5 6 7 4
-> 
-> shadow_price(prices, price ~ period + product + ea)
+>
+> impute_prices(
++   prices,
++   price ~ period + product,
++   ea = ea,
++   method = "overall-mean"
++ )
 [1]  1  2  3  4  5  6  7 12
+>
+> # Can also be done with current price-back price formulation.
+> prices$back_price <- with(
++   prices,
++   price[gpindex::back_period(period, product)]
++ )
+> 
+> impute_prices(
++   prices,
++   cbind(price, back_price) ~ period + product,
++   ea = ea,
++   method = "overall-mean"
++ )
+     price back_price
+[1,]     1          1
+[2,]     2          2
+[3,]     3          3
+[4,]     4          4
+[5,]     5          1
+[6,]     6          2
+[7,]     7          3
+[8,]    12          4
 > 
 > 
 > 
diff --git a/tests/test-making-price-indexes.R b/tests/test-making-price-indexes.R
index d2c679c..73c484d 100644
--- a/tests/test-making-price-indexes.R
+++ b/tests/test-making-price-indexes.R
@@ -217,7 +217,12 @@ chain(stack(ms_index1, ms_index2))
 ## -----------------------------------------------------------------------------
 ms_elementary2 <- ms_prices |>
   transform(
-    imputed_price = carry_forward(price, period = period, product = product)
+    imputed_price = impute_prices(
+      price,
+      period = period,
+      product = product,
+      method = "carry-forward"
+    )
   ) |>
   elementary_index(
     price_relative(imputed_price, period = period, product = product) ~
diff --git a/tests/test-making-price-indexes.Rout.save b/tests/test-making-price-indexes.Rout.save
index 8c93279..5206616 100644
--- a/tests/test-making-price-indexes.Rout.save
+++ b/tests/test-making-price-indexes.Rout.save
@@ -457,7 +457,12 @@ levels 202001    202002    202003    202004
 > ## -----------------------------------------------------------------------------
 > ms_elementary2 <- ms_prices |>
 +   transform(
-+     imputed_price = carry_forward(price, period = period, product = product)
++     imputed_price = impute_prices(
++       price,
++       period = period,
++       product = product,
++       method = "carry-forward"
++     )
 +   ) |>
 +   elementary_index(
 +     price_relative(imputed_price, period = period, product = product) ~
diff --git a/tests/testthat/test-extract-index.R b/tests/testthat/test-extract-index.R
index a6698ac..a0a3c97 100644
--- a/tests/testthat/test-extract-index.R
+++ b/tests/testthat/test-extract-index.R
@@ -39,7 +39,7 @@ test_that("head and tail work", {
 
 test_that("subscripting methods work", {
   expect_equal(epr[], epr)
-  expect_equal(index[], index)
+  expect_equal(index[1:5], index)
   expect_equal(
     epr[c(TRUE, FALSE, TRUE, TRUE), 2:1],
     elementary_index(
diff --git a/tests/testthat/test-impute-prices.R b/tests/testthat/test-impute-prices.R
index 9cc28ea..9a92ee7 100644
--- a/tests/testthat/test-impute-prices.R
+++ b/tests/testthat/test-impute-prices.R
@@ -6,42 +6,69 @@ pias <- with(
   )
 )
 
-sp <- shadow_price(ms_prices, price ~ period + product + business, pias = pias)
+sp <- impute_prices(
+  ms_prices,
+  price ~ period + product,
+  ea = business,
+  pias = pias,
+  method = "overall-mean"
+)
 
 test_that("a length 0 inputs returns a length 0 output", {
   expect_length(
-    carry_forward(
+    impute_prices(
       integer(0),
       period = factor(integer(0), 1:5),
-      product = integer(0)
+      product = integer(0),
+      method = "carry-forward"
     ),
     0
   )
   expect_length(
-    shadow_price(
+    impute_prices(
       integer(0),
       period = integer(0),
       product = integer(0),
       ea = integer(0),
-      pias = pias
+      pias = pias,
+      method = "overall-mean"
     ),
     0
   )
+  expect_equal(
+    impute_prices(
+      cbind(integer(0), integer(0)),
+      period = integer(0),
+      product = integer(0),
+      ea = integer(0),
+      pias = pias,
+      method = "overall-mean"
+    ),
+    matrix(numeric(0), 0, 2)
+  )
 })
 
 test_that("imputing shadow prices does noting", {
   expect_equal(
     sp,
-    shadow_price(ms_prices, sp ~ period + product + business, pias = pias)
+    impute_prices(
+      ms_prices,
+      sp ~ period + product,
+      ea = business,
+      pias = pias,
+      method = "overall-mean"
+    )
   )
 
   expect_equal(
     sp,
-    shadow_price(
+    impute_prices(
       ms_prices,
-      sp ~ period + product + business,
+      sp ~ period + product,
+      ea = business,
       pias = pias,
-      weights = price
+      weights = price,
+      method = "overall-mean"
     )
   )
 })
@@ -81,7 +108,13 @@ test_that("imputing with an improper pias does nothing", {
 
   expect_equal(
     ms_prices$price,
-    shadow_price(ms_prices, price ~ period + product + business, pias = pias2)
+    impute_prices(
+      ms_prices,
+      price ~ period + product,
+      ea = business,
+      pias = pias2,
+      method = "overall-mean"
+    )
   )
 })
 
@@ -130,7 +163,13 @@ test_that("jumbling prices does nothing", {
   )
   ms_prices <- ms_prices[jumble, ]
   expect_equal(
-    shadow_price(ms_prices, price ~ period + product + business, pias = pias),
+    impute_prices(
+      ms_prices,
+      price ~ period + product,
+      ea = business,
+      pias = pias,
+      method = "overall-mean"
+    ),
     sp[jumble]
   )
 })
@@ -141,6 +180,36 @@ test_that("carrying forward/backwards imputation works", {
     period = gl(5, 1),
     product = gl(1, 5)
   )
-  expect_equal(carry_forward(df, price ~ period + product), c(NA, 1, 2, 2, 3))
-  expect_equal(carry_backward(df, price ~ period + product), c(1, 1, 2, 3, 3))
+  expect_equal(
+    impute_prices(df, price ~ period + product, method = "carry-forward"),
+    c(NA, 1, 2, 2, 3)
+  )
+  expect_equal(
+    impute_prices(df, price ~ period + product, method = "carry-backward"),
+    c(1, 1, 2, 3, 3)
+  )
+})
+
+test_that("imputing with a matrix works", {
+  ms_prices$back_price <- with(
+    ms_prices,
+    price[gpindex::back_period(period, product)]
+  )
+  sp2 <- impute_prices(
+    ms_prices,
+    cbind(price, back_price) ~ period + product,
+    ea = business,
+    pias = pias,
+    method = "overall-mean"
+  )
+  expect_equal(sp2[, 1L], sp)
+  expect_equal(
+    impute_prices(
+      matrix(c(1, NA, NA, 4, 5, 6, 3, 2, 1, 1, NA, 3), ncol = 2),
+      rep(1:2, each = 3),
+      rep(1:3, 2),
+      method = "carry-forward"
+    ),
+    matrix(c(1, 2, 1, 4, 5, 6, 3, 2, 1, 1, 2, 3), ncol = 2)
+  )
 })
diff --git a/tests/testthat/test-known-values.R b/tests/testthat/test-known-values.R
index b1c200d..4d54f1b 100644
--- a/tests/testthat/test-known-values.R
+++ b/tests/testthat/test-known-values.R
@@ -45,10 +45,12 @@ test_that("matched-sample index works", {
     )
   )
 
-  sp <- shadow_price(
+  sp <- impute_prices(
     ms_prices,
-    price ~ period + product + business,
-    pias = pias
+    price ~ period + product,
+    ea = business,
+    pias = pias,
+    method = "overall-mean"
   )
 
   rel <- price_relative(ms_prices, sp ~ period + product)