Introduction

The jumps package implements the method described in Maranzano and Pelagatti (2025), “A Hodrick-Prescott filter with automatically selected breaks”, Economic Modelling 150, 107132, doi:10.1016/j.econmod.2025.107132. The terms jumps and breaks are used interchangeably in the package documentation. Indeed, in the initial version of the article, we used the term jumps, and one of the referees suggested that this term is generally devoted to the jump processes used in mathematical finance, while breaks is generally used in econometrics. We agreed with the referee and decided to use the term breaks in the article’s final version. However, we decided to keep the term jumps in the package name and in the function names because we had already written most of it.

We introduced an innovation to the well-known Hodrick-Prescott filter (HPF) that allows for a small number of discontinuities in the otherwise smooth filter. The number, positions and magnitudes of these discontinuities, which we name breaks or jumps, are automatically estimated from the data. The method is based on the minimization of the sum of squared residuals of the HPF subject to a penalty on the number of breaks. The penalty is chosen by the user and can be set to zero, in which case the method reduces to the standard HPF. The technique is implemented in the hpj function, which is the main function of the package. The function hpj is a wrapper for other functions that implement various variations of the technique.

For efficiency, the package’s computational engine is written in C++ through the Rcpp package. The engine implements the Kalman filter and smoother for the state-space form underlying the HPF. All formulae can be found in the article and in the vignette titled Formulae.

Users needing higher control over the estimation process can use the functions with names hpfj* and auto_hpfj*. However, the wrapper function hpj makes the process simpler and more stable. Indeed, good starting values are essential since the technique is based on a complex numerical optimisation. Inside the hpj function, we linearly transform the time series so that the starting values in the hpfj* functions are expected to work fine, and then we apply the anti-transform to all the results so that they refer to the original time series. The functions hpfj* and auto_hpfj* allow full control over the starting values and penalties.

While hpfj* and auto_hpfj* accept only time series stored in a numerical vector (any type of time series object is cast to numerical vector), the hpj function accepts all the main time series objects used in R (ts, zoo, xts, and timeSeries). All the time series generated by the function receive the same class and dates as the original time series.

Basic usage

The function is called as

hpj(y, maxsum = NULL, lambda = NULL, xreg = NULL, ic = c("bic", "hq", "aic", "aicc"), scl = 1000)

where y is the time series to be filtered, maxsum is inverse penalty, when zero no jumps are allowed, the larger maxsum the more and larger breaks are allowed, lambda is the smoothing parameter of the HPF, xreg is a matrix of regressors (we are still woring on this, so, at the moment it is not yet implemented), and ic is the information criterion used to select the penalty when maxsum = NULL. When lambda = NULL, the smoothing parameter is estimated by quasi maximum likelihood. The parameter lambda can be a positive number or one of the following strings (in parenthesis the value of the smoothing parameter corresponding to the string): daily (110,930,628,906), weekly (45,697,600), monthly (129,600), quarterly (1,600), annual (6.25).

The function returns an object of class hpj, which is a list with the following slots:

y the input time series;
maxsum the maximum sum of additional standard deviations;
lambda the smoothing constant of the HP filter;
pars vector of estimated parameters (sigma_slope, sigma_noise, gamma);
hpj: the time series of the HP filter with jumps;
hpj_std the time series of the HP filter with jumps standard deviations;
breaks vector of indices of the detected breaks;
sigmas vector of additional standard deviations of the level disturbance;
xreg matrix of regressors;
df model’s degrees of freedom;
loglik value of the log-likelihood at maximum;
ic vector of information criteria (aic, aicc, bic, hq);
opt the output of the optimization function (nloptr);
call the call to the function.

The methods print and plot are available for the hpj class. The plot method can be customised:

plot(x, prob = NULL, show_breaks = TRUE, main = "original + filter", use_ggplot = TRUE, ...)

where x is the object of class hpj, prob is the coverage of the confidence interval for the filter, which is not plotted when prob = NULL, show_breaks is a logical indicating whether to show the breaks in the plot, main is the title of the plot, use_ggplot is a logical indicating whether to use ggplot2, and ... are additional arguments passed to the plot function when the standard plot is used (that is when use_ggplot = FALSE).

Examples

Simulated time series

library(jumps)
set.seed(2025)
n <- 100

# simulated smooth trend
mu <- 100*cos(3*pi/n*(1:n)) - ((1:n) > 50)*n - c(rep(0, 50), 1:50)*10
# simulated time series
y <- mu + rnorm(n, sd = 20)

# HP filter with jumps with estimated lambda and fixed penalty (maxsum = 50)
hpj_sim <- hpj(y, maxsum = 50)

