This is a major release: the package now covers smoothing splines with discontinuities alongside the original Hodrick-Prescott filter with jumps, which is why the version jumps from 1.0 to 2.0.0.
plot.ssj() gains warp_frac, the share of the plotting height given to the
lower (warp) panel, default 0.3. The panel holds one value per interval and
most of them are zero, so the equal panels the method drew before spent half
the figure on it. Base graphics now use layout() rather than par(mfrow=),
which can only divide the height equally; the ggplot2 figure uses the
panel.heights theme element and so needs ggplot2 >= 4.0.0, warning if a
non-default warp_frac cannot be honoured. The figure is still a plain
ggplot, so it can be extended with + as before. The strip label of that
panel is shortened from "Domain warp" to "Warps", which no longer crowds the
axis beside it at the smaller height.
New family of functions for cubic smoothing splines with discontinuities
(CSSD), implementing the dynamic-programming algorithm of Storath and
Weinmann (2024): ssj_cssd_dp(), auto_ssj_cssd() (K-fold
cross-validation of p and gamma) and predict_cssd_dp().
ssj_cssd() provides a fast piecewise-constant (total-variation)
approximation for larger problems.
New family of smoothing splines with jumps based on time warping:
ssj_mle()/auto_ssj_mle() (smoothing parameter estimated by maximum
likelihood), ssj_kfs()/auto_ssj_kfs() (fixed smoothing parameter) and
ssj_atw()/auto_ssj_atw() (fast alternating time-warping algorithm).
All of them support the extended BIC (ebic) for selecting the jump budget.
The ssj_* functions now return an object of class ssj, with print and
plot methods, mirroring what hpj() already offered. The object is still a
plain list underneath, so existing code that reaches into its slots keeps
working; it gains a method code, the response y as the estimator saw it,
and the call. plot() reproduces Figure 4 of the paper: the fitted curve
over the data, vertical lines at the declared discontinuities and a lower
panel with the estimated domain warps; pass prob for a confidence band.
No title is drawn unless one is passed through main.
New predict() method for ssj objects. For the domain-warping estimators
it evaluates the cubic the model actually implies -- the state carries both
the level and the slope at every sampling site, so the fit is a cubic Hermite
in the latent coordinate -- instead of interpolating the fitted values. It was
validated against the exact model-implied interpolation, obtained by inserting
the query points as missing observations and re-running the Kalman smoother:
the two agree to 1e-12. Beyond the range of the data the fit is extended
linearly. CSSD objects are passed to predict_cssd_dp(), and the
piecewise-constant TV fit gets a step interpolation.
The ssj_* fitters gained a smoothed_slope slot, the derivative of the fit
with respect to the latent (warped) coordinate, which is what predict()
uses.
ssj_mle(), ssj_kfs(), ssj_fix() and ssj_atw(), and their auto_*
wrappers, gained a weights argument. The observation variance becomes
lambda*sigma^2/w_i (for the alternating algorithm the data term becomes
sum w_i (y_i-f_i)^2), so a point with weight w_i counts as the mean of
w_i independent observations. This makes it correct to collapse tied values
of x to their mean and pass the replicate counts as weights, which is the
right way to handle a covariate with many repeated values: verified to
reproduce the fit on the raw replicates exactly for ssj_atw() and to 1e-9
and 1e-5 for ssj_fix() and ssj_kfs(). Unit weights reproduce the previous
results exactly, so nothing existing changes. See the caveat in ?ssj_mle:
because that function estimates the noise level, collapsing loses the
within-group scatter that identifies it.
ssj_mle(), ssj_kfs() and ssj_fix() now collapse tied values of x to
their weighted mean automatically, as ssj_atw() already did, turning the
replicate counts into observation weights and -- for ssj_mle(), which
estimates the noise level -- feeding the within-group scatter back into the
likelihood. Previously they fitted the raw replicated data, where the
intervals between replicates have zero length and a warp placed on one of
them is a discontinuity between two observations at the same x, which is
not identified. On the electricity-demand application that made the fit
depend on the order of the rows, with effective degrees of freedom anywhere
between 8 and 79. The x, smoothed_level and gamma slots are now
returned on the distinct sites, as they already were for ssj_atw().
llt() and llt_delta() returned NA for the whole likelihood when the
FIRST observation was missing: the t = 0 initialisation lacked the
missing-value guard that the loop applies to every later observation. Any
series starting with NA was affected.
ssj_mle() and auto_ssj_mle() gained ss_within, the pooled within-group
sum of squares discarded when replicates are collapsed to their means.
Supplying it alongside the replicate counts in weights restores the term of
the likelihood that collapsing removes, so the noise level -- and hence
lambda -- is estimated as if the raw replicates had been supplied: on a test
case the estimated observation variance goes from 0.086 back to 0.021 against
a raw-data value of 0.021, at three different noise levels and to four or five
significant digits. It changes nothing at a fixed smoothing parameter, and
nothing when omitted.
ssj_atw() already collapsed tied values of x, but gave every collapsed
point unit weight, silently treating a mean of six observations as no more
informative than a single one. It now carries the replicate counts as
weights, which changes the effective degrees of freedom on data with ties
(3.52 rather than 2.84 on a test case).
New function ssj_breaks(), which applies the detection rule of the paper --
a warp counts as a discontinuity when it clears both 5% of the total warping
and 30% of the mean sampling interval -- and returns the declared breaks with
their location, warp and jump size. Breaks are dated at the first observation
of the new regime. An absolute threshold can be given through thresh, as in
Section 4.2 of the paper.
New dataset load_temp with the hourly electricity load of Northern Italy
and the temperature of Milan, and demand_temp, the filtered analysis
sample used in the demand-versus-temperature application. Both are shipped
as .rda, so neither the 11 MB raw market file nor the Open-Meteo API call
is needed to reproduce the analysis. The script that derives demand_temp
is in data-raw/demand_temp.R.
The introduction vignette gained two sections: one on the smoothing splines with jumps, and one working through the electricity demand versus temperature application with base graphics only.
The ebic_xi defaults are now calibrated rather than taken from theory,
and they differ by estimator: 0.3 for ssj_mle(), 0.15 for ssj_kfs() and
0.25 for ssj_atw(), against 1, 1 and 0.3 before. Each is the value that
minimises worst-case regret in root mean squared error over a factorial
design of 48 cells -- four sample sizes, three noise levels and four counts
of discontinuities, 17,400 replications -- shipped as
inst/scripts/xi_calibration.R.
Fits selected by auto_ssj_mle(), auto_ssj_kfs() and auto_ssj_atw() will
therefore change: the smaller penalties admit more discontinuities. Pass
ebic_xi explicitly to reproduce earlier results.
The three differ for real reasons rather than by accident. ssj_kfs() sits
below ssj_mle() because auto_ssj_kfs() holds lambda fixed across the
budget grid while auto_ssj_mle() re-estimates it, and being able to re-tune
the smoothness gives more scope to overfit. For ssj_atw() the optimum falls
steadily with the sample size -- roughly 0.32, 0.22, 0.15, 0.12 at
n = 50, 100, 200, 400 -- while the total charge per declared discontinuity
stays flat at about 3.9, which yields the rule
xi(n) = (3.93 - 0.419*log(n)) / (2*log(n-1)). That rule is documented and
is worth using when n is far from a few hundred; the constant is the
default because it behaves acceptably everywhere and a default that changes
with the data is easy to overlook.
All three sit below the Chen and Chen (2008) consistency floor of about 0.5,
which is expected rather than a defect: df is the trace of the smoothing
matrix and grows with each activated warp, so the plain BIC term already pays
part of the search cost -- about 0.08, 0.17 and 0.21 on the xi scale for
the three estimators.
plot() for ssj objects gains use_ggplot and defaults to it, as
plot.hpj() already did. The figure is built with ggplot2 and returned
rather than drawn, so it can be extended, themed or saved with ggsave();
the fitted curve and the warp profile become facets of a single object rather
than two separate plots, so the x axes stay aligned and one value comes back.
Pass use_ggplot = FALSE for base graphics, which draws directly and returns
NULL.
One consequence worth knowing: the plot is now a value, so inside a loop or a
script it has to be printed to appear, exactly as any ggplot2 figure does.
New auto_ssj_atw_gcv() and auto_ssj_atw_cv(), which choose the smoothing
constant for the alternating algorithm without recourse to ssj_mle(), then
hand it to auto_ssj_atw() for the warping budget as before. The first
minimises generalized cross-validation over lambda at zero budget, at one
fit per candidate; the second cross-validates over lambda and the budget
jointly and keeps only lambda. Both accept weights, warn when the
selected value is an endpoint of the grid, and return the search path.
The two are not equivalent. Cross-validation is the reliable one: the
lambda it selects tracks the maximum-likelihood value closely and its
boundary bound in 0 to 23% of samples on the calibration design. Generalized
cross-validation is around two orders of magnitude cheaper but degenerates
exactly where the method is aimed -- with three clear jumps and
sigma = 0.05 its floor bound in 83% of samples at n = 100 and in every
sample at n = 200, because at zero budget the spline can only represent a
discontinuity by bending, so the criterion improves without limit as lambda
falls. Accuracy is nevertheless within a few percent of the alternatives,
since the floor is a serviceable value; but it is a floor, not a minimum.
The default lambda grid is built in the normalised coordinates
ssj_atw() works in and mapped back, so it follows the scale of x -- a
grid fixed in absolute terms would suit only one scaling, lambda having the
units of x^3.
ssj_atw() gains a pars slot, holding the unbiased residual standard
deviation and the lambda it was given, so the smoothing constant can be
recovered from a fit as it already could from the other estimators.
auto_ssj_cssd() searched for the jump penalty on a grid fixed in absolute
terms, so its answer depended on the units of the response. gamma is a
charge per discontinuity levied against the data-fit and smoothness terms,
both of which scale with the square of the units of y, but the coarse search
grid was exp(seq(log(0.5), log(100), length.out = 20)) whatever the data. On
the electricity-demand example (demand_temp, sd 1523 MWh, noise variance
about 1.3e6) the whole grid lay below 1e-4 of a single squared residual, where
a discontinuity costs nothing: the fit reported 146 discontinuities among 295
sites, at a cross-validated error 2.8 times worse than the optimum the search
never reached, and dividing y by 1000 gave none at all. On a signal of unit
scale the same grid is the opposite extreme -- at sigma = 0.05 on the design of
the simulation study the cross-validated optimum falls below the grid's floor,
so the estimator was held to fewer discontinuities than cross-validation
wanted. The search is now conducted on y standardised to unit standard
deviation, with gamma mapped back to the units of y on return, and the
grid spans six orders of magnitude. Fits are now invariant to the units of the
response (agreeing to 4e-11 MWh across a 1e6 change of units). A selected
penalty at the edge of the search range now warns, and the range is exposed as
the new gamma_grid argument. gamma_init is interpreted in the units of y
and cv_score is reported in the units of y^2.
Scaling the search grid alone would not have sufficed. The search is
parameterised as gamma = p*q/(1-q) with the optimiser working in logit(q),
so the reachable range of gamma depends on the units of y even when the
grid is scaled with them: on the demand data in MWh that route gives the right
answer, but at a thousandfold larger scale it drifts by two orders of
magnitude and at a millionfold larger scale 1 - q underflows and the search
fails outright. Standardising holds logit(q) at the same value for every
scaling. q is now also clamped below 1, so an extreme gamma_grid returns
the classical spline rather than failing.
This changes results. On the electricity-demand data auto_ssj_cssd() now
declares no discontinuities rather than 146; on the simulation design at
sigma = 0.05 it declares markedly more than before. Any CSSD result obtained
with an earlier version should be recomputed.
auto_ssj_atw_cv() cross-validated optimistically on designs with tied
covariate values. Folds were drawn over the observations, so replicates of
the same x landed in different folds. Since ssj_atw() collapses tied x
inside whatever sample it is handed, the training set then contained siblings
of every held-out point and pinned the fit at that very covariate value: the
criterion measured the error of predicting a fresh observation at an
already observed x, not at a new one. On tied designs every held-out
observation was affected, and the cross-validated error understated the
quantity of interest by about a third under heavy replication (median 0.043
against 0.064 over six designs with 25 sites and 10-30 replicates each).
Folds are now drawn over the distinct values of x, so all replicates at a
site are held out together; nfolds is bounded by the number of distinct
values rather than by length(y), and a user-supplied folds that splits a
site now warns. Designs with no ties are unaffected -- the two partitions
coincide there -- so no previously published result changes.
ssj_atw() and the criteria behind it used the wrong effective sample size
when the observation weights were not all one, which includes every fit in
which tied values of x were collapsed, since the replicate counts then
become the weights. mse, gcv, aic, aicc, bic, hq, the reported
log-likelihood and the reported sigma all divided the weighted residual sum
of squares by sum(weights), as though it held one degree of freedom per
underlying observation. It does not: y_i is the weighted mean of its
replicates and has variance sigma^2/w_i, so w_i*(y_i - f_i)^2 has
expectation sigma^2 and the weighted rss carries one degree of freedom per
fitted point. The within-group scatter of the replicates is removed by the
collapse and never reaches the routine at all.
Two consequences, both now fixed by using n. The reported sigma was far
too small -- 0.050 against a true 0.200 on a 400-observation, 40-site test,
where it now returns 0.193. And the data-fit term of every criterion was
inflated relative to its penalty by sum(weights)/n, so the penalty was that
much too weak: on the electricity-demand data (1691 observations at 295
distinct temperatures, a ratio of 5.7) the extended BIC fell monotonically to
the largest budget on the grid and auto_ssj_atw_cv() declared 20
discontinuities, against 4 for auto_ssj_mle() on the same data. It now
declares 2, with an estimated noise standard deviation of 1096 MWh against
1120 MWh from the within-site scatter, and its discontinuity at 27 degrees is
a subset of the four that auto_ssj_mle() finds.
Designs with no ties and unit weights are unaffected, n and sum(weights)
being equal there, so no previously published result changes. The reported
log-likelihood now also carries the 0.5*sum(log(w)) Jacobian of the
collapse, which is free of the budget and so cannot affect selection.
ssj_kfs() and ssj_fix() frequently failed to converge, returning a
scale parameter many orders of magnitude too large and a log-likelihood far
below the attainable one. Both started the optimiser at a fixed sdy/10,
which bears no relation to the optimum: when lambda is small the optimal
sigma is large -- 37 against a start of 0.1 on the test case -- and the
optimiser diverged. Against ssj_mle() run on the same data at the same
lambda, where the two fit exactly the same model and must agree,
ssj_kfs() reached a worse optimum in 37 of 40 test replications, with
gaps of 800 log-likelihood units; ssj_fix() was worse still, at around
7,000,000 units. auto_ssj_kfs() propagated the failure along the whole
grid through its warm start.
Both now start from the profile solution for sigma, which needs a single
filter pass: with lambda fixed, sigma is a pure scale, so the profile
log-likelihood is maximised at mean(v^2/F) evaluated at sigma = 1. That
lands essentially on the optimum, and it is followed by a short multi-start
in the same spirit as the grid of lambda starts already used by
ssj_mle(). ssj_kfs() now matches ssj_mle() to 1e-10 at the same
lambda, and no longer falls short in any of the 40 replications.
Users of ssj_kfs(), auto_ssj_kfs(), ssj_fix() or auto_ssj_fix()
should re-run anything fitted with an earlier version.
ssj_mle(), ssj_kfs() and ssj_atw() reported the log-likelihood of an
internally rescaled copy of y, not of the data as supplied, and with it
the absolute values of aic, aicc, bic, hq and ebic. All three
divide y by its standard deviation before fitting, and the Jacobian of that
map was never added back, so the reported value sat below the true one by
nobs * log(sd(y)) -- 12.5 units on a 60-point test case. The symptom was
that multiplying y by 1000 left the reported log-likelihood unchanged,
where it should fall by nobs * log(1000).
All three now report the likelihood on the scale of the data supplied.
ssj_mle() uses nobs + df_within factors when the within-group scatter has
been restored, that being the number of raw observations. ssj_fix() and the
hpfj*() family never rescaled y and are unchanged.
This does not change any model selection. The offset is identical at
every maxsum for a given data set, so it cancels in the comparisons the
auto_* functions make; the budget they choose, the fitted curve and df
are all unaffected. What changes is the reported loglik and the absolute
level of the criteria, which are now comparable across estimators: on the
identical fitted model ssj_kfs() and ssj_fix() now report a bic of
41.088 against 41.060, where before it was 66.157 against 41.060.
A caveat that remains, and is not a defect of the above: these fitters
approximate the diffuse initial conditions with a large proper prior
(1e5 times the data variance), and the reported likelihood depends on that
arbitrary multiplier -- by about 0.5*log(10) per decade per diffuse state.
Absolute criterion values therefore carry an additive constant fixed by that
convention. It is common to all of ssj_mle(), ssj_kfs() and ssj_fix(),
but ssj_fix() does not rescale x, so its slope prior differs slightly
from the others'; that is the residual 0.03 above. Exact diffuse
initialisation would remove the constant altogether.
The effective degrees of freedom were wrong for every Kalman-based
estimator (ssj_mle(), ssj_kfs(), ssj_fix(), hpfj(), hpfj_fix(),
hpfjx()), and with them aic, aicc, bic, hq and ebic, and hence the
jump budget selected by the auto_* functions. df is the trace of the
smoothing matrix, obtained by summing the observation weights computed inside
llt() and llt_delta() with the recursion of Koopman and Harvey (2003).
Those weights carried the wrong time index on the backward variance smoother:
N was taken one step too early. Since N grows going backwards, the
subtracted terms overshot and the weights could fall well outside (0, 1] and
sum to a negative number -- reported df values such as -579 on a smooth fit
of 60 observations. The weights are now exactly the diagonal of the hat
matrix, verified against a central-difference computation of it for equally
and unequally spaced designs, with and without observation weights, and with
missing values (which correctly receive weight zero). hpj objects return
the corrected weights in their weights slot. ssj_atw() was never
affected: its df comes from the trace of A^-1 S in the alternating
algorithm.
Note for users reproducing published results: the formula in Appendix B.5 of Maranzano and Pelagatti (2025) carries the same index error, and information criteria computed with earlier versions of this package will differ.
hpfj_fix() reported the wrong value for pars["gamma"]: it squared the
first additional standard deviation instead of the estimated sqrt(gamma).
plot.hpj() drew an inverted confidence band when use_ggplot = TRUE and
prob was given (ymin and ymax were swapped, because the negative
quantile from qnorm() was not taken in absolute value).
plot.hpj() now falls back to base graphics, with a warning, when
ggplot2 is not installed, and only uses xts when that package is
available.
BIC.hpj() failed with "argument to 'which' is not logical" when one of the
objects passed through ... had more than one class. It now uses
inherits().
predict_cssd_dp() returned its values in the order of sort(xnew)
rather than in the order of xnew, silently mis-aligning predictions when
the query points were not sorted.
predict_cssd_dp() returned exactly 0 for query points falling between
the last data site of one segment and the first data site of the next, that
is, precisely around each detected discontinuity. Segment membership is now
determined by cutting at the midpoint of each gap.
auto_ssj_cssd() could not be called without an explicit x argument,
although its body handled x = NULL; x now defaults to NULL. It also
sorts and de-duplicates x before cross-validation, as the other CSSD
functions do.
ssj_atw() produced an all-NaN fit when y contained missing values,
which its documentation claimed to allow. It now raises an informative error
pointing to ssj_kfs()/ssj_mle().
hpfj_fix() used var(y) without na.rm = TRUE and initialised the filter
at y[1], both of which failed on series with missing values.
The examples of ssj_atw() and auto_ssj_atw() failed with "'x' and 'y'
lengths differ", because those functions collapse tied x values. All the
spline examples now plot against the x slot of the returned object.
The paper behind the smoothing splines with discontinuities is now cited
throughout: Pelagatti, Maranzano and Gatto Ronchieri (2026), "Smoothing
splines with some discontinuities", arXiv:2610.33137,
doi:10.48550/arXiv.2610.33137. It appears in DESCRIPTION, the package
help page, inst/CITATION, the README, the introduction vignette and the
@references of every function implementing one of these methods --
ssj_mle(), ssj_kfs(), ssj_atw(), ssj_fix() and their auto_*
wrappers, auto_ssj_atw_gcv(), auto_ssj_atw_cv(), ssj_breaks(), the
print, plot and predict methods for ssj objects, and the CSSD
functions, where it sits alongside Storath and Weinmann (2024).
The reference for the Hodrick-Prescott filter with jumps -- Maranzano and
Pelagatti (2025), Economic Modelling 150, 107132,
doi:10.1016/j.econmod.2025.107132 -- is now also given in the
@references of every function of that family, which previously cited it
only at package level: hpj(), hpfj(), hpfj_fix(), hpfjx(), their
auto_* wrappers and the print, plot, logLik, nobs and BIC
methods for hpj objects. The two families are now documented
symmetrically, each function pointing at the paper introducing it.
inst/CITATION referred to ssj_cssd_potts(), which was renamed
ssj_cssd_dp() in this release, and gained a header for the
Hodrick-Prescott entry so that each of the three references says which
functions it covers.
The NAMESPACE file is now generated by roxygen2, so that exports and
imports can no longer drift apart from the source (hpfj_fix() was exported
by hand but had no @export tag).
Argument and return-value documentation has been made uniform across all
functions, and several mismatches with the code have been fixed: the
maxit/maxit_nm argument of auto_ssj_cssd(), the std_devs slot of
hpj() objects (actually named sigmas), the documented default of the
ic argument of auto_ssj_atw(), and the dangling link to a
non-existent cssd_cv_lambda().
predict_cssd_dp() extrapolated beyond the end knots of a segment with
that segment's end cubic. The reference implementation of Storath and Weinmann
extends each segment linearly from its end knots (linext_pp), which is the
convention CSSD is defined with; a cubic extrapolation diverges quickly across
the gaps around a discontinuity. It now extends linearly.
predict() on a domain-warping fit extrapolated beyond the data using the
physical slope just inside the end interval, that is the latent slope times
that interval's dilation tau/delta. The warps live strictly between
observations, so outside the sample the dilation is 1. On a boundary interval
carrying a warp the factor can be several tens -- one simulated replication
had tau/delta = 40, turning a slope of 6.6 into 267 and predicting -11.7
where the signal was 0. It now extrapolates with the bare latent slope.
The documentation of ssj_cssd() no longer claims that the difference
operator is weighted by 1/dx; the implementation deliberately uses the
unweighted operator.
hpj() now raises an explicit error when xreg is supplied, instead of
documenting an argument that was not implemented.
The default maxsum of ssj_mle() was sd(y)/mean(diff(x)), a formula
inherited from ssj_fix(), where the budget constrains disturbance standard
deviations. For the domain-warping estimators the budget is expressed in the
units of x, so the old default was on the wrong scale: on faithful it
evaluated to about 1053 against an x-range of 3.5, leaving the warping
effectively unconstrained. It is now max(x) - min(x), matching both
ssj_kfs() and the upper end of the grid searched by auto_ssj_mle().
Calls that pass maxsum explicitly are unaffected.
ssj_mle(), ssj_kfs() and ssj_fix() gave different results depending on
the order in which the observations were supplied, even when the covariate
had no ties. They computed var(diff(y)) on the unsorted response and used
it to scale the diffuse prior on the slope, so a permutation of the rows
changed the starting point of the Kalman filter (by a factor of 6.5 on the
demand_temp data). It is now computed after sorting by x. Results for
input that was already sorted by x are unchanged.
ssj_kfs() silently ignored its own default maxsum. The default,
max(x) - min(x), refers to x, and R forces such a default lazily, at its
first use; that happened only after x had been rescaled to [0, 1], so the
budget evaluated to 1 instead of the range of x. The defaults of
ssj_kfs(), ssj_mle() and ssj_atw() are now forced while x is still on
its original scale. Calls that pass maxsum explicitly were never affected.
The examples of ssj_atw(), ssj_kfs() and ssj_mle() used maxsum = 150
on the faithful data. For the domain-warping estimators maxsum is a
budget in the units of x, whose range there is about 3.5, so the budget was
some forty times the whole domain and the resulting fits were degenerate (the
Kalman variant even reported negative effective degrees of freedom). The
examples now use a budget of 1.
The reference for the HP filter with jumps has been updated to the published
article: Maranzano and Pelagatti (2025), Economic Modelling 150, 107132,
doi:10.1016/j.econmod.2025.107132, replacing the FEEM working paper in
DESCRIPTION, inst/CITATION, the package help page, the README and the
introduction vignette.
Package-level documentation is now generated from R/jumps-package.R, and
the reference papers, replication data and rendered vignettes are no longer
shipped inside the source tarball.