# The HMTS method (Ishaak, Ouwehand & Remøy 2024) — step-by-step description for benchmark G

Source (primary). Ishaak, F., Ouwehand, P. & Remøy, H. (2024). "Constructing Limited-Revisable and Stable CPPIs for Small Domains." Journal of Official Statistics 40(3), 380–408. DOI 10.1177/0282423X241246617. CC BY-NC. Full text obtained on 2026-10-21 from the TU Delft repository (final published version, https://repository.tudelft.nl/file/File_f298b447-720b-44eb-84ab-8437ee06d026); SAGE blocks automated access. Page numbers below are the journal page numbers printed in that PDF. Received 9 March 2023, accepted 18 December 2023 (p. 407).

Source (implementation). CRAN package REPS 1.1.1 (Statistics Netherlands; authors include Ishaak and Ouwehand; maintainer V. Gajadhar), files R/calculate_hmts.R, R/hmts_helpers.R, R/ts_helpers.R, obtained from https://cran.r-project.org/src/contrib/REPS_1.1.1.tar.gz on 2026-10-21. The package cites the JOS 2024 article as the source of HMTS and (in roxygen comments) an "Ishaak, Ouwehand, Remoy & De Haan (2023)" version that could not be located bibliographically.

Everything below is paraphrased or quoted from these two sources. Where a detail is not stated, this is said explicitly.

# 0. Purpose and context (pp. 380–386)

Target: a constant-quality commercial property price index (CPPI) for small domains with heterogeneous objects, for official statistics. The authors state that hedonic time-dummy pooling fixes small samples but "involves backward revision", and that the rolling time dummy method (RTD) "partly resolves this issue by pooling the data in a specified time window" but applied to commercial property "still often results in a volatile index" (p. 381).
Three practical properties sought (§2.2, pp. 384–386): (1) a stable index ("not to remove market volatility, but only to remove volatility that is intrinsic to the index method or non-existing volatility", p. 384); (2) limited revision — they note EU member states "are by legislation bound to only one preliminary period" for the HPI (p. 385); (3) early detection of turning points, defined as a "structural" change in price development from positive to negative or vice versa (p. 385).
Three methodological properties examined (§2.1, pp. 383–384): identity test, time-reversal test, circularity test.
Acronym: HMTS = "Hedonic Multilateral Time series Splice" (p. 395). The package expands it as "Hedonic Multilateral Time series re-estimation Splicing".

# 1. Overview: a four-step procedure (§3, p. 386)

"(1) hedonic imputations, (2) multilateral imputations, (3) time series re-estimated imputations, and (4) window splicing." The authors write "This is the first study that combines above techniques to construct a price index." Figure 1 (p. 386) is a schematic of the four steps.

# 2. Step 1 — Hedonic imputations (§3.1, pp. 386–388)

(a) Observation/hedonic equation. A log-linear hedonic regression estimated separately for each period t by OLS (Equation (1), p. 387, "based on Eurostat (2013b, 50)"):

ln p_{t,i} = α_t + Σ_{k=1}^{K} β_t^k x_{t,i}^k + ε_{t,i},   ε_{t,i} ~ N(0, σ²)

with p_{t,i} the price of property i sold in period t, α_t the period-t intercept, β_t^k the period-t coefficient of characteristic k, x_{t,i}^k the characteristic, ε_{t,i} the error.

The paper then recalls Laspeyres and Paasche hedonic double-imputation (HDI) indexes, Equations (2)–(3) (p. 387), and defines the average imputed price of period t with base period b, Equation (4) (p. 389):

p̂_{t(b)} = exp(α̂_t) · exp( Σ_{k=1}^{K} β̂_t^k · x̄_b^k )

i.e. the period-t coefficients applied to the mean characteristics of base period b (x̄_b^k = Σ_i x_{b,i}^k / N_b, p. 388). In words: "the mix of sold real estate in period [b] is tracked through time. Its prices are imputed in the next periods as if they were sold in another period" (p. 388). Only the imputed prices are used; "Equations (2) (Laspeyres) and (3) (Paasche) are not performed to construct the HMTS index" (p. 388).

Implementation detail (REPS calculate_hedonic_imputationmatrix, hmts_helpers.R): the model is log(price) ~ numerical_variables + categorical_variables, fitted with stats::lm.fit on the period-t rows (X_dyn, y_dyn); coefficients that are NA are set to 0; the imputation is X_base %*% fitcoef on the individual rows of the base period b and the series value is exp(mean(preds)) — i.e. the geometric mean of the individual imputed prices of base-period properties, which equals the exp of the fitted value at the base-period means for a linear-in-parameters model. Rows with missing values in any model variable are dropped (na.omit). Factor levels are fixed on the whole dataset so that the design matrix is comparable across periods. Numerical variables are not log-transformed by the package; the vignette recommends the user log floor area beforehand. In the paper's application the dependent variable and the numerical regressors (floor area, distance to train station) are in natural logs (Appendix 1 note, p. 407).

What the paper does not say: no minimum number of observations per period is given for Step 1; no treatment of periods where OLS is rank-deficient is specified beyond the package's NA → 0 rule. The authors state as a limitation that "in cases where there are (almost) no observations, the HMTS index in its current form does not work. Each period should have a few observations with explanatory variables that show variation" and suggest, untested, a "multiperiod regression" fallback (§5.1, p. 404).

# 3. Step 2 — Multilateral imputations (§3.2, pp. 388–390)

Every reporting period t is imputed with every possible base period b, filling the full T×T matrix of p̂_{t(b)} (Table 3, p. 390: rows = base period, columns = reporting period). "This results in multiple series and all of these cover the complete time span, hence the term 'multilateral'" (p. 389). The idea is "inspired by the multilateral GEKS ... method as described by Willenborg (2017, 2018)" (p. 389).

A single period-t imputation would be the geometric mean over base periods, Equation (5), p. 389:

P̄_t = ( Π_{b=0}^{T} p̂_{t(b)} )^{1/T},   t = 0, 1, …, T

and a circular index Equation (6): I_{0,t} = P̄_t / P̄_0. But "Equations (5) and (6) are not performed at this stage for the final index. The imputations resulting from Equation (4) are used in the next step" (p. 389).

Implementation detail: in REPS the outer loop is over base periods (current_period), the inner loop over reporting periods; for a base period b the series p̂_{·(b)} spans reporting periods 1…T (or, in production mode, only up to the periods allowed by the update scheme — see §5 below).

# 4. Step 3 — Time-series re-estimation of each base-period series (§3.3, pp. 390–392)

(b) The state-space model. Each base-period series p̂_{t(b)}, t = 0…T, is treated as a noisy observed series and smoothed with a univariate structural time-series model. The paper first presents the random-walk-plus-noise (local level) model, Equations (7)–(8), p. 391 (citing Commandeur & Koopman 2007, p. 9), then states (p. 391): "many models were tested and the best model was chosen by looking at the most suitable fit. The model showed, among which, the lowest AIC (presented in Appendix 2). The model of choice is a local linear trend model. ... Furthermore, the model sets a hyperparameter of the level component equal to zero, while the slope component remains stochastic. This model is also known as a smooth trend model, and since the slope disturbances affect the level component indirectly, generally yield a smoother estimate of the level component."

Equations (9)–(11), p. 392 (per base period b):

p̂_{t(b)}   = μ_{t(b)} + ε_{t(b)},        ε_{t(b)} ~ NID(0, σ_ε²)      (9)   observation
μ_{t+1(b)}  = μ_{t(b)} + ν_{t(b)}                                     (10)  level (no disturbance)
ν_{t+1(b)}  = ν_{t(b)} + ζ_{t(b)},        ζ_{t(b)} ~ NID(0, σ_ζ²)      (11)  slope

with μ_{t(b)} the unobserved level and ν_{t(b)} the slope. The paper writes "μ_{t(b)} = the unobserved level in period t for the series with base period b; ν_{t(b)} = the slope component in period t (deterministic)". So there are two hyperparameters per series: σ_ε² (observation noise) and σ_ζ² (slope noise); the level disturbance variance is fixed at 0.

Hyperparameters: "The model can be estimated by applying a two-step procedure. First, the hyperparameters, that is, the variances of the error terms, are estimated via a Maximum Likelihood procedure. Then the Kalman filter (Kalman 1960) is applied to estimate the state variables (level and slope)" (p. 392). Filtering is forward, then "a backward pass is applied ... This is called smoothing. In our case we are interested in these smoothed estimates" (p. 392). The re-estimated imputation is the smoothed level, Equation (12), p. 392:

p̃_{t(b)} = p̂_{t(b)} − ε̂_{t(b)}     ( = smoothed μ̂_{t(b)} )

Pooling across domains / across base-period series: none. Each base-period series of each domain is re-estimated independently with its own ML hyperparameters (one model fit per base period; Appendix 2, p. 408, reports log-likelihood and AIC per base period). There is no cross-domain or cross-series hierarchy; the "borrowing" in HMTS is purely temporal (adjacent periods of the same series) plus the averaging across base periods in Step 4.

Implementation detail (REPS ts_helpers.R, exact):

The series is log-transformed before modelling (original_series <- log(original_series)), so the state-space model is fitted to ln p̂_{t(b)}, and the trend is exponentiated back (exp(signalsubconf[,1])). (The paper's Equations (9)–(12) are written in price levels; the package works in logs.)
Model: KFAS::SSModel(series ~ SSMtrend(2, Q = list(0, matrix(NA))), H = matrix(NA)) — a degree-2 trend (local linear trend) with level variance fixed at 0 and slope variance NA (estimated), observation variance H estimated. This matches the paper's smooth-trend model.
Starting values: all hyperparameters start at log(0.1); a one-pass Kalman run computes a "scale factor" (1/(n−d)) Σ v_t' F_t^{-1} v_t (d = number of diffuse initial states) and the log-start values are shifted by 0.5·log(scale) (citing Commandeur & Koopman eq. 11.7 / SsfPack manual 9.4).
Estimation: KFAS::fitSSM(..., method = "BFGS", hessian = TRUE) (maximum likelihood); then KFAS::KFS(filtering = "state", smoothing = "state"); the returned trend is the smoothed level signal (KFAS::signal(states = c("level","regression"))), with 95% bands ± 1.96·sqrt(var) computed but not used for the index.
Diagnostics stored per base period: log-likelihood, log-likelihood/n, AIC, AIC/n, number of diffuse states q, number of estimated variances w, σ and σ² for "meas" and "slope" with t-values, and first/last period level and slope estimates with standard errors. (Appendix 2 of the paper, p. 408, prints exactly these LLM and AIC values per base period for 2008Q1–Q4 and 2022Q1–Q4; e.g. 2022Q4: LL 84.3, AIC −160.6, n = 60.)
The roxygen comments say the time-series step "can be turned off with a parameter" (and mention state_space_model = NULL), but in REPS 1.1.1 as shipped calculate_hedonic_imputationmatrix() calls calculate_trend_line_kfas() unconditionally and calculate_hmts() exposes no such switch; the imputation argument of calculate_hedonic_index() applies to Laspeyres/Paasche only. The HMTS-relevant arguments of calculate_hedonic_index(method = "hmts") are number_preliminary_periods (package default 3 if omitted), production_since (default NULL) and resting_points; window_length is the rolling-time-dummy argument (package default 5), not an HMTS argument.

What is not stated: no seasonal component is used (start-value slots for "seas" and "scaling" exist in the code but the selected model has none); no diffuse-initialisation details beyond KFAS defaults; no statement on how many periods are needed for ML estimation to be reliable (the empirical application starts the state-space re-estimation "as if 2015 Q1 was the first moment of index construction" with 28 quarters available, Appendix 2 note, p. 408, and Table 5 footnote says the scheme "is implemented in practice with around twenty starting periods", p. 394).

# 5. Step 4 — Window splice (§3.4, pp. 392–395)

(c) Why: both the multilateral step (adding a reporting period adds a base-period series, changing the geometric mean for every period) and the smoothing step (smoothed values use preceding and succeeding periods) revise the whole series each period (p. 393). "The main idea of window splicing is to keep a limited number of periods provisional instead of the entire time series" (p. 393).

Two-part rule (pp. 393–395):

Update scheme for the re-estimated series (Table 5, p. 394). Fix a window of w periods. When period t is first computed it is provisional; "In this period and the next two periods [for w = 3], the figure for period t is subject to revision. In period t+1, the series with base period t is extended with reporting period t+1. In period t+2, the series with base period t is extended again with reporting period t+2. In both calculations, the figure for period t is revised. This will be the last time, since the series with base period t is fixed from that point on" (p. 393). So each base-period series is re-estimated only until it is w periods long, then frozen; the number of provisional periods equals w − 1.
Index transformation, Equation (17), p. 394: within each base-period series, I_{0,t(b)} = p̃_{t(b)} / p̃_{0(b)}, so that series with different base periods become comparable.
Window splice, Equation (18), p. 394: the final index for reporting period t is the geometric average of the w index values with base periods t−(w−1), …, t:
```
I^{HMTS,w}_{0,t} = ( Π_{b = t−(w−1)}^{t} I_{0,t(b)} )^{1/w},   t = w−1, w, …, T
```
e.g. for w = 3, I_{0,t} is the geometric mean of I_{0,t(t−2)}, I_{0,t(t−1)}, I_{0,t(t)} (p. 394). For the first periods (t ≤ w−2) the window is truncated at the start, Equation (19), p. 395: I^{w}{0,t} = ( Π{b=0}^{w−1} I_{0,t(b)} )^{1/w}. Table 6 (p. 395) shows the three-period example.

Limited-revisable — how many periods revise: with window length w, a period is revised in the w−1 publications following its first release; "The HMTS index with three preliminary periods (four period window splice)" is the configuration used in the empirical section (§4.2, p. 396), i.e. w = 4, three preliminary periods, every period published four times (first estimate, then after +1, +2, +3 periods; §4.3, p. 397). In the Dutch office-market simulation around the 2015 turning point, first estimates were on average 0.82 index points from the final, second estimates 0.29, third 0.03; "two or three preliminary periods seem to be appropriate" (p. 397). Note: the paper's "three preliminary periods" is the package parameter number_preliminary_periods (default 3 in the package doc; the vignette example uses 1 with window_length = 4).

Choice of w: not fixed by the paper. "The choice of window length is rather important ... greater transitivity goes along with a longer window. On the other hand, a shorter window minimizes the number of revisions" (p. 395); they refer to Hill et al. (2022) for a window-length criterion and announce "Another way to look for an optimal window length—and specifically applied to the HMTS—is presented in section 4.2" (p. 395) — in the published §4.2 the only window-related content is the statement above about revisions by estimate number; no formal selection criterion for w is given.

Implementation detail (REPS calculate_hmts_index): number_preliminary_periods controls the update scheme; when production_since is NULL all periods are treated as final (no splice — the whole matrix is averaged: window <- imputations[, -1]). With production_since set, the package rebuilds the series period by period ("To simulate a series, where 1 period a time expires (as in production)"), keeping for each reporting period the index values of the base-period series that are still within the window, then takes calculate_geometric_average(na.omit(...)) per period and rescales to reference_period = 100. The window splice is "a specified (parameter) window ... of index figures that continues in the calculation" and "can be turned off".

# 6. (d) What "small domains" means operationally

The paper gives no numerical threshold defining a small domain. The empirical domain is Dutch office-building transactions 2008–2022, quarterly, "approximately 150 to 1,100 per quarter" (§4.1, p. 396); a robustness run uses "a randomly selected subset of the data with fifty observations per quarter" (§4.3, p. 398), which still produced a usable, "somewhat less stable" index. Appendix 1 (p. 407) shows n = 461 (2008Q1) and n = 785 (2022Q4) in the period regressions.
The only stated lower bound is qualitative (§5.1, p. 404): each period needs "a few observations with explanatory variables that show variation", otherwise the method "does not work".
Hedonic variables in the application (Appendix 1, p. 407): log floor area, log distance to train station, neighbourhood segment dummies (B–E), "In Amsterdam?" dummy, building-age class dummies; adjusted R² 0.85 (2008Q1) and 0.76 (2022Q4). The authors report that "the HMTS method performs well even with the most basic hedonic variables: floor area, location, and building age" (p. 396).

# 7. Properties and evaluation reported (§4.3–4.5, pp. 397–403)

Bootstrap 95% intervals (500 resamples with replacement) show "improved reliability compared to a Fisher HDI" (p. 398). Leave-one-out influence: individual-observation shares range −0.4 to +0.2 index points for HMTS vs −0.67 to +0.27 for Fisher HDI; mean absolute share 0.03 vs 0.06 (pp. 398–399).
Identity test: not satisfied (the index returns to 100 only slowly after prices are reset; "The risk is that [it] responds with a lag to extreme price changes", p. 400). Time reversal: satisfied to seven decimals without the splice, violated with it (p. 401). Circularity: lost after Steps 3–4, but indices with different start years "differ the most at the start of the series" and are otherwise similar (pp. 401–402).
Versus RTD (window 3, movement splice, as in de Haan 2015): the RTD "still shows volatility similar to the Fisher HDI index" and "shows drift, caused by 'weak' periods in 2013"; longer RTD windows move it toward Fisher (pp. 402–403).
Lead/lag: pairwise Granger tests between HMTS and (i) HDI with first base, (ii) multilateral HDI, (iii) Fisher HDI — "there is no clear lag of the HMTS. If any, the HMTS leads the other indices" (p. 403).
Turning-point detection in the 2015 simulation: the first two estimates of the turning-point quarter pointed down, the third was corrected up — "the desired accuracy lags two periods" (p. 397).

# 8. Checklist for implementing benchmark G faithfully (derived strictly from the above)

Per period t: OLS of ln price on the hedonic variables (period-specific coefficients).
For every base period b and reporting period t: p̂_{t(b)} = exp(mean of X_b β̂_t) (log scale mean of base-period rows under period-t coefficients).
For every b: fit a smooth-trend model (level variance = 0, slope variance and observation variance by ML; KFAS SSMtrend(2, Q=list(0,NA)), H=NA, BFGS) to ln p̂_{·(b)}, take the smoothed level, exponentiate → p̃_{t(b)}.
Index each series on its own first period: I_{0,t(b)} = p̃_{t(b)} / p̃_{0(b)}.
Production/vintage mode: a base-period series b is extended only while it has < w periods; afterwards frozen. Final index for t = geometric mean of I_{0,t(b)} for b = t−(w−1)…t (truncated at the start). Paper configuration: w = 4 (three preliminary periods). Rescale to the reference period.
Report: revisions by estimate number (1st…wth), bootstrap CI (500 draws), leave-one-out influence, identity/time-reversal/circularity checks, Granger lead/lag vs Fisher HDI — as in §4.

Open points that the sources do not settle (must be chosen and documented by us): minimum observations per period for Step 1; how to handle a period with no valid regression (paper suggests but does not test a multi-period regression); the window length w (paper: 4 in the application, no criterion); whether to run the state-space step on logs (package) or levels (paper equations) — the package behaviour (logs) is what Statistics Netherlands ships; the number of "starting periods" before the splice regime begins (paper: "around twenty", p. 394; 28 quarters in the application).