Introduction to wqrr: Wavelet Quantile Regression Toolbox

Dr Merwan Roudane

2026-06-01

1 Overview

wqrr bundles eight estimators that combine the Maximal Overlap Discrete Wavelet Transform (MODWT) with quantile regression so that you can study how the relationship between two (or more) time series changes both across the frequency spectrum and across the conditional distribution. All visualisations default to the MATLAB Parula colour map.

Function Method
wavelet_qr() Wavelet Quantile Regression — band-by-quantile slopes
multivariate_wqr() Multivariate WQR — multiple regressors per band
wavelet_qqr() Wavelet QQR — (θ, τ) coefficient and p-value surfaces
np_quantile_causality() Nonparametric Causality-in-Quantiles
wavelet_np_causality() Wavelet variant of the above
wavelet_mediation() Direct, interaction, paths a / b, indirect
wavelet_quantile_correlation() Wavelet Quantile Correlation with CIs
wavelet_quantile_density() Wavelet nonparametric quantile density

For plain bivariate Quantile-on-Quantile regression please use the companion CRAN package QuantileOnQuantile.

library(wqrr)

2 Example data

The package ships a 360-month synthetic dataset of macro-financial returns:

dat <- read.csv(system.file("extdata", "wqrr_data.csv", package = "wqrr"))
str(dat)
#> 'data.frame':    360 obs. of  6 variables:
#>  $ date        : chr  "1990-01-01" "1990-02-01" "1990-03-01" "1990-04-01" ...
#>  $ oil_return  : num  0.497 0.628 1.966 1.227 1.035 ...
#>  $ sp500_return: num  -0.536 -0.5738 -0.4122 -0.2486 0.0603 ...
#>  $ epu         : num  -0.944 0.203 -1.2 -0.356 -0.638 ...
#>  $ gold_return : num  -0.333 0.485 -0.224 -0.398 0.634 ...
#>  $ gdp_growth  : num  1.0198 0.0981 -0.1177 0.3648 0.6882 ...
head(dat)
#>         date oil_return sp500_return     epu gold_return gdp_growth
#> 1 1990-01-01     0.4974      -0.5360 -0.9436     -0.3328     1.0198
#> 2 1990-02-01     0.6283      -0.5738  0.2028      0.4847     0.0981
#> 3 1990-03-01     1.9663      -0.4122 -1.2003     -0.2239    -0.1177
#> 4 1990-04-01     1.2274      -0.2486 -0.3565     -0.3983     0.3648
#> 5 1990-05-01     1.0349       0.0603 -0.6377      0.6337     0.6882
#> 6 1990-06-01    -1.3055       0.1736 -0.4139     -0.1301     0.9846

For speed in this vignette we use the last 192 observations.

dat <- tail(dat, 192)

3 1. Wavelet Quantile Regression (WQR)

Decompose sp500_return and oil_return with the MODWT, aggregate into Short / Medium / Long bands, and run a quantile regression at each quantile within each band.

wqr_fit <- wavelet_qr(
  y         = dat$sp500_return,
  x         = dat$oil_return,
  quantiles = seq(0.1, 0.9, by = 0.1),
  wavelet   = "la8",
  J         = 5,
  verbose   = FALSE
)
print(wqr_fit)
#> 
#> Wavelet Quantile Regression (WQR)
#> ================================================
#>   wavelet = la8   J = 5   N = 192   q-grid = 9 
#>   Short       mean beta = -0.3563   sig(5%) = 9/9
#>   Medium      mean beta = -0.2925   sig(5%) = 9/9
#>   Long        mean beta = +0.3830   sig(5%) = 7/9

The band x quantile slope matrix:

wqr_to_matrix(wqr_fit)
#>          Short     Medium       Long
#> 0.1 -0.5184306 -0.2133165 0.60440120
#> 0.2 -0.3550541 -0.1721389 0.59839761
#> 0.3 -0.3040764 -0.2476786 0.53092488
#> 0.4 -0.3245760 -0.2781106 0.44896588
#> 0.5 -0.3086035 -0.2873292 0.37525825
#> 0.6 -0.3281262 -0.3496506 0.34851198
#> 0.7 -0.3631803 -0.4554679 0.34025097
#> 0.8 -0.3823241 -0.3418528 0.18649830
#> 0.9 -0.3226421 -0.2872026 0.01401732

