Package 'exuber'

Description

Augment and then join objects.

Usage

augment_join(x, y, ...)

Arguments

Value

A tibble joining the augmented statistics of x with the augmented critical values of y; see the methods for the exact columns.

Description

Augment and then join the full statistic sequences of a radf_obj with the critical-value sequences of a radf_cv, with one row for each observation. This is the table autoplot.radf_obj is built on.

Usage

## S3 method for class 'radf_obj'
augment_join(x, y = NULL, trunc = TRUE, ...)

Arguments

Value

A tibble::tibble()

Examples

rsim_data <- radf(sim_data, minw = 20)
cv <- radf_wb_cv(sim_data, minw = 20)

# Join the statistic path and the critical-value path. This is the table autoplot() is built on
aj <- augment_join(rsim_data, cv)
aj

# Reproduce (a simplified version of) autoplot()'s own bsadf-vs-crit line plot
library(ggplot2)
aj %>%
  dplyr::filter(sig == 95, stat == "bsadf") %>%
  tidyr::pivot_longer(c(tstat, crit), names_to = "series") %>%
  ggplot(aes(index, value, col = series)) +
  geom_line() +
  facet_wrap(~id, scales = "free")

Description

Returns the full critical-value sequences (badf_cv, bsadf_cv) of a radf_cv, one row per observation; see tidy.radf_cv for the scalar critical values.

Usage

## S3 method for class 'radf_cv'
augment(x, format = c("wide", "long"), trunc = TRUE, ...)

Arguments

Value

A tibble::tibble()

Description

Returns the full test-statistic sequences (badf, bsadf) of a radf_obj, one row per observation; see tidy.radf_obj for the scalar statistics.

Usage

## S3 method for class 'radf_obj'
augment(x, format = c("wide", "long"), panel = FALSE, trunc = TRUE, ...)

Arguments

Value

A tibble::tibble()

Description

Plots the two input series over the sample, with the test statistic and the lag in the title.

Usage

## S3 method for class 'cobubble_test_obj'
autoplot(object, ...)

Arguments

Value

A ggplot

See Also

cobubble_test

Description

Plots the time-varying contagion coefficient $\delta_2(r)$ against the sample fraction $r$, with a zero reference line.

Usage

## S3 method for class 'contagion_reg_obj'
autoplot(object, ...)

Arguments

Value

A ggplot

See Also

contagion_reg

Description

Plots the CUSUM-type statistic path against its critical values, with one panel for each series. For the two-sided CUSUM-of-squares tests, the sup path is shown against the upper critical value and the inf path against the lower one.

Usage

## S3 method for class 'cusum_test_obj'
autoplot(object, ...)

Arguments

Value

A ggplot

See Also

cusum_test

Description

Plots each series with vertical markers at the estimated origination, collapse and recovery dates.

Usage

## S3 method for class 'dating_hls_obj'
autoplot(object, ...)

Arguments

Value

A ggplot

See Also

dating_hls

Description

Plots each series with vertical markers at every estimated origination, collapse and recovery date.

Usage

## S3 method for class 'dating_hlw_obj'
autoplot(object, ...)

Arguments

Value

A ggplot

See Also

dating_hlw

Description

Plots each series with vertical markers at the estimated origination and collapse dates.

Usage

## S3 method for class 'dating_knp_obj'
autoplot(object, ...)

Arguments

Value

A ggplot

See Also

dating_knp

Description

Plots each series with vertical markers at the estimated origination, collapse and, where estimated, recovery dates.

Usage

## S3 method for class 'dating_pdc_obj'
autoplot(object, ...)

Arguments

Value

A ggplot

See Also

dating_pdc

Description

Takes a ds_radf object and returns a ggplot2 object with a geom_segment layer.

Usage

## S3 method for class 'ds_radf'
autoplot(object, trunc = TRUE, ...)

Arguments

Value

ggplot

Examples

sim_data_wdate %>%
  radf() %>%
  datestamp() %>%
  autoplot()

# Change the colour manually
sim_data_wdate %>%
  radf() %>%
  datestamp() %>%
  autoplot() +
  ggplot2::scale_colour_manual(values = rep("black", 4))

Description

Bar chart of the LBI statistic for each series against its critical value. Series that exceed it are flagged as detected.

Usage

## S3 method for class 'lbi_test_obj'
autoplot(object, ...)

Arguments

Value

A ggplot

See Also

lbi_test

Description

Plots the CUSUM detector path against its boundary, with one panel for each series and vertical markers at the end of the training sample and at the alarm date.

Usage

## S3 method for class 'monitor_cusum_obj'
autoplot(object, ...)

Arguments

Value

A ggplot

See Also

monitor_cusum

Description

Plots the LBI CUSUM detector path against its boundary, with one panel for each series and a vertical marker at the alarm date.

Usage

## S3 method for class 'monitor_lbi_obj'
autoplot(object, ...)

Arguments

Value

A ggplot

See Also

monitor_lbi

Description

Plots the monitoring statistic against its boundary, with one panel for each series and vertical markers at the end of the training sample and at the alarm date.

Usage

## S3 method for class 'monitor_obj'
autoplot(object, ...)

Arguments

Value

A ggplot

See Also

monitor

Description

Plots the quantile monitoring statistic against its boundary, with one panel for each series and a vertical marker at the alarm date.

Usage

## S3 method for class 'monitor_quantile_obj'
autoplot(object, ...)

Arguments

Value

A ggplot

See Also

monitor_quantile

Description

Bar chart of the quantile-DF statistic for each series against its critical value. Series that exceed it are flagged as detected.

Usage

## S3 method for class 'quantile_test_obj'
autoplot(object, ...)

Arguments

Value

A ggplot

See Also

quantile_test

Description

Takes a radf_distr object and returns a ggplot2 object.

Usage

## S3 method for class 'radf_distr'
autoplot(object, ...)

Arguments

Value

A ggplot2::ggplot()

Examples

# Monte Carlo distribution
mdist <- radf_mc_distr(n = 100, nrep = 1000)
autoplot(mdist)

# Wild bootstrap distribution (one facet per series)
wdist <- radf_wb_distr(sim_data)
autoplot(wdist)

# Panel sieve bootstrap distribution
sdist <- radf_sb_distr(sim_data, nboot = 500)
autoplot(sdist)

Description

autoplot.radf_obj takes a radf_obj and a radf_cv and returns a faceted ggplot object. shade is used as an input to shade_opt. It modifies the geom_rect layer that marks the exuberance periods.

Usage

## S3 method for class 'radf_obj'
autoplot(
  object,
  cv = NULL,
  sig_lvl = 95,
  option = c("gsadf", "sadf"),
  min_duration = 0L,
  select_series = NULL,
  nonrejected = FALSE,
  shade_opt = shade(),
  trunc = TRUE,
  include_negative = "DEPRECATED",
  ...
)

## S3 method for class 'radf_obj'
autoplot2(
  object,
  cv = NULL,
  sig_lvl = 95,
  option = c("gsadf", "sadf"),
  min_duration = 0L,
  select_series = NULL,
  nonrejected = FALSE,
  trunc = TRUE,
  shade_opt = shade(),
  ...
)

shade(
  fill = "grey55",
  fill_negative = fill,
  fill_ongoing = NULL,
  opacity = 0.3,
  ...
)

Arguments

Value

A ggplot

Examples

rsim_data <- radf(sim_data_wdate)

autoplot(rsim_data)

# Modify facet_wrap options through ellipsis
autoplot(rsim_data, scales = "free_y", dir = "v")

# Modify the shading options
autoplot(rsim_data, shade_opt = shade(fill = "pink", opacity = 0.5))

# Allow for nonrejected series to be plotted
autoplot(rsim_data, nonrejected = TRUE)

# Remove the shading completely (2 ways)
autoplot(rsim_data, shade_opt = NULL)
autoplot(rsim_data, shade_opt = shade(opacity = 0))

# Plot only the series with the shading options
autoplot2(rsim_data)
autoplot2(rsim_data, trunc = FALSE) # keep the minw period

# We will need ggplot2 from here on out
library(ggplot2)

# Change (overwrite) color, size or linetype
autoplot(rsim_data) +
  scale_color_manual(values = c("black", "black")) +
  scale_linewidth_manual(values = c(0.9, 1)) +
  scale_linetype_manual(values = c("solid", "solid"))

# Change names through labeller (first way)
custom_labels <- c("psy1" = "new_name_for_psy1", "psy2" = "new_name_for_psy2")
autoplot(rsim_data, labeller = labeller(.default = label_value, id = as_labeller(custom_labels)))

# Change names through labeller (second way)
custom_labels2 <- series_names(rsim_data)
names(custom_labels2) <- custom_labels2
custom_labels2[c(3, 5)] <- c("Evans", "Blanchard")
autoplot(rsim_data, labeller = labeller(id = custom_labels2))

# Or change names before plotting
series_names(rsim_data) <- LETTERS[1:5]
autoplot(rsim_data)

# Change Theme options
autoplot(rsim_data) +
  theme(legend.position = "right")

Description

Plots each series with vertical markers at the estimated collapse and recovery dates.

Usage

## S3 method for class 'radf_recovery_obj'
autoplot(object, ...)

Arguments

Value

A ggplot

See Also

radf_recovery

Description

Plots the supDF, supBZ and union statistics of each series against their bootstrap critical values at the chosen significance level.

Usage

## S3 method for class 'radf_sbz_union_obj'
autoplot(object, sig_lvl = 95, ...)

Arguments

Value

A ggplot

See Also

radf_sbz_union

Description

Plots the estimated root and its confidence interval for every episode, one panel per series.

Usage

## S3 method for class 'rootstamp_episodes'
autoplot(object, ...)

Arguments

Value

A ggplot

See Also

rootstamp

Description

Plots the sub-sample against the fitted explosive path implied by the estimated root, $y_1 \rho^{t-1}$.

Usage

## S3 method for class 'rootstamp_est'
autoplot(object, ...)

Arguments

Value

A ggplot

See Also

rootstamp

Description

Plots the recursive SSU statistic path against its critical value, with one panel for each series.

Usage

## S3 method for class 'ssu_test_obj'
autoplot(object, ...)

Arguments

Value

A ggplot

See Also

ssu_test

Description

cobubble_test tests whether two series that each contain an explosive episode are co-explosive. That is, it tests whether the linear combination y_t - alpha - beta * x_{t-lag} is stationary, so that the explosive dynamics in y and x are the same underlying phenomenon, possibly migrating from one series to the other with a lead or lag, and not independent explosive episodes.

Usage

cobubble_test(
  y,
  x,
  lag = NULL,
  lag_grid = -6:6,
  nboot = 499L,
  sig_lvl = 95,
  seed = NULL
)

Arguments

Details

radf is a right-tailed ADF-family test for the presence of explosiveness. This function is a stationarity (KPSS-type) test instead. The null hypothesis is co-explosivity, that is, that the residuals of y regressed on a constant and x_{t-lag} are I(0). The null limiting distribution of the statistic depends on the pattern of heteroskedasticity in the errors (Evripidou, Harvey, Leybourne & Sollis 2022, Theorem 1). The critical values therefore come from a wild bootstrap that reproduces this pattern of heteroskedasticity in the bootstrap samples (Theorem 2).

Value

An object of class cobubble_test_obj: a list with the observed statistic S, the (given or estimated) lag, the bootstrap critical value cv at sig_lvl, the bootstrap p-value p_value and reject, which is TRUE if S exceeds cv, that is, if co-explosivity is rejected.

Status

Note

The critical value comes from a wild bootstrap of the residuals (Theorem 2), which is computed internally in each call. There is no separate, reusable cv function for this test.