print(hpj_sim)
#> Hodrick-Prescott filter with jumps
#> Call:
#>   hpj(y, maxsum = 50) 
#> Parameters:
#>   sd(slope) = 1.848931
#>   sd(noise) = 20.67477
#>   gamma     = 0.03659672
#>   lambda    = 125.0375
#>   maxsum    = 50
#> Model's degrees of freedom:  12.44784 
#> Log-likelihood:  -492.3048 
#> Information criteria:
#>   AIC  = 1009.505
#>   AICc = 1013.373
#>   BIC  = 1041.934
#>   HQ   = 1022.63
#> Break dates:    
#> 1 51

plot(hpj_sim)

plot(hpj_sim, prob = 0.95)

plot(hpj_sim, use_ggplot = FALSE)

Nile time series

# HP filter with jumps with estimated lambda and automatically selected penalty
hpj_nile <- hpj(Nile)

print(hpj_nile)
#> Hodrick-Prescott filter with jumps
#> Call:
#>   hpj(Nile) 
#> Parameters:
#>   sd(slope) = 1.692275e-07
#>   sd(noise) = 127.634
#>   gamma     = 3.428068e-14
#>   lambda    = 5.688412e+17
#>   maxsum    = 270.764
#> Model's degrees of freedom:  2.949825 
#> Log-likelihood:  -644.0627 
#> Information criteria:
#>   AIC  = 1294.025
#>   AICc = 1294.268
#>   BIC  = 1301.71
#>   HQ   = 1297.135
#> Break dates:      
#> 1 1899

plot(hpj_nile, main = "Nile river flow")
#> Don't know how to automatically pick scale for object of type <ts>. Defaulting
#> to continuous.

plot(hpj_nile, prob = 0.95, main = "Nile river flow")
#> Don't know how to automatically pick scale for object of type <ts>. Defaulting
#> to continuous.

plot(hpj_nile, use_ggplot = FALSE, main = "Nile river flow")

Employment in Italy

data("employed_IT")
y <- window(employed_IT[, "Y25.29"], start = c(2009, 1))
hpj_emp <- hpj(y, scl = "original")

print(hpj_emp)
#> Hodrick-Prescott filter with jumps
#> Call:
#>   hpj(y, scl = "original") 
#> Parameters:
#>   sd(slope) = 3.347636
#>   sd(noise) = 21.13837
#>   gamma     = 17.3066
#>   lambda    = 39.87185
#>   maxsum    = 76.59153
#> Model's degrees of freedom:  12.26254 
#> Log-likelihood:  -301.9286 
#> Information criteria:
#>   AIC  = 628.3823
#>   AICc = 635.1959
#>   BIC  = 654.267
#>   HQ   = 638.5267
#> Break dates:         
#> 1 2010.25
#> 2 2012   
#> 3 2013.75
#> 4 2019.75
#> 5 2020.25
plot(hpj_emp, main = "Millions of employed in Italy: age 25-29")

Smoothing splines with jumps

Besides the HP filter, the package fits cubic smoothing splines that allow a small number of discontinuities at unknown locations, following Pelagatti, Maranzano and Gatto Ronchieri (2026).

The idea is to leave the classical roughness penalty untouched and instead deform the domain over which it is computed. Every pair of consecutive observations is granted a non-negative domain warp \(\gamma_i \ge 0\), an amount of fictitious domain length inserted between them, so that the latent spacing becomes \(\tau_i = \delta_i + \gamma_i\) with \(\delta_i = x_{i+1} - x_i\). The spline is fitted, and its roughness penalised, on this warped domain. Where the signal is smooth, warping buys nothing and the optimal \(\gamma_i\) is zero, so the fit coincides with an ordinary smoothing spline; where the signal jumps, the same vertical change can be spread over a longer warped interval, which lowers its curvature cost dramatically. A budget \(\sum_i \gamma_i \le M\) keeps the warping concentrated on the few locations where a jump is actually present.

Four estimators are provided. ssj_mle() casts the model in state-space form and estimates \((\sigma, \lambda, \boldsymbol{\gamma})\) by maximum likelihood through the Kalman filter; ssj_kfs() does the same with \(\lambda\) held fixed; ssj_atw() is the fast alternating domain-warping algorithm, of complexity \(O(n \log n)\) per iteration; and ssj_fix() implements the earlier variant in which the additional disturbance variances, rather than the spacings, are the free parameters. The auto_* wrappers select the budget \(M\) over a grid by information criteria, EBIC by default. Two further wrappers, auto_ssj_atw_gcv() and auto_ssj_atw_cv(), also choose \(\lambda\), so that the alternating algorithm does not need a pilot fit from ssj_mle().