Pretty heatmap with significance stars (MATLAB Parula default colour scale is provided package-wide; for WQR / mediation the Green-Orange-Red sequential scale is conventional):

plot_wqr_heatmap(wqr_fit, colorscale = "Parula")

Or just call plot():

plot(wqr_fit)

A formatted coefficient table with stars:

results_table(wqr_fit, digits = 3)
#>            Q0.10     Q0.20     Q0.30     Q0.40     Q0.50     Q0.60     Q0.70
#> Short  -0.518*** -0.355*** -0.304*** -0.325*** -0.309*** -0.328*** -0.363***
#> Medium  -0.213**  -0.172** -0.248*** -0.278*** -0.287*** -0.350*** -0.455***
#> Long    0.604***  0.598***  0.531***  0.449***  0.375***  0.349***  0.340***
#>            Q0.80     Q0.90
#> Short  -0.382*** -0.323***
#> Medium -0.342*** -0.287***
#> Long       0.186     0.014

4 2. Multivariate WQR

Same idea but with several regressors in the same QR per band:

mwqr_fit <- multivariate_wqr(
  y         = dat$sp500_return,
  X_list    = list(oil = dat$oil_return, epu = dat$epu),
  quantiles = c(0.10, 0.25, 0.50, 0.75, 0.90),
  wavelet   = "la8",
  J         = 5,
  verbose   = FALSE
)
print(mwqr_fit)
#> 
#> Multivariate Wavelet Quantile Regression
#> ================================================
#>   Y = Y   X = oil,epu 
#>   wavelet = la8   J = 5   N = 192   q-grid = 5 
#>   oil        | Short     mean beta = -0.3726  sig(5%) = 5/5
#>   oil        | Medium    mean beta = -0.1742  sig(5%) = 3/5
#>   oil        | Long      mean beta = +0.3057  sig(5%) = 4/5
#>   epu        | Short     mean beta = -0.0517  sig(5%) = 1/5
#>   epu        | Medium    mean beta = +0.3648  sig(5%) = 2/5
#>   epu        | Long      mean beta = -0.2099  sig(5%) = 2/5
plot(mwqr_fit, variable = "oil", colorscale = "Parula")

5 3. Wavelet Quantile-on-Quantile Regression (WQQR)

For the long-horizon band of the response and predictor, build the full (θ, τ) coefficient + p-value surface:

wqqr_fit <- wavelet_qqr(
  y             = dat$sp500_return,
  x             = dat$oil_return,
  quantile_step = 0.10,
  wavelet       = "la8",
  J             = 5,
  band          = "long",
  verbose       = FALSE
)
print(wqqr_fit)
#> 
#> Wavelet Quantile-on-Quantile Regression (WQQR)
#> ================================================
#>   band = Long   wf = la8   J = 5   N = 192 
#>   q-grid = 9  cells = 81 
#>   Coef range  : [ 0.014 , 0.6044 ]
#>   Sig at 5% : 63/81
plot(wqqr_fit, type = "3d", colorscale = "Parula")
plot(wqqr_fit, type = "pvalue")
plot(wqqr_fit, type = "compare")

6 4. Nonparametric Causality-in-Quantiles

Tests whether oil_return Granger-causes sp500_return in each quantile of the conditional distribution, using a Gaussian-kernel local-constant quantile estimator and the Song-Whang-Shin statistic (standard-normal critical values at 1.96 / 1.645 for 5% / 10%).

cause_fit <- np_quantile_causality(
  x         = dat$oil_return,
  y         = dat$sp500_return,
  test_type = "mean",
  q         = seq(0.1, 0.9, by = 0.1)
)
print(cause_fit)
#> 
#> Nonparametric Quantile Causality (mean)
#> ================================================
#>   N = 191   bandwidth = 0.2684   q-grid = 9 
#>   sig at 5% : 1/9   sig at 10%: 5/9
#>   max |t| = 2.007 at tau = 0.70
plot(cause_fit)

The wavelet variant (decomposes both series first):