The function returns its own class and not radf_obj, so it does not work with summary(), ⁠\link{datestamp}⁠ and tidy. It has its own print() and autoplot() methods instead. print() shows the statistic, the critical value and the p-value. See vignette("naming-and-analysis", package = "exuber") for which functions fit the shared pipeline and which do not.

References

Evripidou, A. C., Harvey, D. I., Leybourne, S. J., & Sollis, R. (2022). Testing for co-explosive behaviour in financial time series. Oxford Bulletin of Economics and Statistics, 84(3), 624-650.

See Also

Other multivariate: contagion_reg(), radf_common()

Examples

# A co-explosive pair (the process of Evripidou et al.), which is not rejected
xy <- sim_coexplosive(n = 100, seed = 123)
res <- cobubble_test(xy$y, xy$x, nboot = 199L, seed = 1)
print(res)

# Force a specific lead or lag instead of estimating it
res_lag0 <- cobubble_test(xy$y, xy$x, lag = 0L, nboot = 199L, seed = 1)
print(res_lag0)

# Two independent bubbles: co-explosivity correctly rejected
cobubble_test(sim_data$psy1, sim_data$psy2, nboot = 199L, seed = 1)

# Plot the two series that are tested for co-explosivity
autoplot(res)

Description

contagion_reg estimates the time-varying contagion coefficient of Greenaway-McGrevy & Phillips (2016). It computes a fixed-window rolling AR(1) coefficient sequence for a "core" series and for a "satellite" series y and relates them with a functional (Nadaraya-Watson kernel) regression at a chosen delay d. The coefficient shows how strongly, and how it varies over time, the local persistence of the core series transmits to y, d periods later.

Usage

contagion_reg(
  y,
  core,
  S = NULL,
  d = 0L,
  h = NULL,
  r_grid = seq(0, 1, length.out = 100)
)

Arguments

Details

This is a minimal subset of the procedure in the paper. It contains the fixed-window AR(1) coefficient sequence (their eq. 1), the Nadaraya-Watson regression at a single supplied d (eq. 6) and leave-one-out cross-validated bandwidth selection (eq. 7). Their eq. 8, the automatic search over d, is not implemented. If you need a search, call contagion_reg once for each candidate d and compare the fit.

The paper performs no formal inference on the contagion coefficient itself, with no confidence bands and no hypothesis test. As in the paper, this function is a tool for point estimation and visualization and not a test.

Value

An object of class contagion_reg_obj: a list with the fixed-window AR(1) coefficient sequences (beta_core and beta_j), the selected or supplied bandwidth (h) and the estimated time-varying contagion coefficient (delta2, aligned with r_grid).

Status

Note

The function is not a hypothesis test. It performs no formal inference on the contagion coefficient, with no confidence bands and no significance test, so there is no critical value for it.

The function returns its own class and not radf_obj, so it does not work with summary(), ⁠\link{datestamp}⁠ and tidy. It has its own print() and autoplot() methods instead. print() shows the coefficient path, because the function performs no formal inference. See vignette("naming-and-analysis", package = "exuber") for which functions fit the shared pipeline and which do not.

References

Greenaway-McGrevy, R., & Phillips, P. C. B. (2016). Hot property in New Zealand: Empirical evidence of housing bubbles in the metropolitan centres. New Zealand Economic Papers, 50(1), 88-113.

See Also

cobubble_test for a different, symmetric bivariate bubble relationship that uses a hypothesis test.

Other multivariate: cobubble_test(), radf_common()

Examples

# A co-explosive pair: the AR coefficient of y follows that of x almost one for one
xy <- sim_coexplosive(n = 100, seed = 123)
res <- contagion_reg(xy$y, xy$x, d = 0L)
print(res)

# Plot the estimated time-varying contagion coefficient
autoplot(res)

# Compare a one-period lead (d = 1) with the contemporaneous case
res_d1 <- contagion_reg(xy$y, xy$x, d = 1L)
autoplot(res) +
  ggplot2::geom_line(data = data.frame(r = res_d1$r_grid, delta2 = res_d1$delta2),
    ggplot2::aes(r, delta2), color = "red", inherit.aes = FALSE)

Description

cusum_test implements the retrospective CUSUM ("cs"), generalized CUSUM ("gcs"), CUSUM-of-squares ("cssq") and generalized CUSUM-of-squares ("gcssq") tests of Kurozumi & Nishi (2025). These are the parameter-constancy statistics of Brown et al. (1975), applied to the first differences. The generalized versions take the supremum over every window start as well as every end point.

Usage

cusum_test(data, sig_lvl = 95, type = c("cs", "gcs", "cssq", "gcssq"))

Arguments

Details

CS and GCS reject when the cumulated increments become too large (right tail). CSSQ and GCSSQ are two-sided. They reject when the cumulated squared increments drift too far above or below their full-sample average, with half the level in each tail. The paper finds that the CUSUM-type tests lose almost all their power once the explosive coefficient is stochastic, while the CUSUM-SQ type keeps it. See ssu_test for the more powerful statistics of the paper.

Value

An object of class cusum_test_obj: a list with the statistic path stat (one value for each end point, and for the generalized versions the sup over window starts), stat_inf (the inf path, for CSSQ and GCSSQ only), the statistic sup (and inf), the critical values crit and detected.

Status

Note

All critical values are published asymptotic values (Table I of Kurozumi & Nishi 2025). No minimum window is needed, because the paper finds the statistics insensitive to it.

References

Kurozumi, E., & Nishi, M. (2025). Bubble testing with stochastically varying explosive coefficient. Journal of Time Series Analysis, 46(5), 945-965.

Brown, R. L., Durbin, J., & Evans, J. M. (1975). Techniques for testing the constancy of regression relationships over time. Journal of the Royal Statistical Society B, 37(2), 149-192.

See Also

ssu_test, and monitor_cusum for real-time CUSUM monitoring.

Other volatility-robust tests: radf_kp(), radf_sbz(), radf_sbz_union(), radf_sign(), radf_sign_dm(), radf_tt(), ssu_test()

Examples

y <- sim_psy1(n = 150, te = 75, tf = 150, c = 3, alpha = 1, seed = 2001,
  coef_noise = rnorm(149), coef_a = 4)
cusum_test(y, type = "cssq")
autoplot(cusum_test(y, type = "gcssq"))

Description

Computes the origination, termination and duration of the episodes in which a time series shows explosive dynamics.

Usage

datestamp(object, cv = NULL, min_duration = 0L, ...)

## S3 method for class 'radf_obj'
datestamp(
  object,
  cv = NULL,
  min_duration = 0L,
  sig_lvl = 95,
  option = c("gsadf", "sadf", "svadf"),
  nonrejected = FALSE,
  ...
)

Arguments

Details

datestamp also stores a vector that takes the value 1 when there is a period of explosive behavior and 0 otherwise. You can use it as a dummy variable for the occurrence of exuberance.

Value

A table with the following columns:

• Start:

• Peak:

• End:

• Duration:

• Signal:

• Ongoing:

A list with the estimated origination and termination dates of the episodes of explosive behavior and their duration.

Caveats

option = "svadf": Sarkar & Wells (2026) is a preprint that has not been peer reviewed, which is a weaker standard of evidence than for every other source this package implements. The function emits the same note as a message when you call it with this option. It detects at most one origination and collapse pair per series, as the procedure of the paper does, and it does not find every recurring episode as "gsadf" and "sadf" do.

References

Phillips, P. C. B., Shi, S., & Yu, J. (2015). Testing for Multiple Bubbles: Historical Episodes of Exuberance and Collapse in the S&P 500. International Economic Review, 56(4), 1043-1078. doi:10.1111/iere.12132

Sarkar, A., & Wells, M. T. (2026). Is there an AI bubble? Robust date-stamping for periods of exuberance. arXiv:2604.12062.

Examples

rsim_data <- radf(sim_data)

# SV-ADF asymmetric-threshold dating (no critical values needed)
datestamp(rsim_data, option = "svadf")


# The default `cv` is fetched from the shared critical-value store
# (network on first use); pass `cv = radf_mc_cv(nrow(sim_data))` to stay offline
ds_data <- datestamp(rsim_data)
ds_data

# Choose minimum window
datestamp(rsim_data, min_duration = psy_ds(nrow(sim_data)))

autoplot(ds_data)

Description

dating_hls dates a single bubble episode by fitting four candidate regime-dummy regressions of Delta y_t on y_{t-1} (unit-root-to-end, unit-root-bubble-unit-root, unit-root-bubble-collapse and unit-root-bubble-collapse-unit-root). It fits each by minimizing the residual sum of squares over candidate break fractions and selects among them by BIC.

Usage

dating_hls(data, trim = 0.05)

Arguments

Details

datestamp uses threshold crossing on the recursive BSADF statistic, and dating_pdc uses a fixed structure of three or four regimes with breaks that are estimated sequentially and not jointly. This function searches for the breakpoints jointly within each of four candidate regime structures and lets the BIC pick the structure itself. It can therefore distinguish a bubble that collapses to a new stationary regime (Model 3) from a bubble that fully reverts to a unit root (Model 4) and from a bubble that is still ongoing at the end of the sample (Model 1), which the fixed regime count of dating_pdc cannot. The cost is a joint grid search, in place of the sequential scan of dating_pdc that finds one break at a time.

Value

An object of class dating_hls_obj: a list with the selected model (model, one of 1:4), its breakpoint date or dates (origination, collapse and recovery, with NA for the breakpoints that the selected model does not have), and the BIC value of every candidate model (bic, which shows how close the selection was).

Status

Note

This is an SSR/BIC model-selection dating procedure and not a hypothesis test, so it needs no critical values.

The function returns its own class and not radf_obj, so it does not work with summary(), ⁠\link{datestamp}⁠ and tidy. It has its own print() and autoplot() methods instead. print() shows the dating table (model, origination, collapse, recovery). See vignette("naming-and-analysis", package = "exuber") for which functions fit the shared pipeline and which do not.

References

Harvey, D. I., Leybourne, S. J., & Sollis, R. (2017). Improving the accuracy of asset price bubble start and end date estimators. Journal of Empirical Finance, 40, 121-138.

See Also

dating_pdc for the cheaper sequential-splitting alternative that this function complements, and datestamp for the original threshold-crossing rule of PSY.

Other dating: dating_hlw(), dating_knp(), dating_pdc(), radf_recovery(), rootstamp()

Examples

res <- dating_hls(sim_data$psy1, trim = 0.05)
print(res)

# Plot the series with the selected model's breakpoint(s) overlaid
autoplot(res)

# A whole panel at once, with one subplot for each series
autoplot(dating_hls(sim_data, trim = 0.05))

Description

dating_hlw extends dating_hls to series with more than one explosive episode. It first runs the existing detection and dating of PSY (radf and datestamp) to locate a preliminary start and end for each episode. It then splits the sample into disjoint date windows around them and dates each window again with the SSR/BIC fitting of dating_hls, restricted to Models 2 and 4 for every window except the last.

Usage

dating_hlw(
  data,
  cv = NULL,
  minw = NULL,
  trim = 0.1,
  min_duration = NULL,
  nboot = 199L,
  seed = NULL,
  join = 3L
)

Arguments

Details

When exactly one episode is detected, the function reduces to dating_hls applied to the whole series, which is a stated property of the paper. The single window then runs over [1, n] and fits all four models.

Value

An object of class dating_hlw_obj: a list with one element for each series. Each element is a data frame with one row for each detected episode (model, origination, collapse, recovery). A series with no step-1 detected episode gets a data frame with zero rows.

Status

Note