The strength of the EBIC search penalty, ebic_xi, has a calibrated default that differs by estimator – 0.3, 0.15 and 0.25 for ssj_mle(), ssj_kfs() and ssj_atw(). They differ because the three are different procedures, not because they fit different models; see the “Choosing ebic_xi” section of each help page.

A simulated example

We use the test signal of Section 3 of the paper: a Bessel function of the first kind carrying three discontinuities, at \(x = 0.3\) (of size \(+0.3\)), \(x = 0.4\) (\(-0.4\)) and \(x = 0.6\) (\(-0.6\)), observed at 100 random locations with Gaussian noise.

g1 <- function(x) {
  besselJ(20 * x, nu = 1) + x * (x >= 0.3 & x <= 0.4) - x * (x >= 0.6)
}
true_disc <- c(0.3, 0.4, 0.6)

set.seed(2026)
n <- 100
x <- sort(runif(n))
y <- g1(x) + rnorm(n, sd = 0.05)

Two of the estimators are shown here, because they are the two a user would actually choose between, and they reach the same answer by different routes. auto_ssj_mle() estimates \((\sigma, \lambda)\) by maximum likelihood through the Kalman filter and lets EBIC pick the warping budget; auto_ssj_atw_cv() runs the fast alternating algorithm, taking \(\lambda\) from a cross-validation of its own and again leaving the budget to EBIC. The second is self-contained: unlike the two-step recipe used in the paper’s simulation study it needs no pilot fit from the likelihood-based estimator. The ebic_xi defaults are calibrated and differ by estimator; ?ssj_atw explains why, and nothing needs overriding here.

fit_mle <- auto_ssj_mle(y = y, x = x, ic = "ebic")
fit_mle
#> Smoothing spline with jumps
#> Method: SSJ-MLE: domain warping, state-space maximum likelihood
#> Call:
#>   auto_ssj_mle(y = y, x = x, ic = "ebic")
#> 
#> Parameters:
#>   sigma    = 27.27
#>   lambda   = 3.534e-06
#>   budget   = 0.09764   (maximum sum of the domain warps)
#> 
#> Observations: 100
#> Effective degrees of freedom: 26.17
#> Log-likelihood: 93.58
#> Information criteria:
#>   AIC   = -134.8
#>   AICC  = -115.3
#>   BIC   = -66.66
#>   HQ    = -107.2
#>   EBIC  = -39.09
#> 
#> Declared discontinuities: 5
#>         x     jump     warp
#> 32 0.2984  0.07775 0.008003
#> 33 0.3178  0.23885 0.006113
#> 44 0.4107 -0.25526 0.030401
#> 45 0.4124 -0.03780 0.006045
#> 64 0.6017 -0.42911 0.039998

fit_atw <- auto_ssj_atw_cv(y = y, x = x, ic = "ebic")
fit_atw
#> Smoothing spline with jumps
#> Method: SSJ-ATW: alternating domain-warping algorithm
#> Call:
#>   auto_ssj_atw_cv(y = y, x = x, ic = "ebic")
#> 
#> Parameters:
#>   sigma    = 0.05091
#>   lambda   = 2.943e-06
#>   budget   = 0.09764   (maximum sum of the domain warps)
#> 
#> Observations: 100
#> Effective degrees of freedom: 27.3
#> Log-likelihood: 169.5
#> Information criteria:
#>   AIC   = -572.8
#>   AICC  = -551.3
#>   BIC   = -501.7
#>   HQ    = -544.1
#>   EBIC  = -478.8
#> 
#> Declared discontinuities: 5
#>         x     jump     warp
#> 32 0.2984  0.08228 0.008554
#> 33 0.3178  0.23757 0.005578
#> 44 0.4107 -0.25368 0.027743
#> 45 0.4124 -0.04819 0.007752
#> 64 0.6017 -0.43060 0.037672

The two smoothing constants agree closely, which is what makes the cross-validated route a usable replacement for the pilot:

c(lambda_mle = fit_mle$pars[["lambda"]],
  lambda_cv  = fit_atw$pars[["lambda"]])
#>   lambda_mle    lambda_cv 
#> 3.533619e-06 2.943219e-06

Both find all three discontinuities: SSJ-MLE at 0.298, 0.318, 0.411, 0.412, 0.602 and SSJ-ATW-CV at 0.298, 0.318, 0.411, 0.412, 0.602, against true locations 0.3, 0.4, 0.6. Recall that a break is dated at the first observation belonging to the new regime, so the reported location is the first sampled point after the discontinuity rather than the discontinuity itself.

