2026-11-05
\[ \newcommand{\vy}{\mathbf{y}} \newcommand{\vx}{\mathbf{x}} \newcommand{\vb}{\mathbf{b}} \newcommand{\ve}{\mathbf{e}} \newcommand{\vone}{\mathbf{1}} \newcommand{\vbeta}{\boldsymbol{\beta}} \newcommand{\vmu}{\boldsymbol{\mu}} \newcommand{\mX}{\mathbf{X}} \newcommand{\mZ}{\mathbf{Z}} \newcommand{\mV}{\mathbf{V}} \newcommand{\mI}{\mathbf{I}} \newcommand{\mJ}{\mathbf{J}} \newcommand{\mSigma}{\boldsymbol{\Sigma}} \newcommand{\vg}{\mathbf{g}} \newcommand{\vc}{\mathbf{c}} \newcommand{\vz}{\mathbf{z}} \newcommand{\vlambda}{\boldsymbol{\lambda}} \newcommand{\mA}{\mathbf{A}} \newcommand{\mB}{\mathbf{B}} \newcommand{\mK}{\mathbf{K}} \newcommand{\mD}{\mathbf{D}} \newcommand{\mL}{\mathbf{L}} \newcommand{\mQ}{\mathbf{Q}} \newcommand{\E}{\mathrm{E}} \newcommand{\Var}{\mathrm{Var}} \newcommand{\Cov}{\mathrm{Cov}} \newcommand{\tr}{\operatorname{tr}} \]
This vignette sets out what the three filters compute: the HP filter as a penalised least-squares problem, the criterion of the modified HP filter and the exact cross-validation it approximates, the weights, gains and end-point loss that measure the end-point bias, and the end-point schemes of the fully modified HP filter. It also says how the package computes them, in \(O(T)\) operations per smoothing parameter.
With the series \(\vx = (x_1, \dots, x_T)^\top\) and the \((T - 2) \times T\) matrix \(\mK\) of second differences (rows \((1, -2, 1)\)), the trend minimises \(\lVert\vx - \vg\rVert^2 + \lambda\lVert\mK\vg\rVert^2\), so \[ \hat{\vg} = (\mI + \lambda\mA)^{-1}\vx = \mB\vx, \qquad \mA = \mK^\top\mK , \] and \(\hat{\vc} = \vx - \hat{\vg}\). Because \(\mK\vone = \mathbf{0}\), every row of \(\mB\) sums to one: the trend is a weighted average of the observations.
Computation. \(\mI + \lambda\mA\) is symmetric, positive definite and pentadiagonal. The package factors it as \(\mL\mD\mL^\top\) with a unit lower-triangular \(\mL\) of bandwidth two and solves in \(O(T)\) operations; version 0.1.0 solved a dense \(T \times T\) system, \(O(T^3)\), for every \(\lambda\) of the search.
McDermott (1997) chooses \(\lambda\) by cross-validation: leave out one observation, predict it from the trend of the others, and sum the squared errors. Generalised cross-validation (Craven and Wahba 1979) replaces the leave-one-out errors by \[
\mathrm{GCV}(\lambda) = \frac{T^{-1}\sum_t\left(x_t - \hat{g}_t(\lambda)\right)^2}{\left(1 - T^{-1}\tr\mB(\lambda)\right)^2}
\] (Choudhary et al. 2014, eq. 5), and the papers and their code use the approximation \[
\mathrm{GCV}(\lambda) \approx \frac{1}{T}\left(1 + \frac{2T}{\lambda}\right)\sum_t\left(x_t - \hat{g}_t(\lambda)\right)^2
\] (their eq. 6, attributed to Silverman 1984). mhp_filter() minimises the approximation over \(\lambda = 1, \dots,\) max_lambda by default, and the exact criterion with criterion = "exact": the eigenvalues \(\mu_i\) of \(\mA\), computed once, give \(\tr\mB(\lambda) = \sum_i (1 + \lambda\mu_i)^{-1}\) and the residual sum of squares for every \(\lambda\) in \(O(T)\) operations.
data(us_macro)
x <- log(us_macro$gdp)
ap <-
mhp_filter(
x = x
, max_lambda = 20000
, as_dt = FALSE
, log_scale = FALSE
, criterion = "approximate"
)
ex <- suppressWarnings(
mhp_filter(
x = x
, max_lambda = 20000
, as_dt = FALSE
, log_scale = FALSE
, criterion = "exact"
)
)
c(approximate = ap$lambda, exact = ex$lambda)approximate exact
683 1 d <- data.table::data.table(
lambda = rep(seq_len(length.out = 20000), 2)
, criterion = c(ap$gcv_path / min(ap$gcv_path), ex$gcv_path / min(ex$gcv_path))
, type = rep(c("approximate (eq. 6)", "exact GCV (eq. 5)"), each = 20000)
)
ggplot2::ggplot(data = d, mapping = ggplot2::aes(x = lambda, y = criterion, colour = type)) +
ggplot2::geom_line() +
ggplot2::scale_x_log10(breaks = scales::breaks_log(n = 6), labels = scales::label_comma()) +
ggplot2::labs(x = "lambda (log scale)", y = "criterion / its minimum", colour = NULL) +
mhp_theme(
base_size = 11
, base_family = "serif"
, ticks_inward = TRUE
, legend_position = "bottom"
)
The exact criterion falls steadily as \(\lambda\) falls, down to \(\lambda = 1\), where the trend almost interpolates the series: generalised cross-validation is known to undersmooth series whose deviations from the trend are autocorrelated, as business cycles are. The approximation adds a much heavier penalty for small \(\lambda\) (at \(\lambda = 683\) the factor \(1 + 2T/\lambda\) is 1.59, against 1.17 for the exact one), and that penalty, not cross-validation as such, gives the modified HP filter its interior minimum. The package keeps the published approximation as the default.
Row \(t\) of \(\mB\) holds the weights of \(\hat{g}_t\). In the middle of the series they are symmetric; near the ends they are not, and the end observation gets a large weight in its own trend (Hanif et al. 2017, Figure 1).
w <-
hp_weights(
x = 50
, lambda = 100
, form = "rows"
, rows = NULL
)
plot(x = w)
The gain of \(\hat{g}_t\) at frequency \(\omega\), \[ P_t(\omega) = \left|\sum_{j=1}^{T} b_{tj}\,e^{i\omega(j - t)}\right| , \] is the factor by which the trend scales an oscillation of that frequency (eq. 6 of Hanif et al. 2017). In the middle the gain falls to zero at high frequencies; at the ends it does not, so the trend keeps part of the cycle.
g <-
hp_gain(
x = 50
, lambda = 100
, form = "rows"
, theta = seq(from = 0, to = pi, length.out = 101)
, rows = NULL
)
plot(x = g)
The package computes \(\sum_j b_{tj}e^{i\omega j}\) for all \(t\) at once as \(\mB\cos(\omega\cdot)\) and \(\mB\sin(\omega\cdot)\): two solves per frequency, without forming \(\mB\).
Bloechl (2014) measures the distortion by the distance of each gain from the gain at the centre \(c\) of the series, \[ l(t) = \theta\sum_{\omega}\left[P_c(\omega) - P_t(\omega)\right]^2, \qquad L = \sum_{t=1}^{T} l(t), \] with \(\omega = 0, 0.1, \dots, 3.1\) and \(\theta = 0.1\) (eqs. 7 and 8 of Hanif et al. 2017). It depends on \(T\) and the smoothing parameters only.
hl <-
hp_loss(
x = 50
, lambda = 100
, form = "rows"
, theta = seq(from = 0, to = 3.1, by = 0.1)
, weight = 0.1
)
plot(x = hl)
Hanif et al. (2017) raise the smoothing parameter of the \(k\) values at each end: \(\lambda_t = \lambda\) in the middle and \(\lambda + \alpha j^{i}\) for the \(j\)-th value from the middle, linear (\(i = 1\)) or quadratic (\(i = 2\)).
hp_schedule(
T = 50
, lambda = 100
, alpha = 10
, k = 20
, power = 1
, form = "rows"
) [1] 300 290 280 270 260 250 240 230 220 210 200 190 180 170 160 150 140 130 120
[20] 110 100 100 100 100 100 100 100 100 100 100 110 120 130 140 150 160 170 180
[39] 190 200 210 220 230 240 250 260 270 280 290 300The published MATLAB code multiplies row \(t\) of \(\mA\) by \(\lambda_t\): \[ \hat{\vg} = (\mI + \mathrm{diag}(\lambda_t)\mA)^{-1}\vx \qquad \text{(\texttt{form = "rows"}, the default)}. \] This matrix is not symmetric, but \(\mI + \mathrm{diag}(\lambda_t)\mA = \mathrm{diag}(\lambda_t) (\mathrm{diag}(1/\lambda_t) + \mA)\), so the package solves the symmetric pentadiagonal system \((\mathrm{diag}(1/\lambda_t) + \mA)\hat{\vg} = \mathrm{diag}(1/\lambda_t)\vx\). The penalised least-squares problem with one smoothing parameter per second difference, \(\lVert\vx - \vg\rVert^2 + \sum_r\lambda_r(\mK\vg)_r^2\), gives instead \[ \hat{\vg} = (\mI + \mK^\top\mathrm{diag}(\lambda_r)\mK)^{-1}\vx \qquad \text{(\texttt{form = "differences"})}, \] the form in which Bloechl (2014) and eq. 9 of the paper write the scheme. Both are \(O(T)\); for a constant \(\lambda\) they coincide.
y <- us_macro$gdp
rows <-
fmhp_filter(
x = y
, max_lambda = 3000
, k_linear = 15
, k_quadratic = NULL
, alpha_linear = seq(from = 0, to = 2000, by = 10)
, alpha_quadratic = seq(from = 0, to = 2000, by = 5)
, schemes = c("linear", "quadratic")
, form = "rows"
, lambda = NULL
, log_scale = TRUE
, as_dt = FALSE
)
diffs <-
fmhp_filter(
x = y
, max_lambda = 3000
, k_linear = 15
, k_quadratic = NULL
, alpha_linear = seq(from = 0, to = 2000, by = 10)
, alpha_quadratic = seq(from = 0, to = 2000, by = 5)
, schemes = c("linear", "quadratic")
, form = "differences"
, lambda = NULL
, log_scale = TRUE
, as_dt = FALSE
)
rbind(
rows = rows$schemes[scheme == rows$scheme]
, differences = diffs$schemes[scheme == diffs$scheme]
) scheme k alpha loss lambda gcv
<char> <int> <num> <num> <int> <num>
1: quadratic 22 15 0.6163442 1448 0.0003335623
2: quadratic 22 95 1.1645828 1410 0.0003353531The defaults are those of the MATLAB code: \(\lambda \le 3000\); \(k = 15\) and \(\alpha \in \{0, 10, \dots, 2000\}\) for the linear scheme; \(k = \mathrm{round}(\min(22, T/2 - 5))\) and \(\alpha \in \{0, 5, \dots, 2000\}\) for the quadratic one; logarithms. The text of the paper states \(k = 20\); k_linear and k_quadratic set it. MATLAB rounds halves up, and so does the package (the centre \(c\) of the loss and the quadratic \(k\)); vignette 4 shows that this reproduces the published losses.
Bloechl, A. (2014). Reducing the excess variability of the Hodrick-Prescott filter by flexible penalization. Discussion Paper 2014-1, University of Munich.
Choudhary, M. A., Hanif, M. N. and Iqbal, J. (2014). On smoothing macroeconomic time series using the modified HP filter. Applied Economics, 46(19), 2205–2214.
Craven, P. and Wahba, G. (1979). Smoothing noisy data with spline functions. Numerische Mathematik, 31, 377–403.
Hanif, M. N., Iqbal, J. and Choudhary, M. A. (2017). Fully modified HP filter. SBP Working Paper Series, No. 88, State Bank of Pakistan.
McDermott, C. J. (1997). An automatic method for choosing the smoothing parameter in the HP filter. Unpublished, International Monetary Fund.
Silverman, B. W. (1984). A fast and efficient cross-validation method for smoothing parameter choice in spline regression. Journal of the American Statistical Association, 79(387), 584–589.