The step-2 SSR/BIC dating within each window needs no critical values, as in dating_hls. The step-1 PSY detection and dating pass does use a wild bootstrap critical value (cv, nboot and seed below), but only to locate the preliminary episode windows and not for the dating step.

The function returns its own class and not radf_obj, so it does not work with summary(), ⁠\link{datestamp}⁠ and tidy. It has its own print() and autoplot() methods instead. print() shows the dating table (model, origination, collapse, recovery). See vignette("naming-and-analysis", package = "exuber") for which functions fit the shared pipeline and which do not.

References

Harvey, D. I., Leybourne, S. J., & Whitehouse, E. J. (2020). Date-stamping multiple bubble regimes. Journal of Empirical Finance, 58, 226-246.

See Also

dating_hls for the single-bubble fitting that this function wraps, and datestamp for the multi-bubble threshold-crossing dating of PSY.

Other dating: dating_hls(), dating_knp(), dating_pdc(), radf_recovery(), rootstamp()

Examples

res <- dating_hlw(sim_data$psy1, trim = 0.1, nboot = 199L, seed = 1)
print(res)

# Plot the breakpoints of every detected episode over the series
autoplot(res)

# A two-bubble series, for which dating_hls() alone would fit only one bubble
res2 <- dating_hlw(sim_psy2(n = 200, seed = 123), trim = 0.1, nboot = 199L, seed = 1)
autoplot(res2)

Description

dating_knp dates bubble episodes (origination and collapse) by minimizing a sum of squared residuals that is corrected by omitting a residual. The model has alternating regimes: a unit root, an explosive regime, and a unit root that resumes from a shifted level after an instantaneous collapse, and so on. Plain OLS over this model is provably inconsistent, because the estimate of the origination date converges to the true collapse date and not to the origination date. The default omit = TRUE fixes this by dropping the squared residual at each candidate collapse date from the objective before minimizing.

Usage

dating_knp(data, trim = 0.05, omit = TRUE, breaks = 2L)

Arguments

Details

breaks = 2 (the default) is the single-bubble model. More breaks use the dynamic-programming algorithm of Kejriwal, Nguyen & Perron, which returns the exact global minimizer of the objective in O(breaks * n^2). As in the paper, the number of breaks is taken as given, for example two for each episode that datestamp finds.

Value

An object of class dating_knp_obj: a list with origination and collapse (dates) and delta (the fitted explosive AR coefficient). For a single bubble these are named vectors with one value for each series, and for more bubbles they are matrices with one row for each bubble and one column for each series.

Status

Note

This is a model-selection dating procedure based on the residual sum of squares and not a hypothesis test, so it needs no critical values.

The function returns its own class and not radf_obj, so it does not work with summary(), ⁠\link{datestamp}⁠ and tidy. It has its own print() and autoplot() methods instead. print() shows the dating table (model, origination, collapse, recovery). See vignette("naming-and-analysis", package = "exuber") for which functions fit the shared pipeline and which do not.

References

Kejriwal, M., Nguyen, L., & Perron, P. (2025). An improved procedure for retrospectively dating the emergence and collapse of bubbles. Journal of Time Series Analysis, 46(5), 867-883.

See Also

dating_hls and dating_pdc for related SSR-based dating approaches.

Other dating: dating_hls(), dating_hlw(), dating_pdc(), radf_recovery(), rootstamp()

Examples

res <- dating_knp(sim_data$psy1, trim = 0.05)
print(res)
autoplot(res)

# Compare the bias-corrected estimate with the plain (inconsistent) OLS one,
# adding an extra reference line to the output of autoplot()
res_plain <- dating_knp(sim_data$psy1, trim = 0.05, omit = FALSE)
autoplot(res) +
  ggplot2::geom_vline(xintercept = as.numeric(res_plain$origination), linetype = 3)

# Two bubbles
dating_knp(sim_data$psy2, breaks = 4)

Description

dating_pdc dates a single bubble episode with the sequential sample-splitting method of Pang, Du & Chong (2021) and its four-regime extension by Kurozumi & Skrobotov (2023). The regime structure is fixed: a unit root, an explosive regime, a stationary collapse and, optionally, a final unit-root recovery regime. The breakpoints are estimated one at a time. Each is a closed-form residual-sum-of-squares minimization over a no-intercept AR(1) model, computed in $O(T)$ time with cumulative sums.

Usage

dating_pdc(
  data,
  regimes = 3L,
  trim = 0.05,
  type = c("ols", "wls"),
  kernel = c("gaussian", "uniform"),
  h = NULL
)

Arguments

Details