There are more declarations than discontinuities, and it is worth seeing why. The warping budget is spread over the two intervals either side of a transition rather than concentrated on one, so a single jump can clear the declaration threshold twice at adjacent locations; the two reported sizes then add up to roughly the true one. This is a property of the reporting rule, which treats each interval separately, not of the fit – the curve itself is unaffected, and ssj_breaks() accepts stricter frac_total and frac_delta if you would rather see one declaration per transition. Judge the fit by the curve and the warp profile; treat the table as a guide to where the transitions are.

The plot method draws the fitted curve over the data and the estimated domain warps below, which is where the jumps become visible: spikes standing out of an otherwise flat sequence. Since version 2.0.0 it builds the figure with ggplot2 and returns it, so the two panels are facets of one object and the result can be themed or saved like any other plot; pass use_ggplot = FALSE for base graphics instead.

plot(fit_mle, xlab = "x", ylab = "y")

plot(fit_atw, xlab = "x", ylab = "y")

The two are hard to tell apart, which is the point: the alternating algorithm reproduces the maximum-likelihood fit at a fraction of the cost, and it no longer needs that fit in order to run.

Superimposing the true signal confirms that the fit is smooth between the breaks and sharp at them.

plot(x, y, pch = 16, cex = 0.6, col = "grey55", xlab = "x", ylab = "y")
lines(x, g1(x), col = "grey20", lty = 2)
lines(fit_atw$x, fit_atw$smoothed_level, col = "firebrick", lwd = 2)
abline(v = ssj_breaks(fit_atw)$x, lty = 3, col = "firebrick")
legend("bottomleft", bty = "n", cex = 0.85,
       lty = c(2, 1), lwd = c(1, 2), col = c("grey20", "firebrick"),
       legend = c("true signal", "smoothing spline with jumps"))

Application: electricity demand and temperature

The package ships the data used in the electricity demand versus temperature application of inst/scripts/load_temperature.R. load_temp holds the full hourly series (Northern-Italy load from the GME day-ahead market, Milan temperature from the Open-Meteo archive, 2023-01-01 to 2025-09-30), while demand_temp is the homogeneous analysis sample obtained from it: Tuesdays to Thursdays only, non-vacation months, the central hours 10-15, and no extreme loads or temperatures. Those filters strip out the weekday, holiday and time-of-day effects that would otherwise confound the relationship of interest.

data(demand_temp)
str(demand_temp)
#> 'data.frame':    1691 obs. of  2 variables:
#>  $ Temperature: num  8.8 10.8 12.2 12.9 12.9 12.1 10.9 12.6 12.6 13.6 ...
#>  $ Load       : num  25278 24460 23997 22126 22212 ...

Demand is a U-shaped function of temperature – heating on the left, cooling on the right – but the interesting question is whether the curve is smooth. Thermostat settings, air-conditioning switch-on thresholds and similar behavioural rules suggest that it need not be.

Here we use the alternating algorithm auto_ssj_atw() rather than auto_ssj_mle(), for two reasons. It is far faster: ssj_mle() estimates one warping parameter per observation, which on these 1691 points takes minutes, whereas the alternating algorithm takes a fraction of a second. And it collapses tied values of the covariate, which matters here because the sample has only 295 distinct temperatures: with ties left in place the budget can be spent warping the interval between two observations recorded at the same temperature, which is not identified. Since ssj_atw() needs \(\lambda\) to be supplied, we pass the value obtained from a preliminary auto_ssj_mle() fit – the pilot strategy used in the paper. That pilot must itself be run on the tie-averaged data, for the reason just given: fitted on the raw observations it returns \(\lambda = 154.5\), on the averaged ones \(\lambda = 254.3\).

fit <- auto_ssj_atw(y = demand_temp$Load, x = demand_temp$Temperature,
                    lambda = 254.3, ic = "ebic")
fit
#> Smoothing spline with jumps
#> Method: SSJ-ATW: alternating domain-warping algorithm
#> Call:
#>   auto_ssj_atw(y = demand_temp$Load, x = demand_temp$Temperature, 
#>        lambda = 254.3, ic = "ebic")
#> 
#> Parameters:
#>   sigma    = 1001
#>   lambda   = 254.3
#>   budget   = 32.3   (maximum sum of the domain warps)
#> 
#> Observations: 295
#> Effective degrees of freedom: 12.55
#> Log-likelihood: -12569
#> Information criteria:
#>   AIC   = 20364
#>   AICC  = 20364
#>   BIC   = 20432
#>   HQ    = 20389
#>   EBIC  = 20492
#> 
#> Declared discontinuities: 8
#>        x   jump  warp
#>  10  4.0 -550.6 2.350
#>  46  8.7 -297.2 2.029
#>  77 11.9 -485.8 3.282
#> 109 15.1 -212.7 1.725
#> 179 22.3  134.1 2.075
#> 224 27.0  298.9 1.655
#> 225 27.1  668.6 3.654
#> 288 34.0  880.2 5.041