wcr_fit <- wavelet_np_causality(
  x       = dat$oil_return,
  y       = dat$sp500_return,
  q       = seq(0.1, 0.9, by = 0.1),
  wavelet = "la8",
  J       = 5,
  verbose = FALSE
)
print(wcr_fit)
#> 
#> Wavelet Nonparametric Quantile Causality (mean)
#> ================================================
#>   wavelet = la8   J = 5 
#>   Short       sig(5%) = 0/9   max|t| = 1.711
#>   Medium      sig(5%) = 5/9   max|t| = 3.102
#>   Long        sig(5%) = 8/9   max|t| = 7.928
plot(wcr_fit, colorscale = "Parula")

7 5. Wavelet Quantile Mediation & Moderation

Direct, interaction (moderation), path a, path b, and indirect (a·b) effects for the triplet (sp500_return, oil_return, epu):

med_fit <- wavelet_mediation(
  y         = dat$sp500_return,
  x         = dat$oil_return,
  z         = dat$epu,
  quantiles = c(0.10, 0.25, 0.50, 0.75, 0.90),
  wavelet   = "la8",
  J         = 5,
  dep_name  = "SP500",
  main_name = "OIL",
  mod_name  = "EPU",
  verbose   = FALSE
)
print(med_fit)
#> 
#> Wavelet Quantile Mediation & Moderation
#> ================================================
#>   Y = SP500   X = OIL   Z = EPU 
#>   wavelet = la8   J = 5   N = 192   q-grid = 5   bands = Short,Medium,Long 
#>   direct        sig(5%) = 13 / 15
#>   interaction   sig(5%) = 1 / 15
#>   path_a        sig(5%) = 14 / 15
#>   path_b        sig(5%) = 5 / 15
#>   indirect      sig(5%) = 5 / 15
panels <- plot_mediation_panel(med_fit, colorscale = "Parula")
panels$Direct
panels$Indirect

8 6. Wavelet Quantile Correlation

Quantile correlation between MODWT detail levels of two series with parametric-bootstrap 95% CIs:

wc_fit <- wavelet_quantile_correlation(
  x         = dat$oil_return,
  y         = dat$sp500_return,
  quantiles = c(0.10, 0.25, 0.50, 0.75, 0.90),
  wavelet   = "la8",
  J         = 5,
  n_sim     = 100,
  verbose   = FALSE
)
print(wc_fit)
#> 
#> Wavelet Quantile Correlation (WQC)
#> ================================================
#>   wavelet = la8   J = 5   n_sim = 100 
#>   Significant cells : 11 / 25
head(wc_fit$results)
#>   Level Quantile Estimated_QC   CI_Lower  CI_Upper
#> 1     1     0.10   -0.1160714 -0.1160714 0.1629464
#> 2     1     0.25   -0.1388889 -0.1388889 0.1388889
#> 3     1     0.50   -0.3333333 -0.1359375 0.1250000
#> 4     1     0.75   -0.1666667 -0.1256944 0.1388889
#> 5     1     0.90   -0.1160714 -0.1160714 0.1629464
#> 6     2     0.10   -0.1160714 -0.1160714 0.1364397
plot(wc_fit, colorscale = "Parula")

9 7. Wavelet Quantile Density Estimation

qd_fit <- wavelet_quantile_density(dat$sp500_return, j0 = 4, bandwidth = 0.15)
print(qd_fit)
#> 
#> Wavelet Quantile Density Estimation
#> ================================================
#>   N = 192   coarsest j0 = 4   bandwidth = 0.15
plot(qd_fit)

10 8. Colour palettes

wqrr_colorscales(show_preview = TRUE)
#> 
#> Available wqrr color scales
#> ===========================
#>   Parula          : MATLAB R2014b default (perceptually uniform)
#>   Jet             : Classic MATLAB rainbow
#>   Turbo           : Google Turbo: improved jet
#>   BlueRed         : Diverging (blue low, red high)
#>   Sinha           : Sinha cross-quantile heatmap
#>   GreenOrangeRed  : WQR-style sequential
#>   GreenYellowRed  : Mediation-panel default
#>   Viridis         : Perceptually uniform
#>   Plasma          : High contrast
#>   Cividis         : Optimized for CVD
#>   Inferno         : Dark to bright
#>   Magma           : Dark purple to white
#>   RdBu            : Diverging red/blue
#> 
#>   Default in wqrr plots: "Parula"

The MATLAB Parula colormap is reproduced exactly from MathWorks’ 64 RGB stops:

op <- par(mar = c(0, 0, 0, 0))
image(matrix(1:256, ncol = 1), col = parula_colors(256), axes = FALSE)

par(op)

11 References