datestamp finds where the recursive BSADF statistic crosses a critical value. This function instead fits an explicit regime-switching model directly to the series, and it needs no critical values. PDC prove that the collapse date is identified first, because its effect on the residual sum of squares dominates that of the origination date. This justifies estimating the breaks sequentially and not jointly. The alternative of Harvey, Leybourne & Sollis (2017) selects among models by BIC and fits them jointly, and it is not implemented here (see the package's enhancement notes for the cost and benefit considerations).

type = "wls" adds the correction of Kurozumi & Skrobotov (2023) for time-varying volatility. The function first fits the plain ("ols") model and collects its fitted residuals for each regime. It then smooths their squares nonparametrically, with the same Nadaraya-Watson kernel and leave-one-out bandwidth estimator that exuber already uses in radf_sbz_cv and radf_kp. Finally it reruns the same sequential break search with each squared term weighted by the inverse of the estimated spot variance. This needs no new critical-value theory. Like the OLS version, it is a point estimate and not a threshold-crossing test.

Value

An object of class dating_pdc_obj. It is a data.frame with one row for each series and the columns origination, collapse and, if regimes = 4, recovery, which give the estimated break dates (or observation indices, if no date index is available). It has its own print() and autoplot() methods.

Status

Note

The function returns its own class and not radf_obj, so it does not work with summary(), ⁠\link{datestamp}⁠ and tidy. It has its own print() and autoplot() methods instead. print() shows the dating table (model, origination, collapse, recovery). See vignette("naming-and-analysis", package = "exuber") for which functions fit the shared pipeline and which do not.

References

Pang, T., Du, L., & Chong, T. T. L. (2021). Estimating multiple breaks in the bubble regime with SSR minimization. Journal of Management Science and Engineering.

Kurozumi, E., & Skrobotov, A. (2023). Bubble dating: a sequential testing approach.

Kurozumi, E., & Skrobotov, A. (2023). Improving the accuracy of bubble date estimators under time-varying volatility. arXiv:2306.02977.

See Also

datestamp for the threshold-crossing alternative of PSY.

Other dating: dating_hls(), dating_hlw(), dating_knp(), radf_recovery(), rootstamp()

Examples

# The unit-root, explosive, collapse and recovery process of sim_ps1() is the
# regime structure that dating_pdc() fits (true breaks at 40, 60 and 70)
y <- sim_ps1(n = 100, seed = 1)
res <- dating_pdc(y, regimes = 3L, trim = 0.05)
print(res)
autoplot(res)

# Four-regime extension, which adds a post-collapse recovery breakpoint
res4 <- dating_pdc(y, regimes = 4L, trim = 0.05)
autoplot(res4)

# Volatility-weighted (WLS) variant, robust to time-varying volatility
y_vol <- sim_ps1(n = 100, seed = 1, e = sim_vol_break(99))
autoplot(dating_pdc(y_vol, type = "wls"))

Description

Reports whether the null hypothesis of a unit root is rejected against the alternative of explosive behavior, for each series in a dataset.

Usage

diagnostics(object, cv = NULL, ...)

## S3 method for class 'radf_obj'
diagnostics(object, cv = NULL, option = c("gsadf", "sadf"), sig_lvl = 95, ...)

Arguments

Details

diagnostics also stores a vector that takes the value 1 when there is a period of explosive behavior and 0 otherwise.

Value

A list with the series that reject the null hypothesis (positive) and the series that do not (negative), together with the significance level.

Examples

# The default `cv` is fetched from the shared critical-value store
# (network on first use); pass `cv = radf_mc_cv(nrow(sim_data))` to stay offline
rsim_data <- radf(sim_data)
diagnostics(rsim_data)

diagnostics(rsim_data, option = "sadf")

# Gate on the 90% critical value instead of the 95% default
diagnostics(rsim_data, sig_lvl = 90)

Description

The naming conventions (radf_, _test, dating_ and monitor_, see vignette("naming-and-analysis")) are a guide and not a contract. They are easy to misremember, and they are occasionally traded off deliberately: monitor() uses the ADF family internally but is named after what it does. This function is the queryable source of truth for which of the package's test, dating, monitoring and root-inference functions belong to which family.

Usage

exuber_functions(family = NULL)

Arguments

Value

A tibble with the columns name, family and description, and one row for each function.

See Also

vignette("naming-and-analysis", package = "exuber") for the full naming scheme and for which functions work with summary(), datestamp, tidy and autoplot.

Examples

exuber_functions()
exuber_functions(family = "monitor")
exuber_functions(family = "test")

Description

Retrieve or replace the index of an object.

Usage

index(x, ...)

index(x) <- value

Arguments

Details

If you do not specify an index for the estimation, the function generates a pseudo-index, a sequential numeric series. After the estimation you can use index() to retrieve the index or ⁠index<-()⁠ to replace it. The index can be numeric or Date.

Value

The index of the object (a Date vector when the input data carried one, otherwise an integer sequence). The replacement form returns the modified object.

Examples

# A plain numeric vector gets a pseudo-index (1, 2, 3, ...)
rsim <- radf(sim_data)
head(index(rsim))

# A data.frame with a Date column uses it as the index automatically
rsim_wdate <- radf(sim_data_wdate)
head(index(rsim_wdate))
class(index(rsim_wdate))

# autoplot() uses index() internally for the x-axis
autoplot(rsim_wdate)

# Replace the index, e.g. with a custom Date sequence
index(rsim) <- seq(as.Date("2000-01-01"), by = "month", length.out = length(index(rsim)))
head(index(rsim))
autoplot(rsim)

Description

lbi_test implements the static locally best invariant (LBI) test of Breitung & Diegel (2025) for a bubble that is known, or assumed, to span the entire sample: LBI = (y_T - y_1) / (sigma_tilde * sqrt(T - 1)), where sigma_tilde^2 is the sample variance of the first differences. The test is robust to heteroskedasticity by construction, because the invariance property of the statistic does not depend on the exact form of the innovation variance. The null distribution is standard normal, so the test needs no bootstrap, no simulation and no published table.

Usage

lbi_test(data, sig_lvl = 95)

Arguments

Details

Only the static test, which uses a single full-sample window, is implemented. The main contribution of Breitung & Diegel is a sequential, exponentially weighted extension for monitoring a series when the start date is unknown. Its exact weighting scheme and boundary constant are not pinned down here, so this function does not implement it.

Value

An object of class lbi_test_obj: a list with the test statistic stat, the standard-normal critical value crit and detected (logical, stat > crit).

Status

Note

The critical value is closed-form: it is the standard normal (qnorm) quantile at sig_lvl, so no bootstrap, simulation or table is needed.

The function returns its own class and not radf_obj, so it does not work with summary(), ⁠\link{datestamp}⁠ and tidy. It has its own print() and autoplot() methods instead. print() shows the statistic, the critical value and whether the bubble was detected. See vignette("naming-and-analysis", package = "exuber") for which functions fit the shared pipeline and which do not.

References

Breitung, J., & Diegel, M. (2025). A locally best invariant sequential test for explosive behavior in the presence of nonstationary volatility. Journal of Time Series Analysis.

See Also

radf for the recursive ADF-family alternative that this test complements.

Other alternative tests: quantile_test()

Examples

y <- sim_psy1(n = 60, te = 1, tf = 60, seed = 1) # explosive from the start
res <- lbi_test(y)
print(res)

# Compare the statistic to its critical value
autoplot(res)

Description

monitor implements real-time monitoring. You fix a training window [1, T*] that is assumed free of exuberance and calibrate a critical value on it. The function then compares the running recursive statistic at each subsequent point T*+1, ..., T with that fixed boundary and flags the first date at which the boundary is breached.

Usage

monitor(
  data,
  r_star = 0.5,
  minw = NULL,
  nboot = 500L,
  sig_lvl = 95,
  lag = 0L,
  type = c("fixed", "aic", "bic"),
  seed = NULL,
  boundary = c("bootstrap", "kurozumi", "fluc"),
  s0 = 0
)

Arguments

Details

boundary = "bootstrap" (the default) implements Phillips & Shi (2020). The boundary is a wild-bootstrap quantile of the GSADF-type statistic (see the tb parameter of radf_wb_ps_cv), and it is compared with the bsadf sequence of radf(). The function calibrates on the training window only (data[1:T*]) and not on the full series. The null-model fit inside radf_wb_ps_cv (adf_res()) uses all the data it is given and does not truncate them to tb, so passing data after T*, which may be explosive, directly to it would leak future information into the calibration of the null.

boundary = "kurozumi" implements the closed-form alternative of Kurozumi (2020). It needs no bootstrap and compares a published constant (his Table 1) with the badf sequence of radf(). The default s0 = 0 gives his SADF(k) detector, where the window start is fixed at 1. Setting s0 to 0.4 or 0.8 switches to his GSADF_{s0}(k) generalization. The window start then ranges over [1, floor(T* * s0)] and is not fixed at 1, and the comparison uses his boundary function, which varies with k and is not constant, together with its own published scaling constant. sig_lvl must be one of 90, 95 or 99, the levels that his table tabulates.

boundary = "fluc" implements the FLUC detector of Homm & Breitung (2012). Their DF_{t/n} is also exactly the badf sequence of radf(), and it is compared with a published constant from their Table 7 (the case without detrending) and not with a simulated one. sig_lvl must be one of 90, 95 or 99.

Value

An object of class monitor_obj: a list with the full-sample statistic path (stat, which is bsadf for boundary = "bootstrap" and badf for "kurozumi" and "fluc"), the calibrated boundary (one flat value for each series), the length of the training window T_star, and alarm and alarm_date (the first observation or date in the monitoring period at which stat breaches the boundary, NA if it never does).

Status

Note

The function returns its own class and not radf_obj, so it does not work with summary(), ⁠\link{datestamp}⁠ and tidy. It has its own print() and autoplot() methods instead. print() shows the boundary and the alarm. See vignette("naming-and-analysis", package = "exuber") for which functions fit the shared pipeline and which do not.

References

Phillips, P. C., & Shi, S. (2020). Real time monitoring of asset markets: Bubbles and crises. In Handbook of Statistics (Vol. 42, pp. 61-80). Elsevier.

Kurozumi, E. (2020). Asymptotic properties of bubble monitoring tests. Econometric Reviews, 39(5), 510-538.

Homm, U., & Breitung, J. (2012). Testing for speculative bubbles in stock markets: A comparison of alternative methods. Journal of Financial Econometrics, 10(1), 198-231.

See Also

radf_wb_ps_cv for the underlying wild bootstrap, and datestamp for the existing full-sample dating of origination and collapse, which is not a monitoring procedure.

Other monitoring: monitor_cusum(), monitor_lbi(), monitor_quantile()

Examples

# A bubble-free training window (first half), explosive from t = 150 on
y <- sim_psy1(n = 200, te = 150, tf = 200, seed = 7)
# Default: Phillips & Shi (2020) wild bootstrap boundary
mon <- monitor(y, r_star = 0.5, nboot = 200)
print(mon)
autoplot(mon)

# Closed-form boundary of Kurozumi (2020), which needs no bootstrap
mon_kz <- monitor(y, r_star = 0.5, boundary = "kurozumi")
autoplot(mon_kz)

# Homm & Breitung (2012) FLUC boundary
autoplot(monitor(y, r_star = 0.5, boundary = "fluc"))

Description

monitor_cusum implements the CUSUM real-time monitoring procedure of Homm & Breitung (2012). You fix a training window [1, T*] that is assumed free of exuberance. The function then compares the standardized cumulative sum of the post-training first differences, S_t = (y_t - y_{T*}) / sigma_hat_t, with the closed-form boundary c_t * sqrt(t), where c_t = sqrt(b_alpha + log(t / T*)), and flags the first date at which the boundary is breached.

Usage

monitor_cusum(
  data,
  r_star = 0.5,
  b_alpha = 4.6,
  boundary = c("asymptotic", "finite"),
  sig_lvl = 95,
  type = c("standard", "kernel"),
  h = 20,
  kernel = c("gaussian", "uniform")
)

Arguments

Value

An object of class monitor_cusum_obj: a list with the statistic path in the monitoring region (S) and boundary, the length of the training window T_star, and alarm and alarm_date (the first breach, NA if there is none).

Status

Note

The boundary is closed-form throughout. It is either a fixed asymptotic constant (boundary = "asymptotic", b_alpha = 4.6) or a lookup in the published finite-sample table (boundary = "finite", Table 8 of Homm & Breitung 2012). No simulation and no separate cv function is needed.

monitor (Family A) uses a recursive ADF-family statistic whose boundary needs a wild bootstrap for calibration. This function uses a structurally different statistic, a standardized running sum and not a recursive regression. Its boundary is an asymptotic closed form (the inequality of Chu, Stinchcombe & White 1996, eq. 28 of HB), so it needs no bootstrap and no simulation and does not depend on the data beyond the running variance estimate itself.

type = "kernel" instead uses the volatility-robust modification ("CUSUMV") of Astill, Harvey, Leybourne, Taylor & Zu (2023). Before cumulating, it standardizes each first difference by its own one-sided kernel spot-variance estimate (their eq. 6-7) and not by a single running variance. Their Corollary 1 establishes that the same boundary function gives a controlled asymptotic false-alarm rate even under time-varying volatility. The standard CUSUM statistic does not have this property, because its size-control result requires homoskedasticity.

The function returns its own class and not radf_obj, so it does not work with summary(), ⁠\link{datestamp}⁠ and tidy. It has its own print() and autoplot() methods instead. print() shows the boundary and the alarm. See vignette("naming-and-analysis", package = "exuber") for which functions fit the shared pipeline and which do not.

References

Homm, U., & Breitung, J. (2012). Testing for speculative bubbles in stock markets: A comparison of alternative methods. Journal of Financial Econometrics, 10(1), 198-231.

Chu, C. S. J., Stinchcombe, M., & White, H. (1996). Monitoring structural change. Econometrica, 64(5), 1045-1065.

Astill, S., Harvey, D. I., Leybourne, S. J., Taylor, A. M. R., & Zu, Y. (2023). CUSUM-based monitoring for explosive episodes in financial data in the presence of time-varying volatility. Journal of Financial Econometrics, 21(1), 187-227.

See Also

monitor for the monitoring alternative based on the recursive ADF (Family A).

Other monitoring: monitor(), monitor_lbi(), monitor_quantile()

Examples

# A martingale training window, explosive from t = 150 to the sample end
y <- sim_psy1(n = 200, te = 150, tf = 200, seed = 7)
res <- monitor_cusum(y, r_star = 0.5)
print(res) # alarm should fire soon after t = 150
autoplot(res)

# Volatility-robust "CUSUMV" variant (Astill, Harvey, Leybourne, Taylor & Zu 2023):
# the same bubble, but volatility triples at t = 120, after the training window
y_vol <- sim_psy1(n = 200, te = 150, tf = 200, seed = 7,
  e = sim_vol_break(199, tau = 0.6))
res_kernel <- monitor_cusum(y_vol, r_star = 0.5, type = "kernel")
autoplot(res_kernel)

Description

monitor_lbi implements the sequential, constant-boundary extension of the locally best invariant statistic of lbi_test. It monitors a series in real time when the start date of the bubble is unknown. After a training window [1, T*] that is assumed free of exuberance, it compares the (optionally exponentially weighted) partial sum of the post-training first differences with a constant boundary and flags the first monitoring date at which the boundary is breached.

Usage

monitor_lbi(data, r_star = 0.5, c_bar = 0, sig_lvl = 95)

Arguments

Details

Their eq. 15 shows that this partial sum, normalized by the fixed length of the monitoring horizon, converges to a standard Brownian motion on [0, 1] under the null. The normalization is not sqrt(t), which differs from the Chu-Stinchcombe-White-style boundary of monitor_cusum. A single constant boundary therefore controls size uniformly across the whole monitoring window. The paper shows that this constant-boundary detector ("mCUSUM" at c_bar = 0, "wCUSUM" at c_bar > 0) is more powerful than the classical CUSUM test with a time-varying boundary, to which it is compared.

Value

An object of class monitor_lbi_obj: a list with the statistic path in the monitoring region (stat), the constant boundary, the length of the training window T_star, and alarm and alarm_date (the first breach, NA if there is none).

Status

Note

The critical value is the published constant boundary in Table 1 of Breitung & Diegel (2025), so it is a table lookup and needs no simulation.

The function returns its own class and not radf_obj, so it does not work with summary(), ⁠\link{datestamp}⁠ and tidy. It has its own print() and autoplot() methods instead. print() shows the boundary and the alarm. See vignette("naming-and-analysis", package = "exuber") for which functions fit the shared pipeline and which do not.

References

Breitung, J., & Diegel, M. (2025). A locally best invariant sequential test for explosive behavior in the presence of nonstationary volatility. Journal of Time Series Analysis.

See Also

lbi_test for the static version, which assumes a known bubble window that spans the full sample. monitor_cusum and monitor are monitoring detectors with a structurally different construction.

Other monitoring: monitor(), monitor_cusum(), monitor_quantile()

Examples

# A martingale training window, explosive from t = 150 to the sample end
y <- sim_psy1(n = 200, te = 150, tf = 200, seed = 7)
res <- monitor_lbi(y, r_star = 100)
print(res) # alarm should fire soon after t = 150
autoplot(res)

# wCUSUM: exponentially up-weight later monitoring observations
autoplot(monitor_lbi(y, r_star = 100, c_bar = 2))

Description

monitor_quantile implements the QPWY and QPSY real-time monitoring strategies of Wu, Shi & Wu (2025). They are quantile-regression (QR) analogues of the recursive ADF t-statistics of PWY and PSY, and they test at a chosen conditional quantile tau, in contrast to the single full-sample test in quantile_test. type = "qpwy" uses the expanding window [1, r], which has the shape of badf in radf. type = "qpsy" also takes the supremum over every window start, which has the shape of bsadf.

Usage

monitor_quantile(
  data,
  tau = 0.5,
  minw = NULL,
  nrep = 500L,
  sig_lvl = 95,
  seed = NULL,
  type = c("qpwy", "qpsy"),
  boundary = c("asymptotic", "bootstrap")
)

Arguments

Details

The point statistic needs genuine QR fits, because there is no closed-form recursive update as there is for OLS. QPWY needs O(T) fits and QPSY needs O(T^2), so QPSY takes seconds for each series at n = 200 and the time grows quadratically.

The boundary is a single flat value and not one value for each r. It is the quantile of the supremum of each null path, constructed like the sadf_cv of radf_mc_cv, which controls the first-crossing false-alarm rate. boundary chooses how the null paths are generated.

With boundary = "asymptotic" (the default) the function simulates the limiting null distribution delta * Q_{r1,r2} + sqrt(1 - delta^2) * Z_{r1,r2} in each call. Here delta is a correlation coefficient estimated from the data (as in quantile_test), Q is the Dickey-Fuller t functional and Z is its counterpart driven by an independent Brownian motion. Both are simulated for every window, so no QR fits are needed and the boundary is cheap.

With boundary = "bootstrap" the function applies Algorithm 1 of Wu, Shi & Wu (2025) to the whole path. It resamples the centred first differences of the series with replacement, cumulates them into a null random walk and recomputes the QPWY or QPSY path on it, nrep times. The boundary is the quantile of the nrep path maxima. It follows the finite-sample distribution of the statistic in the data at hand, so it removes most of the size distortion of the asymptotic boundary described below. Each replicate costs a full statistic path, which is O(T) QR fits for QPWY and O(T^2) for QPSY. The paper also discards the first 100 draws of each resample. We omit this step because the statistic has an intercept, so it does not depend on the level of the series, and the draws are i.i.d., so the retained values already have the same distribution. Set options(exuber.parallel = TRUE) to spread the replicates over workers.

Value

An object of class monitor_quantile_obj: a list with the statistic path stat, the flat boundary, the estimated delta (reported for both boundaries, used only by the asymptotic one), and alarm and alarm_date (the first breach, NA if there is none).

Caveats

The asymptotic boundary is well sized near the median in finite samples, with a false-alarm rate of 3.5 to 4.0\ innovations, tau = 0.5). Away from the median the small early windows make both statistics oversized, and QPSY badly so even with Gaussian innovations. The false-alarm rate of QPSY is 35\ innovations it is 21\ (n = 100). With $t_3$ innovations the rate for QPWY is 7.5 to 8.5\ tau = 0.2 and 0.8 and 12.5\ Wu, Shi & Wu advise against extreme quantiles in small samples and use bootstrap critical values for monitoring, which is what boundary = "bootstrap" provides. With it the false-alarm rate was between 4.5 and 6.5\ the cases above (n = 100 for QPWY, n = 60 for QPSY), including QPSY at tau = 0.8 and 0.9. The exception is QPSY with $t_3$ innovations at tau = 0.9, which stays at 10\ less extreme tau help. For type = "qpsy" with tau away from 0.5 and the asymptotic boundary, the function emits a short pointer as a message (see suppressMessages) and stores it as attr(x, "caveat"). The numbers for both boundaries are in docs/alternative-paradigms.md.