All the ssj_* functions return an object of class ssj, with print and plot methods. The declared discontinuities reported above follow the rule of the paper: an interval counts as a discontinuity when its estimated domain warp \(\hat\gamma_t\) clears both 5% of the total warping and 30% of the mean sampling interval. The table itself is available through ssj_breaks(), which also accepts a single absolute threshold through its thresh argument.

ssj_breaks(fit)
#>   index    x     warp      jump
#> 1    10  4.0 2.349993 -550.5671
#> 2    46  8.7 2.028604 -297.1779
#> 3    77 11.9 3.281853 -485.7816
#> 4   109 15.1 1.725359 -212.6548
#> 5   179 22.3 2.074768  134.1489
#> 6   224 27.0 1.654665  298.8849
#> 7   225 27.1 3.653521  668.6330
#> 8   288 34.0 5.040918  880.2336

The plot method reproduces Figure 4 of the paper: the fitted curve over the data, with vertical lines at the declared breaks, and a lower panel with the estimated warps, which measure the strength of each jump. The points shown are the 295 distinct temperatures actually seen by the estimator rather than all 1691 observations. Pass prob to add a confidence band, whenever the estimator returns pointwise variances (ssj_mle(), ssj_kfs() and ssj_fix() do; ssj_atw() does not).

plot(fit, xlab = "Temperature (Celsius)", ylab = "Demand (MWh)")

The declared breaks fall into two groups, and their signs separate them cleanly. Above about 22 degrees they are upward and grow with temperature, peaking either side of 27 degrees, which is consistent with air conditioning being switched on once a comfort threshold is crossed; the paper’s Section 4.1 identifies the jump at 27 degrees as the most relevant one, and it appears here as the largest of the warps. Below 16 degrees they are downward and smaller, tracing the mild role of electric heating in Italy: demand falls as the weather warms towards the comfort zone. A classical smoothing spline would round all of these transitions off into a gentle curve.

Two of the declarations, at 27.0 and 27.1 degrees, are the same transition seen twice, for the reason discussed in the simulated example above. Raising ebic_xi, or the thresholds of ssj_breaks(), returns a shorter list concentrated on the cooling region alone.

The same data can be analysed with the CSSD algorithm of Storath and Weinmann (2024), which penalises the number of discontinuities directly instead of a warping budget. It is a different proposal, not one of ours, and the print method says so.

Its delta argument is the pointwise noise standard deviation, and supplying it is what makes the data-fit and roughness terms commensurate. Leaving delta = 1 on a series measured in MWh would make the data term dominate by several orders of magnitude and return an essentially interpolating fit.

# robust noise scale, immune to the jumps themselves
sig <- mad(diff(demand_temp$Load[order(demand_temp$Temperature)])) / sqrt(2)

fit_cssd <- ssj_cssd_dp(demand_temp$Load, demand_temp$Temperature,
                           p = 0.9, gamma = 10, delta = sig)
fit_cssd
#> Smoothing spline with jumps
#> Method: CSSD: cubic smoothing spline with discontinuities (Storath and Weinmann, 2024)
#> Call:
#>   ssj_cssd_dp(y = demand_temp$Load, x = demand_temp$Temperature, 
#>        p = 0.9, gamma = 10, delta = sig)
#> 
#> Parameters:
#>   p        = 0.9
#>   gamma    = 10
#> 
#> Observations: 295
#> 
#> Declared discontinuities: 2
#>        x   jump
#> 142 18.4  453.1
#> 224 27.0 1119.9

ssj_cssd_dp() returns an ssj object too, so the same print and plot methods apply. There is no lower panel here: CSSD has no domain warps to show, because its discontinuity set is part of the solution of the underlying combinatorial problem rather than a by-product of a continuous penalty.

plot(fit_cssd, xlab = "Temperature (Celsius)", ylab = "Demand (MWh)")

The fitted piecewise cubic spline can be evaluated at arbitrary temperatures with predict_cssd_dp(), which is what makes the size of a discontinuity directly readable:

predict_cssd_dp(fit_cssd, xnew = c(20, 26.9, 27.1, 30))
#> [1] 21617.71 21886.88 23038.07 23944.75

The two methods are quite different – one warps the domain under a budget constraint, the other penalises the number of breaks directly – yet both put a discontinuity at about 27 degrees, which is reassuring.

Introduction to the jumps package