*! version 1.1.0 14feb2026

Title

qardl -- Quantile Autoregressive Distributed-Lag (QARDL) Model

Syntax

qardl
depvar
indepvars
[if] [in],
tau(numlist)
[options]

Model
tau(numlist) quantile levels; required; values in (0,1)

p(#) autoregressive lag order for depvar; default 0 = auto BIC

q(#) distributed lag order for indepvars; default 0 = auto BIC

pmax(#) maximum p for BIC search; default is pmax(7)

qmax(#) maximum q for BIC search; default is qmax(7)

ecm estimate QARDL-ECM (Error Correction Model) form

noconstant suppress constant term

Rolling estimation
rolling(#) rolling window size; 0 = auto (10% of sample)

window(#) alias for rolling()

Simulation
simulate(# [#]) Monte Carlo with reps [and nobs]

Testing
waldtest(string) Wald test specification

Output
graph produce quantile process graphs

notable suppress coefficient tables

level(#) set confidence level; default is level(95)

depvar and indepvars may contain time-series operators; see tsvarlist.

Description

qardl estimates the Quantile Autoregressive Distributed-Lag (QARDL) model
proposed by Cho, Kim & Shin (2015). The QARDL model extends
the traditional Pesaran & Shin (1998) ARDL cointegration
framework into a quantile regression setting, enabling researchers to examine how
the relationship between an I(1) dependent variable and I(1) regressors varies
across different quantiles of the conditional distribution.

qardl provides the following estimation results for each specified quantile:

beta(tau) long-run cointegrating parameters

phi(tau) short-run autoregressive (AR) parameters

gamma(tau) short-run impact (distributed lag) parameters

Along with their asymptotic covariance matrices and Wald tests for parameter
constancy across quantiles.

The key advantage of QARDL over standard ARDL is the ability to detect
asymmetric and heterogeneous long-run and short-run dynamics.
For example, the speed of adjustment or the long-run equilibrium may differ
between bull and bear markets (upper vs. lower quantiles).

Options

Model

tau(numlist) specifies the quantile levels at which to estimate the model.
Values must be strictly between 0 and 1. Multiple quantiles can be specified for
across-quantile testing. At least two quantiles are needed for Wald tests of
parameter constancy. Common choices:

tau(0.25 0.5 0.75) for quartile analysis

tau(0.1 0.25 0.5 0.75 0.9) for comprehensive coverage

tau(0.1(0.1)0.9) for a fine grid of 9 quantiles

p(#) specifies the autoregressive lag order for the dependent variable
(p >= 1). If both p() and q() are left at their defaults (0),
the BIC-based automatic lag selection is performed, displaying a full BIC grid
of all (p,q) combinations with the optimal pair highlighted.

q(#) specifies the distributed lag order for the independent variables
(q >= 1). Works analogously to p().

pmax(#) specifies the maximum autoregressive lag order for the BIC
search; default is pmax(7). The grid searches all combinations from
p=1 to pmax and q=1 to qmax.

qmax(#) specifies the maximum distributed lag order for the BIC search;
default is qmax(7).

ecm requests the QARDL Error Correction Model (ECM) form. In addition
to the standard beta, phi, and gamma, the ECM form estimates:

phi*(tau) — cumulative AR coefficients that capture the short-run dynamics
in the ECM parameterization.

theta(tau) — ECM impact coefficients of the differenced independent
variables.

Rolling estimation

rolling(#) activates rolling-window QARDL estimation with the specified
window size. This produces time-varying parameter estimates and Wald test
statistics, useful for structural break analysis. If set to 0, the window
size is automatically chosen as max(10% of sample, p+q+k+10).

Simulation

simulate(# [#]) runs a Monte Carlo simulation to evaluate finite-sample
properties of the Wald tests. The first number is the number of replications;
the second (optional) is the sample size per replication (default = actual
sample size). Reports empirical rejection rates at 10%, 5%, and 1% levels.

Testing

waldtest(string) specifies the Wald test. By default, qardl
automatically performs constancy tests for all parameters. The null hypothesis
is H0: parameter(tau_i) = parameter(tau_).

Output

graph produces quantile process plots showing parameter estimates and
confidence bands across quantiles. For rolling estimation, also produces
time-varying plots.

notable suppresses the coefficient display tables while still performing
estimation and storing all results.

Remarks

Data requirements. The variables must be time series. Ensure data is
tsset before using qardl. The QARDL framework is designed for
cointegrated I(1) variables, but it does not test for unit roots or
cointegration rank — users should determine the integration order of their
variables beforehand.

Lag selection. When p() and q() are not specified, the BIC
criterion selects the optimal lag orders via OLS (at the conditional mean).
The selected orders are then used for quantile regression at all specified
tau values.

Parameter ordering. The parameter vectors (beta, phi, gamma) are stored
in quantile-first order. For example, with 2 independent variables and
3 quantiles, the beta vector is ordered as:

beta = (beta_x1(tau1), beta_x2(tau1), beta_x1(tau2), beta_x2(tau2), beta_x1(tau3), beta_x2(tau3))'

Similarly, with p=2 lags and 3 quantiles:

phi = (phi_1(tau1), phi_2(tau1), phi_1(tau2), phi_2(tau2), phi_1(tau3), phi_2(tau3))'

This matches the ordering used in the original GAUSS implementation by Cho, Kim & Shin.

Examples

Example 1: Basic QARDL with fixed lags

Setup

. webuse lutkepohl2, clear

. tsset qtr

Estimate QARDL(2,1) at three key quantiles

. qardl dln_inv dln_inc dln_consump, tau(0.25 0.5 0.75) p(2) q(1)

This produces three output tables:

1. Long-run parameters (beta): Equilibrium relationship at each quantile.
Different beta values across quantiles indicate asymmetric long-run effects.

2. Short-run AR parameters (phi): How past depvar values affect the current
level at each quantile. Labeled as L1.dln_inv, L2.dln_inv.

3. Short-run impact parameters (gamma): Contemporaneous effect of x-variables
at each quantile.

Plus a Wald test table checking if parameters are constant across quantiles.

Example 2: Automatic BIC lag selection

Let the BIC choose optimal p and q (searched over 1 to 7)

. qardl dln_inv dln_inc dln_consump, tau(0.1 0.25 0.5 0.75 0.9)

For a narrower search range:

. qardl dln_inv dln_inc dln_consump, tau(0.1 0.25 0.5 0.75 0.9) pmax(4) qmax(4)

Example 3: QARDL-ECM estimation

Estimate the Error Correction form — useful when the ECM term (speed of
adjustment) is of interest

. qardl dln_inv dln_inc dln_consump, tau(0.25 0.5 0.75) p(2) q(1) ecm

This additionally reports:

phi*(tau): Cumulative AR parameters in the ECM parameterization.

theta(tau): Impact coefficients of dx in the ECM form.

Plus ECM-specific Wald tests for constancy of phi* and theta.

Example 4: Fine quantile grid with graphs

Estimate at 9 quantiles and plot the quantile process

. qardl dln_inv dln_inc, tau(0.1(0.1)0.9) p(1) q(1) graph

Example 5: Rolling window estimation

Time-varying QARDL with a 60-observation rolling window

. qardl dln_inv dln_inc dln_consump, tau(0.25 0.5 0.75) p(1) q(1) rolling(60)

Auto window size (10% of sample)

. qardl dln_inv dln_inc, tau(0.5) p(1) q(1) rolling(0)

Example 6: Monte Carlo simulation

Evaluate Wald test size with 500 replications and n=500

. qardl dln_inv dln_inc, tau(0.25 0.5 0.75) p(1) q(1) simulate(500 500)

Example 7: Post-estimation access to stored results

Run estimation and access results

. qardl dln_inv dln_inc dln_consump, tau(0.25 0.5 0.75) p(2) q(1)

View stored beta coefficients

. matrix list e(beta)

View beta covariance matrix

. matrix list e(beta_cov)

View raw quantile regression coefficients

. matrix list e(bt_raw)

Access specific elements

. display "Beta for x1 at tau=0.25: " e(beta)[1,1]

. display "Optimal p = " e(p) " , q = " e(q)

. display "Number of quantiles = " e(ntau)

Example 8: Reproducing the Cho, Kim & Shin demo

The original GAUSS demo uses data with 2 independent variables,
selects p via BIC (gets p=2, overridden to p=3 in the demo), q=1,
and estimates at 9 quantiles (tau = 0.1, 0.2, ..., 0.9).
To replicate with qardl:

. * Load your data (e.g., from the GAUSS qardl_data.dat)

. * Assuming variables are y, x1, x2

. qardl y x1 x2, tau(0.1(0.1)0.9) p(3) q(1)

The output should match the GAUSS demo results:

Phi (with p=3): L1.y (≈0.26), L2.y (≈-0.007), L3.y (≈-0.002) at each quantile

Beta: x1 (≈6.665), x2 (≈6.667) at each quantile

Gamma: x1 (≈4.99), x2 (≈4.99) at each quantile

Example 9: Generate example data with known DGP

Generate data from the Cho, Kim & Shin (2015) DGP for testing

. qardl_makedata, n(500) seed(12345)

Estimate and compare to known true values (beta = 6.6667)

. qardl y x1 x2, tau(0.25 0.5 0.75) p(1) q(2)

Example 10: Post-estimation Wald tests

After estimation, run individual Wald tests

. qardl dln_inv dln_inc dln_consump, tau(0.25 0.5 0.75) p(1) q(1)

Test beta constancy (quantile cointegration)

. _qardl_waldtest, type(beta) tau(0.25 0.5 0.75)

Test phi constancy (short-run AR asymmetry)

. _qardl_waldtest, type(phi) tau(0.25 0.5 0.75)

Test gamma constancy (short-run impact asymmetry)

. _qardl_waldtest, type(gamma) tau(0.25 0.5 0.75)

For ECM: test phi* and theta constancy

. qardl dln_inv dln_inc dln_consump, tau(0.25 0.5 0.75) p(2) q(1) ecm

. _qardl_waldtest, type(phi_ecm) tau(0.25 0.5 0.75)

. _qardl_waldtest, type(theta) tau(0.25 0.5 0.75)

Example 11: Rolling estimation with graphs

Combine rolling window estimation with visualizations

. qardl dln_inv dln_inc dln_consump, tau(0.1 0.25 0.5 0.75 0.9) p(1) q(1) rolling(60) graph

This produces:

1. Quantile process plots for beta, phi, gamma with 95% confidence bands

2. Combined panel of all parameter plots

3. Rolling window plots of beta at the median quantile

4. Rolling Wald statistic plot with 5% critical value line

Example 12: Advanced asymmetry analysis

After estimation, run the full asymmetry diagnostic suite

. qardl dln_inv dln_inc dln_consump, tau(0.1 0.25 0.5 0.75 0.9) p(1) q(1)

. qardl_analysis

This produces:

1. Asymmetry summary tables for beta and gamma (min, max, ratio, index, Wald)

2. Coefficient heatmap across quantiles

3. Pairwise quantile equality test matrix

4. Fan charts showing all variables on a single plot

5. Asymmetry ratio bar chart (max/min for each variable)

6. Tail divergence plot with confidence intervals

7. Quantile gradient plot showing where asymmetry is concentrated

8. Combined dashboard panel

For tables only (no graphs):

. qardl_analysis, nograph

For graphs only (no tables):

. qardl_analysis, nosummary

Example 13: Automatic lag order selection with BIC grid

Let the BIC criterion select optimal (p,q) and display the full grid

. qardl dln_inv dln_inc dln_consump, tau(0.1 0.25 0.5 0.75 0.9)

This displays a table of BIC values for all combinations p=1,...,7 and q=1,...,7,
with the optimal pair marked by *. You can limit the search range:

. qardl dln_inv dln_inc dln_consump, tau(0.1 0.5 0.9) pmax(4) qmax(4)

The BIC grid follows Cho et al. (2015) and
the original MATLAB/GAUSS pqorder() function. It computes
BIC = n*ln(mean(u^2)) + k*ln(n) for each (p,q) via OLS at the conditional
mean, then selects the combination with smallest BIC value.

Output interpretation

Output tables. The qardl output displays three main coefficient
tables, each organized by quantile. Within each quantile block, all
lags/variables are listed:

Long-Run Parameters: beta(tau)

Variable Quantile Estimate Std.Err. t-stat p-value

---- tau = 0.25

x1 0.25 6.6646 0.5184 12.856 0.0000
x2 0.25 6.6669 0.1954 34.114 0.0000
---- tau = 0.50

x1 0.50 6.6660 0.4900 13.604 0.0000
x2 0.50 6.6667 0.1847 36.085 0.0000

Interpreting beta (long-run parameters). The long-run equilibrium
relationship at quantile tau is:

y = alpha(tau) + beta_1(tau)*x1 + beta_2(tau)*x2

If beta varies across quantiles, the long-run relationship is asymmetric:
the equilibrium impact of x differs depending on whether y is in the upper or
lower tail of its distribution.

Interpreting phi (short-run AR parameters). phi_i(tau) captures the
persistence of the dependent variable at quantile tau. The sum of all phi
coefficients indicates the total short-run persistence. If the sum is close to 1,
shocks are highly persistent at that quantile.

Interpreting gamma (short-run impact parameters). gamma_j(tau) measures the
immediate impact of x on y at quantile tau. If gamma varies across quantiles, the
short-run response is asymmetric.

Interpreting the Wald test. The Wald test for constancy tests:

H0: parameter(tau_i) = parameter(tau_) for all adjacent quantile pairs.

Rejection (p < 0.05) means the parameter varies significantly across
quantiles — evidence of quantile heterogeneity. This is the key result
of the QARDL approach: if the standard ARDL (conditional mean) misses these
quantile differences, it provides an incomplete picture.

Interpreting ECM parameters. With the ecm option:

phi*(tau): Cumulative autoregressive parameters after reparameterization.
These correspond to the coefficients of the differenced lagged dependent
variable in the ECM form.

theta(tau): The short-run impact coefficients in the ECM form,
combining the contemporaneous and lagged effects of dx.

Stored results

qardl stores the following in e():

Scalars

e(N) number of observations

e(p) autoregressive lag order

e(q) distributed lag order

e(k) number of independent variables

e(ntau) number of quantile levels

Matrices

e(tau) (ntau x 1) vector of quantile levels

e(beta) (k*ntau x 1) long-run parameters

e(beta_cov) (k*ntau x k*ntau) covariance of beta

e(phi) (p*ntau x 1) short-run AR parameters

e(phi_cov) (p*ntau x p*ntau) covariance of phi

e(gamma) (k*ntau x 1) short-run impact parameters

e(gamma_cov) (k*ntau x k*ntau) covariance of gamma

e(bt_raw) (ncols x ntau) raw quantile regression coefficients

e(fh) (ntau x 1) kernel density estimates

With ecm option, additionally:

e(phi_ecm) ((p-1)*ntau x 1) ECM cumulative AR parameters

e(phi_ecm_cov) covariance of phi_ecm

e(theta) (q*k*ntau x 1) ECM impact parameters

e(theta_cov) covariance of theta

With rolling() option, additionally:

e(rolling_beta) (nwindows x k*ntau) time-varying beta

e(rolling_phi) (nwindows x p*ntau) time-varying phi

e(rolling_gamma) (nwindows x k*ntau) time-varying gamma

e(rolling_wald_beta) (nwindows x 1) time-varying beta Wald stats

e(rolling_wald_phi) (nwindows x 1) time-varying phi Wald stats

e(rolling_wald_gamma) (nwindows x 1) time-varying gamma Wald stats

e(rolling_window) rolling window size

Macros

e(cmd) qardl

e(model) qardl or qardl-ecm

e(depvar) name of dependent variable

e(indepvars) names of independent variables

e(title) estimation title

e(author) Dr Merwan Roudane

e(email) merwanroudane920@gmail.com

Methods and formulas

Notation. Let y_t be the dependent variable, x_t = (x_1t,...,x_kt)' a
k-vector of I(1) regressors, t = 1,...,n. Let tau denote a quantile level
in (0,1) and s the number of specified quantiles.

1. QARDL(p,q) Model (Cho et al., 2015, eq. 1)

The conditional quantile function is:

Q_y(tau | F_t-1) = alpha(tau)
+ SUM_j=0^q-1 delta_j(tau)' * Dx_t-j
+ gamma(tau)' * x_t
+ SUM_i=1^p phi_i(tau) * y_t-i

where Dx_t-j = x_t-j - x_t-j-1 are first differences.

Parameters:

- alpha(tau): intercept at quantile tau

- delta_j(tau): (k x 1) response to the j-th lagged difference

- gamma(tau): (k x 1) coefficient on the level of x_t

- phi_i(tau): autoregressive coefficient on y_t-i

Estimation is by quantile regression (IRLS algorithm initialized with OLS).

2. Long-Run Parameters: beta(tau) (Cho et al., 2015, eq. 4)

The long-run cointegrating coefficient at quantile tau:

beta(tau) = gamma(tau) / (1 - SUM_i phi_i(tau))

If constant across quantiles, beta(tau) = beta for all tau, recovering the
standard Pesaran and Shin (1998) ARDL result.

Covariance (Theorem 1):

V_beta = Omega (kron) M^-1

where (kron) denotes the Kronecker product.

Omega is (s x s) with element:

Omega(j,i) = [min(tau_j, tau_i) - tau_j*tau_i] * b(tau_j) * b(tau_i)

where b(tau) = 1 / [(1 - SUM phi_i(tau)) * f(tau)] and f(tau) is the
conditional density of the regression error.

M is the (k x k) data-dependent matrix:

M = [X'X - X'W * inv(W'W) * W'X] / (n-q)^2

where X = x_q+1:n (level regressors) and W = (1, barw) (constant and
cumulated lagged differences).

3. Conditional Density Estimation

The density f(tau) is estimated nonparametrically:

f_hat(tau) = (1/n) * SUM_t phi(-u_t / h_B) / h_B

where phi(.) is the standard normal density and h_B is the Bofinger (1975)
bandwidth:

h_B = [4.5 * phi(z_tau)^4 / (n * (2*z_tau^2 + 1)^2)]^(1/5)

with z_tau = Phi^-1(tau). This follows the original MATLAB and GAUSS codes.

4. Short-Run AR Parameters: phi(tau) (Theorem 2)

phi(tau) = (phi_1(tau),...,phi_p(tau))' are from the quantile regression.
Joint covariance across quantiles (nuisance parameter approach):

V_phi: block (j,i) = c(tau_j, tau_i) * Psi(j,i)

where:

c(tau_j, tau_i) = [min(tau_j,tau_i) - tau_j*tau_i] / [f(tau_j)*f(tau_i)]

Psi(j,i) = inv(L_jj) * L_ji * inv(L_ii)

L is the (s*p x s*p) matrix of auxiliary regression residual products from
regressions of each lagged y_t-l on (1, x, Dx-lags).

5. Short-Run Impact Parameters: gamma(tau) (Corollary 1)

gamma(tau) are the level coefficients on x_t. Covariance by delta method:

V_gamma = Lambda * V_phi * Lambda'

Lambda is (k*s x s*p) block-diagonal with Lambda_j = beta(tau_j) * iota_p'.

6. Wald Tests for Parameter Constancy (Section 3)

Tests whether parameters are constant across quantiles:

H0: param(tau_1) = param(tau_2) = ... = param(tau_s)

Implemented via (s-1) pairwise equality restrictions R*param = 0.

Wald statistics:

W_beta = (n-1)^2 * (R*b)' * inv[R * V_beta * R'] * (R*b)
W_phi = (n-1) * (R*p)' * inv[R * V_phi * R'] * (R*p)
W_gamma = (n-1) * (R*g)' * inv[R * V_gamma * R'] * (R*g)

Under H0: W ~ chi2(df) with df = (s-1)*d. W_beta uses (n-1)^2 because
V_beta contains (n-q)^-2 from M; W_phi and W_gamma use (n-1).

Rejection (p < 0.05) = quantile heterogeneity = evidence of
quantile cointegration.

7. ECM Speed of Adjustment: rho(tau)

The error correction coefficient at each quantile is:

rho(tau) = SUM phi_i(tau) - 1

Negative rho(tau) implies convergence to equilibrium.
The SE is computed via the delta method from the phi covariance.

8. Pairwise Equality Tests (by variable)

Variable-specific pairwise tests for each (tau_i, tau_j) pair:

W_v = scale * [param_v(tau_i) - param_v(tau_j)]^2 / Var(diff)

where Var(diff) = V_ii + V_jj - 2*V_ij, and scale = (n-1)^2 for beta,
(n-1) for gamma. Under H0, W_v ~ chi2(1).

9. BIC Lag Selection (p. 290)

BIC(p,q) = n * ln(mean(u_hat^2)) + k_total * ln(n)

u_hat = OLS residuals, k_total = 1 + q*k + k + p. Grid search over
p = 1,...,pmax and q = 1,...,qmax (default 7 each).

Companion commands

The qardl package includes the following companion commands:

qardl_analysis advanced post-estimation analysis: quantile
cointegration summary table (per-variable Wald tests with QC verdict),
asymmetry diagnostics (min/max/ratio/index), coefficient heatmap,
pairwise equality tests, fan charts, asymmetry ratio bars, tail
divergence, gradient plots, ECM speed-of-adjustment bar chart, pairwise
p-value dot plots (beta & gamma), and two combined dashboards;
options: nosummary nograph nopairwise

qardl_makedata generate example data from the Cho, Kim & Shin (2015)
DGP with known true parameter values for validation

_qardl_waldtest post-estimation Wald test for individual parameter
types (beta, phi, gamma, phi_ecm, theta); returns r(wald), r(df), r(pval)

qardl_graph standalone graphing command (called automatically with
graph option, but can be called separately after estimation)

Usage examples:

. qardl_analysis

. qardl_analysis, nograph

. qardl_makedata, n(1000) seed(42)

. _qardl_waldtest, type(beta) tau(0.25 0.5 0.75)

References

Cho, J. S., Kim, T., & Shin, Y. (2015).
Quantile cointegration in the autoregressive distributed-lag modeling framework.
Journal of Econometrics, 188(2), 281-300.
https://doi.org/10.1016/j.jeconom.2015.05.003

Pesaran, M. H., & Shin, Y. (1998).
An autoregressive distributed-lag modelling approach to cointegration analysis.
Econometric Society Monographs, 31, 371-413.

Xiao, Z. (2009).
Quantile cointegrating regression.
Journal of Econometrics, 150(2), 248-260.

Koenker, R., & Bassett, G. (1978).
Regression quantiles.
Econometrica, 46(1), 33-50.

Author

Dr Merwan Roudane

Email: merwanroudane920@gmail.com

qardl v1.1.0 — February 2026

The GAUSS and MATLAB implementations were developed by Jin Seo Cho, Tae-hwan Kim
& Yongcheol Shin. This Stata implementation was translated and extended by
Dr Merwan Roudane.

Please cite this package as:

Roudane, M. (2026). qardl: Stata module for Quantile Autoregressive
Distributed-Lag estimation. Based on Cho, Kim & Shin (2015).