Status

Note

The function simulates the critical value (the boundary) internally in each call, with an unexported helper, quantile_boundary_sim. There is currently no reusable exported cv counterpart for this function. This gap is tracked separately and is not addressed here.

The function returns its own class and not radf_obj, so it does not work with summary(), ⁠\link{datestamp}⁠ and tidy. It has its own print() and autoplot() methods instead. print() shows the statistic, the boundary and delta. See vignette("naming-and-analysis", package = "exuber") for which functions fit the shared pipeline and which do not.

References

Wu, R., Shi, S., & Wu, J. (2025). Quantile analysis for financial bubble detection and surveillance. Journal of Time Series Analysis, 46(5), 908-931.

See Also

quantile_test for the static, full-sample version of this test, and monitor for the monitoring alternative based on OLS.

Other monitoring: monitor(), monitor_cusum(), monitor_lbi()

Examples

# Heavy-tailed (t3) innovations, explosive from t = 150 to the sample end
y <- sim_psy1(n = 200, te = 150, tf = 200, seed = 7,
  e = sim_innov(199, dist = "t", df = 3))
res <- monitor_quantile(y, tau = 0.5, nrep = 100, seed = 1)
print(res)
autoplot(res)

# Monitoring an upper quantile is typically more powerful for right-tailed
# bubbles, but see the Caveats section on extreme quantiles
autoplot(monitor_quantile(y, tau = 0.8, nrep = 100, seed = 1))

# QPSY: supremum over window starts too (O(n^2) QR fits, slower)
monitor_quantile(y[101:200], tau = 0.5, nrep = 100, seed = 1, type = "qpsy")

# Bootstrap boundary (Algorithm 1 of Wu, Shi & Wu): one full statistic path per
# replicate, so keep nrep small in a first run
monitor_quantile(y, tau = 0.8, nrep = 49, seed = 1, boundary = "bootstrap")

Description

This function finds the number of observations in the window over which size is to be controlled.

Usage

ps_tb(n, freq = c("monthly", "quarterly", "annual", "weekly"), size = 2)

Arguments

Value

A single integer, the length of the training window in observations, suitable for the r_star argument of monitor.

References

Phillips, P. C., & Shi, S. (2020). Real time monitoring of asset markets: Bubbles and crises. In Handbook of Statistics (Vol. 42, pp. 61-80). Elsevier.

Shi, S., Hurn, S., Phillips, P.C.B., 2018. Causal change detection in possibly integrated systems: Revisiting the money-income relationship.

Examples

# Training window that controls size over a 2-year span of monthly data
tb <- ps_tb(100, freq = "monthly", size = 2)
tb


# Use it directly as the training window of monitor()
monitor(sim_data, r_star = tb, boundary = "kurozumi")

Description

psy_minw and psy_ds use the rules of thumb proposed by Phillips et al. (2015) to compute the minimum window size and the minimum duration of an episode of exuberance, respectively.

Usage

psy_minw(n)

psy_ds(n, rule = 1, delta = 1)

Arguments

Details

For the minimum duration period, psy_ds allows the user to choose from two rules:

$\text{rule}_1 = \delta \log(n) \quad \text{and} \quad \text{rule}_2 = \delta \log(n)/n$

delta depends on the frequency of the data and the minimal duration condition.

Value

psy_minw: a single integer, the minimum window length. psy_ds: a single integer, the minimum episode duration in observations.

References

Phillips, P. C. B., Shi, S., & Yu, J. (2015). Testing for Multiple Bubbles: Historical Episodes of Exuberance and Collapse in the S&P 500. International Economic Review, 56(4), 1043-1078. doi:10.1111/iere.12132

Examples

psy_minw(100)
psy_ds(100)

Description

quantile_test implements the "global test" of Wu, Shi & Wu (2025). It is a quantile-regression (QR) analogue of the Dickey-Fuller t-ratio, and it tests for explosive behavior at a chosen conditional quantile tau of y_t given y_{t-1}, and not at the conditional mean. It is a single static test and not a recursive scan. The comparable statistic in radf is the single-shot adf statistic and not the recursive bsadf.

Usage

quantile_test(
  data,
  tau = "optimal",
  tau_grid = seq(0.2, 0.8, by = 0.05),
  nrep = 1000L,
  sig_lvl = 95,
  seed = NULL
)

Arguments

Details

tau = "optimal" (the default) selects the quantile that minimizes the asymptotic variance of the QR estimator (their eq. 33) by grid search over tau_grid. The grid excludes the extreme quantiles, which the paper itself recommends avoiding at practical sample sizes.

The function simulates the critical value in each call and does not use a fixed table. There is currently no reusable exported cv counterpart for this function. This gap is tracked separately and is not addressed here. The limiting null distribution of the statistic is sqrt(1 - delta^2) * z + delta * Q, where z ~ N(0, 1), delta is a correlation coefficient estimated from the data, and Q is the standard demeaned Dickey-Fuller t-statistic distribution. The function simulates Q with the same random-walk-plus-OLS t-statistic construction used elsewhere in this package (see radf_mc_cv).

Value

An object of class quantile_test_obj: a list with the test statistic tstat, the selected tau, the estimated correlation delta, the simulated crit value and detected (logical, tstat > crit).

Status

Note

The function returns its own class and not radf_obj, so it does not work with summary(), ⁠\link{datestamp}⁠ and tidy. It has its own print() and autoplot() methods instead. print() shows the statistic, the critical value and delta. See vignette("naming-and-analysis", package = "exuber") for which functions fit the shared pipeline and which do not.

References

Wu, R., Shi, S., & Wu, J. (2025). Quantile analysis for financial bubble detection and surveillance. Journal of Time Series Analysis, 46(5), 908-931.

See Also

radf for the family of tests based on mean regression (ADF, SADF and GSADF) that this test complements.

Other alternative tests: lbi_test()

Examples

# Heavy-tailed (t3) innovations, where a quantile test is more useful than a test of the mean
y <- sim_psy1(n = 100, seed = 1, e = sim_innov(99, dist = "t", df = 3))
res <- quantile_test(y, nrep = 100, seed = 1)
print(res)
autoplot(res)

# Test at a fixed upper quantile instead of the optimal one
autoplot(quantile_test(y, tau = 0.9, nrep = 100, seed = 1))

Description

radf computes the recursive univariate and panel augmented Dickey-Fuller test statistics (ADF, SADF, GSADF and their backward sequences).

Usage

radf(data, minw = NULL, lag = 0L)

Arguments

Details

radf() is vectorized: it handles several series at once. The panel statistics are computed as a by-product of the univariate estimations, at almost no additional cost.

Value

A list that contains the unit root test statistics (sequence):

And attributes:

References

Phillips, P. C. B., Wu, Y., & Yu, J. (2011). Explosive Behavior in The 1990s Nasdaq: When Did Exuberance Escalate Asset Values? International Economic Review, 52(1), 201-226.

Phillips, P. C. B., Shi, S., & Yu, J. (2015). Testing for Multiple Bubbles: Historical Episodes of Exuberance and Collapse in the S&P 500. International Economic Review, 56(4), 1043-1078. doi:10.1111/iere.12132

Pavlidis, E., Yusupova, A., Paya, I., Peel, D., Martínez-García, E., Mack, A., & Grossman, V. (2016). Episodes of exuberance in housing markets: in search of the smoking gun. The Journal of Real Estate Finance and Economics, 53(4), 419-449. doi:10.1007/s11146-015-9531-2

Examples

# We will use simulated data that are stored as data
sim_data

rsim <- radf(sim_data)

str(rsim)

# We would also use data that contain a Date column
sim_data_wdate

rsim_wdate <- radf(sim_data_wdate)

tidy(rsim_wdate)

augment(rsim_wdate)

tidy(rsim_wdate, panel = TRUE)

head(index(rsim_wdate))

# For lag = 1 and minimum window = 20
rsim_20 <- radf(sim_data, minw = 20, lag = 1)

Description

radf_common tests for a bubble that is common to a panel of series (Chen, Phillips & Shi, 2023). It extracts the first principal component of the panel and runs the ordinary radf test on it. The output is an ordinary radf_obj, so every downstream method (tidy(), autoplot(), datestamp() and so on) works on it without further effort.

Usage

radf_common(data, minw = NULL, r = 1)

Arguments

Details

Theorem 4.3 of the paper claims that the null limiting distribution of the resulting statistic is asymptotically identical to the standard PSY/GSADF distribution, which would let radf_mc_cv apply directly. An independent validation found that this identity does not hold at practical panel widths N. At N = 100 the true critical value is more than double that of radf_mc_cv, and the gap grows as N increases. PCA on a panel of independent (non-cointegrated) I(1) series does not behave like a single random walk once there are more series from which transient co-movement can arise. Use radf_common_cv for the critical values and not radf_mc_cv. The critical values of radf_mc_cv do not depend on the panel width, and they are badly undersized here once N grows past a handful of series.

Value

A radf_obj (see radf) computed on the first principal component of the panel, with the fitted prcomp object attached as an attribute (attr(x, "prcomp")).

Status

Note

The function returns the output of radf() computed on the extracted factor, and radf_common_cv() computes the full time-varying boundary alongside the scalar critical values. The full summary(), datestamp, tidy and autoplot pipeline therefore works (see vignette("naming-and-analysis", package = "exuber")).

References

Chen, Y., Phillips, P. C. B., & Shi, S. (2023). Common Bubble Detection in Large Dimensional Financial Systems. Journal of Financial Econometrics, 21(4), 989-1063.

See Also

radf for the underlying test, which is unmodified, and radf_common_cv for its critical values, which are specific to the panel width.

Other multivariate: cobubble_test(), contagion_reg()

Examples

# A panel of 5 series driven by one shared latent bubble factor
x <- sim_common(n_series = 5, n = 100, seed = 123)
res <- radf_common(x, minw = 20)
print(res)

# radf_common_cv() is needed here and radf_mc_cv() does not apply (see Details)
cv <- radf_common_cv(n = 100, N = ncol(x), minw = 20)
summary(res, cv = cv)

# The result is an ordinary radf_obj, so autoplot() and datestamp() work directly
autoplot(res, cv = cv)
datestamp(res, cv = cv)

Description

radf_common_cv simulates critical values for radf_common under its own null of no common explosive factor. The null is a panel of N independent random walks, from which one principal component is extracted and tested exactly as in radf_common. Independent validation showed that radf_mc_cv, whose critical values do not depend on the panel width, is badly undersized as a stand-in for this null once N grows past a handful of series. The null distribution here does depend on N, so N must match the panel on which radf_common was run.

Usage

radf_common_cv(n, N, minw = NULL, nrep = 1000L, seed = NULL)

Arguments

Value

A list with adf_cv, sadf_cv, gsadf_cv, badf_cv and bsadf_cv. It has the same shape as the return value of radf_mc_cv, so you can use it as the cv argument of datestamp, tidy and autoplot on a radf_common result.

Status

References

Chen, Y., Phillips, P. C. B., & Shi, S. (2023). Common Bubble Detection in Large Dimensional Financial Systems. Journal of Financial Econometrics, 21(4), 989-1063.

See Also

radf_common, radf_mc_cv

Other critical values: radf_mc_cv(), radf_recovery_cv(), radf_sb_cv(), radf_sbz_cv(), radf_sign_cv(), radf_sign_dm_cv(), radf_tt_cv(), radf_wb_cv(), radf_wb_ps_cv()

Examples

cv <- radf_common_cv(n = 100, N = 5, minw = 20, nrep = 200)
tidy(cv)

Description

radf_kp implements the heteroskedasticity-robust PSY test of Harvey, Leybourne, Taylor & Zu (2024), which needs no bootstrap. It "purges" unconditional heteroskedasticity by dividing each first difference of the series by a kernel spot-volatility estimate (eq. 4-5) and cumulating the result. It then runs the ordinary radf, with an intercept, on the purged series.

Usage

radf_kp(data, minw = NULL, kernel = c("gaussian", "uniform"), h = NULL)

Arguments

Details

The paper proves (Theorem 1 and Remark 3.2) that the null limiting distribution of the purged statistic is identical to the standard homoskedastic GSADF null. radf_mc_cv, the existing and already fast Monte Carlo critical values of exuber, therefore apply directly to the result. Unlike radf_wb_cv and radf_sbz_cv, no new bootstrap or simulation code is needed.

Only the with-intercept variant ($PSY_\sigma$ in the paper) is implemented. The paper also proposes a variant without an intercept and a union-of-rejections test that combines both. They are not implemented here (see the package's enhancement notes for the cost and benefit considerations).

Value

A radf_obj with the same structure as the output of radf, computed on the volatility-purged series, so radf_mc_cv, tidy() and the other methods apply directly.

Status

Note

The function returns the unmodified output of radf, so the full summary(), datestamp, tidy and autoplot pipeline works exactly as it does for plain radf() (see vignette("naming-and-analysis", package = "exuber")).

References

Harvey, D. I., Leybourne, S. J., Taylor, A. M. R., & Zu, Y. (2024). A new heteroskedasticity-robust test for explosive bubbles. Journal of Time Series Analysis. doi:10.1111/jtsa.12784

See Also

radf_mc_cv for the critical values of this test, which are unmodified, radf_wb_cv for a bootstrap-based alternative, and radf_tt for another alternative that needs no bootstrap.

Other volatility-robust tests: cusum_test(), radf_sbz(), radf_sbz_union(), radf_sign(), radf_sign_dm(), radf_tt(), ssu_test()

Examples

# Volatility triples half-way through the sample. This is the case of
# non-stationary volatility that this test is built for, and plain radf()
# over-rejects here
y <- sim_psy1(n = 200, seed = 1, e = sim_vol_break(199))
res <- radf_kp(y, minw = 20)
print(res)

# radf_mc_cv() applies unmodified (see Details)
cv <- radf_mc_cv(n = attr(res, "n"), minw = 20)
summary(res, cv = cv)
autoplot(res, cv = cv)

Description

radf_mc_cv computes Monte Carlo critical values for the recursive unit root tests. radf_mc_distr computes the simulated distribution.

Usage

radf_mc_cv(n, minw = NULL, nrep = 1000L, seed = NULL, lag = 0)

radf_mc_distr(n, minw = NULL, nrep = 1000L, seed = NULL, lag = 0)

Arguments

Value

For radf_mc_cv, a list with the critical values for the ADF, BADF, BSADF and GSADF test statistics. For radf_mc_distr, a list with the ADF, SADF and GSADF distributions.

See Also

radf_wb_cv for wild bootstrap critical values and radf_sb_cv for sieve bootstrap critical values

Other critical values: radf_common_cv(), radf_recovery_cv(), radf_sb_cv(), radf_sbz_cv(), radf_sign_cv(), radf_sign_dm_cv(), radf_tt_cv(), radf_wb_cv(), radf_wb_ps_cv()

Examples

# Default minimum window
mc <- radf_mc_cv(n = 100)

tidy(mc)

# Change the minimum window and the number of simulations
mc2 <- radf_mc_cv(n = 100, nrep = 600, minw = 20)

tidy(mc2)

mdist <- radf_mc_distr(n = 100, nrep = 1000)

autoplot(mdist)

# Apply the critical values to actual data
rsim_data <- radf(sim_data, minw = 20)
autoplot(rsim_data, cv = mc2)

Description

radf_recovery implements the reverse-regression dating of Phillips & Shi (2014). It reverses the series and runs the existing bsadf recursion of radf() on it. It then locates the first up-crossing of a critical-value boundary calibrated for the reversal, which is the market recovery date, and the next down-crossing, which is the crisis (collapse) origination date in the original series. Both are mapped back to the original time index.

Usage

radf_recovery(
  data,
  minw = NULL,
  lag = 0,
  nrep = 1000L,
  sig_lvl = 95,
  seed = NULL
)

Arguments

Details

The function returns two dates for each series. f_c is the crisis origination (collapse-onset) date, a reverse-regression alternative to the collapse date that datestamp already dates from the forward test. f_r is the market recovery date, and f_c <= f_r always holds by construction, because the down-crossing is searched only after the up-crossing. If no up-crossing is found, neither date is identified (NA, detected = FALSE). If an up-crossing is found but no later down-crossing occurs before the reverse-time sample ends, f_c is NA and censored = TRUE, which means that the crisis origination predates the observed sample.

Value

An object of class radf_recovery_obj: a list with f_c and f_r (the estimated dates, NA if not identified), detected (logical, whether an up-crossing was found at all) and censored (logical, whether f_c is left-censored by the start of the reverse-time sample).

Caveats

Validation status (2026-08-10). f_r (the recovery date) validates well against synthetic collapse-then-recovery data. Its bias is in the same range as in the Monte Carlo study of the paper, a few observations early. f_c (the crisis origination date) shows a materially larger residual bias in Monte Carlo checks. The empirical false-detection rate under a pure random-walk null (n = 100, minw = 20, 95\ than the comparable numbers for the forward tests elsewhere in this package. We found and fixed one artifact of the synthetic process during validation, where a level jump at a regime boundary produced a spurious spike. The remaining bias in f_c and the elevated false-detection rate are not fully explained. They may be genuine finite-sample noise from the literal first-down-crossing rule of the paper. The inf in its eq. 9 has no persistence requirement, so a transient dip below the boundary is enough to trigger a premature f_c. We have not ruled out a subtler implementation issue. Treat f_c and the overall detection rate as exploratory until they are validated further, and see docs/dating-and-root-inference.md for the full numbers. The function emits the same short pointer as a message when it is called (use suppressMessages to silence it) and stores it as attr(x, "caveat") on the returned object.

Note

The function returns its own class and not radf_obj, so it does not work with summary(), ⁠\link{datestamp}⁠ and tidy. It has its own print() and autoplot() methods instead. print() shows the origination and recovery dates. See vignette("naming-and-analysis", package = "exuber") for which functions fit the shared pipeline and which do not.

References

Phillips, P. C. B., & Shi, S. (2014). Financial Bubble Implosion and Reverse Regression. Cowles Foundation Discussion Paper No. 1967, Yale University. Published in Econometric Theory.

See Also

datestamp for the forward, non-reversed dating of origination and collapse that this function complements.

Other dating: dating_hls(), dating_hlw(), dating_knp(), dating_pdc(), rootstamp()

Examples

# The expansion, bubble, collapse and recovery process of sim_ps1()
y <- sim_ps1(n = 100, seed = 1)
res <- radf_recovery(y, minw = 15, nrep = 200, seed = 1)
print(res)

# Plot the series with the estimated collapse (f_c) and recovery (f_r) points
autoplot(res)

Description

Computes critical values for the reverse-regression BSADF statistic that radf_recovery uses. They are calibrated to the null limiting distribution of that statistic (Theorem 1 of Phillips & Shi 2014) and not to the standard forward boundary of radf_mc_cv. The simulated null path is reversed before the recursive computation, because reversal induces an endogeneity that has no analogue in the forward regression (see the Details of radf_recovery).

Usage

radf_recovery_cv(n, minw = NULL, nrep = 1000L, seed = NULL, lag = 0)

Arguments

Value

A list of class radf_cv with a single element, bsadf_cv. It is a matrix of critical values (columns 90%, 95%, 99%) with one row for each reverse-time position, aligned in the same way as the bsadf_cv of radf_mc_cv is aligned to radf()$bsadf.

Status

Note

print() and tidy() are not yet implemented for the class of this object. recovery_cv has no tidy_radf_cv method. The objects from radf_sign_cv() and radf_tt_cv() fall back to the method of mc_cv, but that fallback does not apply here, because this object carries only bsadf_cv and not the adf_cv, sadf_cv and gsadf_cv fields that the method expects. Inspect cv$bsadf_cv directly. Calling print(cv) or tidy(cv) currently gives an error.

See Also

radf_recovery, radf_mc_cv

Other critical values: radf_common_cv(), radf_mc_cv(), radf_sb_cv(), radf_sbz_cv(), radf_sign_cv(), radf_sign_dm_cv(), radf_tt_cv(), radf_wb_cv(), radf_wb_ps_cv()

Examples

cv <- radf_recovery_cv(n = 100, minw = 20, nrep = 200)
range(cv$bsadf_cv)

Description

radf_sb_cv computes critical values for the panel recursive unit root test with the sieve bootstrap procedure of Pavlidis et al. (2016). radf_sb_distr computes the distribution.

Usage

radf_sb_cv(
  data,
  minw = NULL,
  lag = 0L,
  nboot = 500L,
  type = c("fixed", "aic", "bic"),
  max_lag = 8L,
  seed = NULL
)

radf_sb_distr(
  data,
  minw = NULL,
  lag = 0L,
  nboot = 500L,
  type = c("fixed", "aic", "bic"),
  max_lag = 8L,
  seed = NULL
)

Arguments

Value

For radf_sb_cv, a list with the critical values for the panel BSADF and panel GSADF test statistics. For radf_sb_distr, a numeric vector with the distribution of the panel GSADF statistic.

References

Pavlidis, E., Yusupova, A., Paya, I., Peel, D., Martínez-García, E., Mack, A., & Grossman, V. (2016). Episodes of exuberance in housing markets: In search of the smoking gun. The Journal of Real Estate Finance and Economics, 53(4), 419-449. doi:10.1007/s11146-015-9531-2

Pedersen, T. Q., & Schütte, E. C. M. (2020). Testing for explosive bubbles in the presence of autocorrelated innovations. Journal of Empirical Finance, 58, 207-225.

See Also

radf_mc_cv for Monte Carlo critical values and radf_wb_cv for wild Bootstrap critical values

Other critical values: radf_common_cv(), radf_mc_cv(), radf_recovery_cv(), radf_sbz_cv(), radf_sign_cv(), radf_sign_dm_cv(), radf_tt_cv(), radf_wb_cv(), radf_wb_ps_cv()

Examples

rsim_data <- radf(sim_data, lag = 1)

# Critical values should have the same lag length as \code{radf()}
sb <- radf_sb_cv(sim_data, lag = 1)

tidy(sb)

summary(rsim_data, cv = sb)

autoplot(rsim_data, cv = sb)

# Simulate distribution
sdist <- radf_sb_distr(sim_data, lag = 1, nboot = 1000)

autoplot(sdist)

# Automatic BIC lag selection instead of a fixed lag
sb_bic <- radf_sb_cv(sim_data, type = "bic")

Description

radf_sbz computes the WLS (kernel-volatility-weighted) recursive sup-ADF statistic of Harvey, Leybourne & Zu (2019), called supBZ in their notation, with wls_dfstat_grid() (internal). It returns the same shape as radf: the scalars adf, sadf and gsadf plus the full recursive paths badf and bsadf. The result therefore carries the radf_obj class, and the full summary(), datestamp, tidy and autoplot pipeline works with it when it is paired with radf_sbz_cv.

Usage

radf_sbz(data, minw = NULL, kernel = c("gaussian", "uniform"), h = NULL)

Arguments

Details

The bundled radf_sbz_union combines this statistic with the classic supDF statistic into a bootstrap-calibrated union test. supBZ alone needs a bootstrap only to be tested and not to be defined, so it splits into a statistic and a critical-value function, as most of exuber does.

Value

An object of class radf_sbz_obj/radf_obj: a list with adf, sadf and gsadf (one value per series) and badf and bsadf (matrices, one column per series).

Status

Note

The test needs the critical values from radf_sbz_cv, and neither radf_wb_cv nor radf_mc_cv applies. The null distribution of supBZ depends on the WLS weighting, so it needs its own critical-value function, which is data-dependent and uses a wild bootstrap, for the same reason as radf() and radf_wb_cv.

References

Harvey, D. I., Leybourne, S. J., & Zu, Y. (2019). Testing explosive bubbles with time-varying volatility. Econometric Reviews, 38(10), 1131-1151.

See Also

radf_sbz_cv for critical values, and radf_sbz_union for the main bootstrap union-of-rejections test of the paper, against the classic supDF statistic.

Other volatility-robust tests: cusum_test(), radf_kp(), radf_sbz_union(), radf_sign(), radf_sign_dm(), radf_tt(), ssu_test()

Examples

# Volatility triples at t = 100, then a strong explosive regime (rho = 1.03)
# runs from t = 120 to the sample end. The kernel-volatility weighting of supBZ
# costs enough power that the milder default bubble of sim_psy1() does not clear it
y <- sim_psy1(n = 200, te = 120, tf = 200, c = 0.03, alpha = 0, seed = 1,
  e = sim_vol_break(199))
res <- radf_sbz(y, minw = 20)
print(res)

cv <- radf_sbz_cv(y, minw = 20, nboot = 200, seed = 1)
summary(res, cv = cv)
tidy(res, cv = cv)
datestamp(res, cv = cv)
autoplot(res, cv = cv)

Description

radf_sbz_cv performs the HLST (2016) wild bootstrap. It is the same algorithm as radf_wb_cv, applied to the WLS/kernel-volatility statistic of radf_sbz instead of the classic supDF statistic. It generates critical values that include the time-varying badf_cv and bsadf_cv boundary that datestamp and autoplot need, and not only the three scalar critical values that summary() uses.

Usage

radf_sbz_cv(
  data,
  minw = NULL,
  nboot = 499L,
  kernel = c("gaussian", "uniform"),
  h = NULL,
  seed = NULL
)

Arguments

Value

An object of class radf_cv/sbz_cv/wb_cv: a list with the critical values adf_cv, sadf_cv and gsadf_cv (one row per series) and badf_cv and bsadf_cv (one array per series, one row per recursion point).

Status

References

Harvey, D. I., Leybourne, S. J., & Zu, Y. (2019). Testing explosive bubbles with time-varying volatility. Econometric Reviews, 38(10), 1131-1151.

See Also

radf_sbz for the statistic that this function pairs with, and radf_sbz_union for the bundled union-of-rejections test against the classic supDF statistic. You cannot obtain that test from this function and radf_wb_cv separately (see the Details of radf_sbz_union for why).

Other critical values: radf_common_cv(), radf_mc_cv(), radf_recovery_cv(), radf_sb_cv(), radf_sign_cv(), radf_sign_dm_cv(), radf_tt_cv(), radf_wb_cv(), radf_wb_ps_cv()

Examples

y <- sim_psy1(150, seed = 1)
res <- radf_sbz(y, minw = 20)
cv <- radf_sbz_cv(y, minw = 20, nboot = 200, seed = 1)
summary(res, cv = cv)
datestamp(res, cv = cv)

Description

radf_sbz_union performs the HLST (2016) wild bootstrap, the same algorithm as radf_wb_cv, jointly on the classic sup-ADF statistic (supDF, that is, the sadf of radf()) and on the WLS/kernel-volatility statistic supBZ of Harvey, Leybourne & Zu (2019). It combines them into the union-of-rejections statistic U of the paper. supBZ can have substantially higher power than supDF under many patterns of time-varying volatility, and lower power under others, for example upward volatility trends. U is designed to capture whichever of the two is more powerful for a given series.

Usage

radf_sbz_union(
  data,
  minw = NULL,
  nboot = 499L,
  kernel = c("gaussian", "uniform"),
  h = NULL,
  seed = NULL
)

Arguments

Details

The value of U, and not only its significance, is defined with a bootstrap-calibrated scaling ratio between the 95\ supDF and supBZ (Section 2.3 of the paper). The union also keeps its size guarantee (Theorem 3 of the paper) only if the joint bootstrap computes supDF and supBZ from the same resampled series in each replication. This coupling is why the function stays a single bundled function, and does not split into a statistic and a critical-value function as most of exuber does. supBZ alone has no such coupling, so it does split. See radf_sbz and radf_sbz_cv for the route that uses only supBZ, with the usual summary(), datestamp, tidy and autoplot pipeline.

Value

A list with bootstrap p-values (p_supDF, p_supBZ, p_U) and critical values (supDF_cv, supBZ_cv, U_cv) for each series.

Status

Note

This function bundles the statistic and its critical values in a single call. Unlike radf() and radf_wb_cv(), there is no separate statistic function without critical values and no other critical-value function to pair it with, because the value of U requires the bootstrap by construction (see Details).

The function returns its own class and not radf_obj, so it does not work with summary(), ⁠\link{datestamp}⁠ and tidy. It has its own print() and autoplot() methods instead. print() shows the test statistics together with their critical values, because the object bundles both. The autoplot() method compares supDF, supBZ and U with their critical values for each series. See vignette("naming-and-analysis", package = "exuber") for which functions fit the shared pipeline and which do not.

References

Harvey, D. I., Leybourne, S. J., & Zu, Y. (2019). Testing explosive bubbles with time-varying volatility. Econometric Reviews, 38(10), 1131-1151.

See Also

radf_wb_cv for the underlying wild bootstrap, which uses supDF only, radf_sbz and radf_sbz_cv for the route that uses supBZ only and has full pipeline support, and radf_tt for a heteroskedasticity-robust alternative that needs no bootstrap.

Other volatility-robust tests: cusum_test(), radf_kp(), radf_sbz(), radf_sign(), radf_sign_dm(), radf_tt(), ssu_test()

Examples

y <- sim_psy1(n = 200, te = 120, tf = 200, c = 0.03, alpha = 0, seed = 1,
  e = sim_vol_break(199))
res <- radf_sbz_union(y, nboot = 200, seed = 1)
print(res)
autoplot(res)

Description

radf_sign computes the sign-based variant of the recursive right-tailed unit root test of Harvey, Leybourne & Zu (2020). Instead of applying the (double-)supremum ADF test to the series itself, it applies it to the cumulated sign of the first differences, C_t = sum(sign(diff(y))). The sign() function removes all information about magnitudes, so the recursive DF statistic of C_t is exactly invariant to the pattern of volatility in the innovations, even when volatility changes over time. Unlike radf, the test needs no wild bootstrap to control its size under heteroskedasticity. The critical values in radf_sign_cv are pivotal, so they are computed once and not for each dataset.

Usage

radf_sign(data, minw = NULL)

Arguments

Details

The price of this invariance is power. The paper finds that the sign-based test outperforms the standard PSY test for many specifications of time-varying volatility and bubbles, but not for all of them, and the standard test can still win in some. The strategy that the paper recommends in practice is a bootstrap-based union of rejections that combines both tests. We have not implemented it (see the package's enhancement notes for the cost and benefit considerations), and this function provides the standalone sign-based test only. sadf is the single-supremum sPWY statistic (r1 = 0 fixed), and gsadf is the double-supremum sPSY statistic.

Value

An object of class radf_sign_obj/radf_obj. It is the same adf/badf/sadf/bsadf/gsadf list as for radf, computed on the sign-transformed series, and it pairs with radf_sign_cv.

Level-shift robustness

Harvey, Leybourne, Tatlow & Zu (2025) show that this test keeps its standard null distribution, the one without level shifts, in the presence of deterministic level shifts, provided that the number of shifts grows strictly more slowly than sqrt(T). The size of the shifts does not matter. This is a materially weaker requirement than the one the standard PSY test needs for size control, which restricts the number and the magnitude of the shifts jointly. In their simulations the standard test is never correctly sized once the number of shifts grows at rate sqrt(T), while this test stays close to its nominal size.

Status

Note

The test needs the critical values from radf_sign_cv, and neither radf_wb_cv nor any other bootstrap applies. The statistic is pivotal (exactly invariant to heteroskedasticity), so its critical values are simulated once and not for each dataset.

The result carries the radf_obj class. Since 2026-08-18 the full summary(), datestamp, tidy and autoplot pipeline works, because radf_sign_cv() now computes the time-varying badf_cv and bsadf_cv boundary that the last two need, and not only the three scalar critical values that summary() uses. See vignette("naming-and-analysis", package = "exuber").

References

Harvey, D. I., Leybourne, S. J., & Zu, Y. (2020). Sign-based unit root tests for explosive financial bubbles in the presence of deterministically time-varying volatility. Econometric Theory, 36(1), 122-169.

Harvey, D. I., Leybourne, S. J., Tatlow, D., & Zu, Y. (2025). Unit root tests for explosive financial bubbles in the presence of deterministic level shifts. Oxford Bulletin of Economics and Statistics, 87(5), 879-901. doi:10.1111/obes.12668

See Also

radf_sign_cv for critical values, radf_sign_dm for the recursively demeaned sign-based analogue, which has the same level-shift robustness, and radf for the standard test, which is not invariant.

Other volatility-robust tests: cusum_test(), radf_kp(), radf_sbz(), radf_sbz_union(), radf_sign_dm(), radf_tt(), ssu_test()

Examples

# Volatility triples half-way through the sample. This is the case of
# non-stationary volatility that this test is built for, and plain radf()
# over-rejects here
y <- sim_psy1(n = 200, seed = 1, e = sim_vol_break(199))
res <- radf_sign(y, minw = 20)
print(res)

cv <- radf_sign_cv(n = 200, minw = 20)
summary(res, cv = cv)
tidy(res, cv = cv)
datestamp(res, cv = cv)
autoplot(res, cv = cv)

Description

Simulates the asymptotic null distribution of the statistic of radf_sign. By Theorem 2 of Harvey, Leybourne & Zu (2020), this distribution does not depend on the volatility process at all (exact invariance). Like radf_tt_cv, and unlike radf_wb_cv, it therefore does not have to be recomputed for each dataset. A large n with the default nrep approximates the T -> Inf limit of the paper.

Usage

radf_sign_cv(n, minw = NULL, nrep = 2000L, seed = NULL)

Arguments

Details

You can check sadf_cv (single-supremum, r1 = 0 fixed) against the asymptotic (T = Inf) sPWY values in Table 1 of the paper. For minw/n = 0.1, the values at (10\ gsadf_cv (double-supremum) corresponds to the sPSY row: (2.933, 3.180, 3.655).

Value

An object of class radf_cv/sign_cv/mc_cv with the same structure as radf_mc_cv.

Status

References

Harvey, D. I., Leybourne, S. J., & Zu, Y. (2020). Sign-based unit root tests for explosive financial bubbles in the presence of deterministically time-varying volatility. Econometric Theory, 36(1), 122-169.

See Also

Other critical values: radf_common_cv(), radf_mc_cv(), radf_recovery_cv(), radf_sb_cv(), radf_sbz_cv(), radf_sign_dm_cv(), radf_tt_cv(), radf_wb_cv(), radf_wb_ps_cv()

Examples

cv <- radf_sign_cv(n = 200, minw = 20)
tidy(cv)

Description

radf_sign_dm computes the second sign-based analogue of the recursive right-tailed unit root test of Harvey, Leybourne & Zu (2020), which the paper denotes $\bar{s}PWY$/$\bar{s}PSY$. The construction is the same as in radf_sign, but it is built on a recursively (expanding-window) demeaned cumulated-sign series, Ctilde_t = sum_{i=2}^{t} (sign(diff(y)_i) - mean(sign(diff(y)_{2:i}))), and not on the raw cumulated sign that radf_sign uses.

Usage

radf_sign_dm(data, minw = NULL)

Arguments

Details

Harvey, Leybourne, Tatlow & Zu (2025) show that this statistic shares the asymptotic level-shift robustness of radf_sign (see the ⁠Level-shift robustness⁠ section of that function). It does not need Assumption 2 of the underlying HLZ (2020) theory, that the median of the innovations is zero, which is a strictly weaker requirement than the one radf_sign needs for its own invariance result. Their finite-sample simulations also find that the recursive demeaning tends to reduce the size distortion under level shifts further than radf_sign does, although both are asymptotically robust to level shifts under the same condition.

Value

An object of class radf_sign_dm_obj/radf_obj. It is the same adf/badf/sadf/bsadf/gsadf list as for radf, and it pairs with radf_sign_dm_cv.

Status

Note

The test needs the critical values from radf_sign_dm_cv, and radf_sign_cv does not apply, because it is calibrated to the non-demeaned radf_sign statistic. The statistic is pivotal like radf_sign, so no bootstrap is needed for each dataset.

The result carries the radf_obj class. Since 2026-08-18 the full summary(), datestamp, tidy and autoplot pipeline works, with the same fix as for radf_sign: radf_sign_dm_cv() now also computes badf_cv and bsadf_cv, and not only the three scalar critical values. See vignette("naming-and-analysis", package = "exuber").

References

Harvey, D. I., Leybourne, S. J., & Zu, Y. (2020). Sign-based unit root tests for explosive financial bubbles in the presence of deterministically time-varying volatility. Econometric Theory, 36(1), 122-169.

Harvey, D. I., Leybourne, S. J., Tatlow, D., & Zu, Y. (2025). Unit root tests for explosive financial bubbles in the presence of deterministic level shifts. Oxford Bulletin of Economics and Statistics, 87(5), 879-901. doi:10.1111/obes.12668

See Also

radf_sign_dm_cv for critical values, and radf_sign for the non-demeaned sign-based analogue.

Other volatility-robust tests: cusum_test(), radf_kp(), radf_sbz(), radf_sbz_union(), radf_sign(), radf_tt(), ssu_test()

Examples

# Volatility triples half-way through the sample. This is the case of
# non-stationary volatility that this test is built for, and plain radf()
# over-rejects here
y <- sim_psy1(n = 200, seed = 1, e = sim_vol_break(199))
res <- radf_sign_dm(y, minw = 20)
print(res)

cv <- radf_sign_dm_cv(n = 200, minw = 20)
summary(res, cv = cv)
tidy(res, cv = cv)
datestamp(res, cv = cv)
autoplot(res, cv = cv)

Description

Simulates the asymptotic null distribution of the statistic of radf_sign_dm. As for radf_sign_cv, this distribution does not depend on the volatility process (exact invariance, the analogue for this variant of Theorem 2 of HLZ 2020), so it does not have to be recomputed for each dataset.

Usage

radf_sign_dm_cv(n, minw = NULL, nrep = 2000L, seed = NULL)

Arguments

Value

An object of class radf_cv/sign_dm_cv/mc_cv with the same structure as radf_mc_cv.

Status

References

Harvey, D. I., Leybourne, S. J., & Zu, Y. (2020). Sign-based unit root tests for explosive financial bubbles in the presence of deterministically time-varying volatility. Econometric Theory, 36(1), 122-169.

See Also

Other critical values: radf_common_cv(), radf_mc_cv(), radf_recovery_cv(), radf_sb_cv(), radf_sbz_cv(), radf_sign_cv(), radf_tt_cv(), radf_wb_cv(), radf_wb_ps_cv()

Examples

cv <- radf_sign_dm_cv(n = 200, minw = 20)
tidy(cv)

Description

radf_tt computes the STADF and GSTADF test statistics of Kurozumi, Skrobotov & Tsarev, a heteroskedasticity-robust alternative to radf that needs no bootstrap. It time-deforms the series with a nonparametric estimate of its variance profile, after which the usual asymptotic recursive sup-ADF critical values for homoskedastic errors apply.

Usage

radf_tt(data, minw = NULL, kernel = c("uniform", "gaussian"), h = NULL)

Arguments

Details

We recommend radf_tt_cv for the critical values. They are pivotal (asymptotically free of the volatility process), so they do not have to be recomputed for each dataset, unlike a bootstrap. The wild bootstrap of Harvey, Leybourne, Sollis & Taylor in radf_wb_cv is a bootstrap-based alternative. Consider it if non-pivotality or the finite-sample robustness of the bootstrap is a specific concern.

Value

An object of class radf_tt_obj/radf_obj: the same adf/badf/sadf/bsadf/gsadf list as for radf, computed on the time-transformed series. It works with summary(), datestamp(), tidy() and autoplot() when paired with radf_tt_cv.

Status

Note

The result carries the radf_obj class. Since 2026-08-18 the full summary(), datestamp, tidy and autoplot pipeline works, because radf_tt_cv() now computes the time-varying badf_cv and bsadf_cv boundary that the last two need, and not only the three scalar critical values that summary() uses. See vignette("naming-and-analysis", package = "exuber").

References

Kurozumi, E., Skrobotov, A., & Tsarev, A. (2024). Time-Transformed Test for Bubbles under Non-stationary Volatility. Journal of Financial Econometrics. doi:10.1093/jjfinec/nbae026

See Also

radf_tt_cv for the pivotal asymptotic critical values that need no bootstrap, and radf_wb_cv for the bootstrap-based alternative (Harvey, Leybourne, Sollis & Taylor).

Other volatility-robust tests: cusum_test(), radf_kp(), radf_sbz(), radf_sbz_union(), radf_sign(), radf_sign_dm(), ssu_test()

Examples

# Volatility triples half-way through the sample. This is the case of
# non-stationary volatility that this test is built for, and plain radf()
# over-rejects here
y <- sim_psy1(n = 200, seed = 1, e = sim_vol_break(199))
res <- radf_tt(y, minw = 20)
print(res)

cv <- radf_tt_cv(n = 200, minw = 20)
summary(res, cv = cv)
tidy(res, cv = cv)
datestamp(res, cv = cv)
autoplot(res, cv = cv)

Description

This is the dedicated critical-value function for radf_tt. It simulates the asymptotic null distribution of the GLS-demeaned recursive sup-ADF statistic that radf_tt uses. By Theorem 1 of Kurozumi, Skrobotov & Tsarev, this distribution is free of the volatility process (pivotal). Unlike radf_wb_cv, it therefore does not have to be recomputed for each dataset. A large n with the default nrep approximates well the T -> Inf limit used in the paper.

Usage

radf_tt_cv(n, minw = NULL, nrep = 2000L, seed = NULL)

Arguments

Details

You can check the sadf_cv column (STADF, that is, r1 = 0 fixed) against the asymptotic values of Whitehouse (2019) quoted in footnote 4 of Kurozumi, Skrobotov & Tsarev. For minw/n = 0.1, the values at (10\ 1\ GSTADF (gsadf_cv). The paper gives its GSTADF critical values not as numbers in the text but only as "easily computed from" the code of the authors in R.

Value

An object of class radf_cv/tt_cv/mc_cv with the same structure as radf_mc_cv: the scalars adf_cv, sadf_cv and gsadf_cv for each level plus the sequences badf_cv and bsadf_cv. You can use it wherever a radf_cv is accepted.

Status

Note

Since 2026-08-18 the function also computes badf_cv and bsadf_cv, a time-varying boundary with one row for each recursion point. datestamp and autoplot therefore work on radf_tt results, and not only summary() and tidy. The bsadf_cv of radf_mc_cv is a cummax() across replicates, which works around the output shape of the base C++ engine. Here no such shortcut is needed, because the bsadf of gls_dfstat_grid() (internal) is already the sup over all window starts at each point, and the boundary is the per-time-point quantile across replicates. We validated the function in three ways. The last row of badf_cv is bit-identical to adf_cv, an exact identity because adf is the last point of badf in every replicate. The empirical false-alarm rate under H0 is conservative relative to the nominal level (3.3\ power on a synthetic bubble matches the established radf() and radf_mc_cv() pipeline almost exactly (18\ series. radf_sign_cv() and radf_sign_dm_cv() now compute the same boundary (see vignette("naming-and-analysis")).

References

Kurozumi, E., Skrobotov, A., & Tsarev, A. (2024). Time-Transformed Test for Bubbles under Non-stationary Volatility. Journal of Financial Econometrics. doi:10.1093/jjfinec/nbae026

See Also

Other critical values: radf_common_cv(), radf_mc_cv(), radf_recovery_cv(), radf_sb_cv(), radf_sbz_cv(), radf_sign_cv(), radf_sign_dm_cv(), radf_wb_cv(), radf_wb_ps_cv()

Examples

cv <- radf_tt_cv(n = 200, minw = 20)
tidy(cv)

# Volatility triples half-way through the sample. This is the case of
# non-stationary volatility that this test is built for, and plain radf()
# over-rejects here
y <- sim_psy1(n = 200, seed = 1, e = sim_vol_break(199))
res <- radf_tt(y, minw = 20)
datestamp(res, cv = cv)
autoplot(res, cv = cv)