*! version 1.0.1 24aug2026

Title

ardldml -- DML-Bounds: ARDL bounds testing for cointegration with
many persistent controls

Syntax

ardldml depvar focalvar [if] [in],
controls(varlist)
[options]

Model
controls(varlist) nuisance controls, in levels; required

integrated(varlist) subset of controls() to treat as I(1)

lags(#) short-run lag order p; default lags(4)

case(#) deterministic case 1-5; default case(3)

dlags also put lags of D.focalvar in the conditional design

Cross-fitting
blocks(#) number of chronological blocks K; default blocks(5)

kfolds(#) synonym for blocks(), for ddml users

buffer(#) buffer width h in observations; default buffer(0)

First stage
mzproj(string) adaptive (default), plain or ols

penalty(string) plugin (default), low, medium, high, or a number

cpen(#) constant c in the plug-in penalty; default cpen(1.1)

ltol(#) LASSO convergence tolerance; default ltol(1e-10)

Bootstrap
breps(#) bootstrap replications B; default breps(999)

seed(#) random-number seed

bscheme(string) system (default) or fixed

nobootstrap skip the bootstrap (no critical value is produced)

nofreeze re-select the stationary support on every draw too

etafile(filename) read the Rademacher weights from a text file

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

Reporting
showfirst report which variables the first stage selected

graph plot the bootstrap null distribution

graphblocks plot the h-block cross-fitting structure

name(string) name for the graph

graphopts(string) other twoway_options

noheader suppress the header block

notable suppress the results table

nolegend suppress the interpretive footer

The data must be tsset as a single time series, and the estimation
sample must be a contiguous span with no gaps.

depvar and focalvar enter in levels. Take logs yourself
beforehand if the relation is in logs.

Description

ardldml implements DML-Bounds, the procedure of Villena (2026): a
test for a long-run (cointegrating) relationship between depvar and
focalvar when the conditioning set is high-dimensional and may itself
carry stochastic trends.

The classical ARDL bounds test of Pesaran, Shin and Smith (2001) brackets one
unknown: whether the tested regressors are I(0) or I(1). With a large,
persistent control set a second unknown appears. Partialling out controls that
carry stochastic trends can absorb part of the long-run variation that
identifies the error-correction relation. What then governs the null is not the
integration order of the original regressors but the effective integrated
count k-tilde: the number of stochastic trends that survive residualisation.
k-tilde = k puts the null at the classical I(1) endpoint, k-tilde = 0 at the
I(0) endpoint, and anything between lands inside the bracket.

The procedure has three stages. A balanced first stage projects the
differenced outcome on stationary regressors only and the lagged levels on the
control levels, cross-fitted over contiguous time blocks separated by an
h-observation buffer. The lagged levels are then residualised, and the F form
of the Wald test of joint insignificance is computed on the orthogonalised
terms. Inference comes from a restricted system wild bootstrap that
regenerates the outcome and the focal regressor jointly under the null.

There is no bounds table in this command. The critical value is generated
from your data. Comparing the statistic with the tabulated 4.94/5.73 values
over-rejects severely once the controls are integrated; see
ardldml methods.

See ardldml methods for the equation-by-equation
derivation and the step-to-equation map, and
ardldml postestimation for the diagnostics you
should run before believing any single fit.

Options

Model

controls(varlist) gives the nuisance controls, in levels. Required.
The set may be large relative to the sample; that is the point of the method.
It must not contain depvar or focalvar.

integrated(varlist) names the subset of controls() to treat as
I(1). This choice is a specification decision, not a pre-test: integrated
controls enter the differenced-target projection in first differences, and only
integrated controls can absorb a stochastic trend. Choose them on economic
grounds. If you leave the option empty every control is treated as stationary,
which forces k-tilde = k and switches the whole trend-absorption mechanism off.

lags(#) sets the short-run lag order p, the number of lagged differences
of depvar entering the conditional design. Default 4, the paper's monthly
specification. The effective sample loses p+1 observations.

case(#) selects the deterministic case, numbered as in Pesaran, Shin and
Smith and as in ardl: 1 no intercept no trend, 2 restricted intercept,
3 unrestricted intercept (default), 4 unrestricted intercept and restricted
trend, 5 unrestricted intercept and trend. This does not enter the
statistic, which is a no-intercept projection on the orthogonalised levels
either way. It sets the deterministic terms of the restricted conditional model
the bootstrap resamples from, which Algorithm 1 requires to match the empirical
specification.

dlags adds lags of D.focalvar to the conditional design, giving the
general ARDL(p,q) short-run structure. Off by default, because equation (3) of
the paper carries only the contemporaneous term. The bootstrap's marginal model
always has its own lags either way.

Cross-fitting

blocks(#) is the number of contiguous chronological blocks K. Folds
cannot be random in a time series: randomly held-out observations are adjacent
to their own training data and the first-stage error is then correlated with
the evaluation innovations. kfolds(#) is accepted as a synonym for users
coming from ddml.

buffer(#) is the buffer width h. Observations within h of an evaluation
window are dropped from that window's training set, which is what decouples the
first-stage error from the evaluation-fold innovations. It should cover the
memory of the process: with monthly data and p = 4, values of 6 to 12 are
reasonable. It costs sample; estat blocks reports
exactly how much.

First stage

mzproj(string) chooses the projection of the lagged levels on the control
levels, which is the only place trend absorption can happen.

adaptive adaptive LASSO with marginal-slope weights on the
integrated block, the paper's default

plain vanilla L1, which over-selects integrated regressors and so
induces spurious absorption

ols unpenalised, the low-dimensional corner; unstable or infeasible
when the control set is large

penalty(string) sets how the penalty for the differenced-target
equation is chosen: plugin (default) uses c*sqrt(log d/n)*sigma-hat;
low, medium and high are the three rolling-origin
cross-validation rules of the paper's Section 7.5 (the profile minimum, the
geometric midpoint, and the one-standard-error rule); a number fixes it
directly. The level projection always uses the plug-in rule regardless.

cpen(#) is the constant c in the plug-in penalty. 1.1 is the paper's
value. Lower values select more aggressively.

ltol(#) is the convergence tolerance of the coordinate-descent LASSO.
The default is tight enough that the selected support is determinate and does
not depend on the machine's floating-point arithmetic. Loosen it only to probe
how close a borderline control sits to its selection boundary.

Bootstrap

breps(#) sets the number of bootstrap replications. The paper uses 999.

seed(#) seeds Stata's random-number generator, making the run
reproducible within Stata. It cannot reproduce the reference Python
package's draws, which come from NumPy's PCG64 stream; use etafile() for
that.

bscheme(string) chooses the resampling scheme. system is Algorithm
1: one Rademacher sequence multiplies the stacked residual pair, so every
draw carries the empirical contemporaneous covariance between the conditional
and marginal innovations, and the focal regressor is regenerated and
re-cumulated. fixed holds the focal regressor at its realised path and
reweights the conditional residual alone. fixed is valid only under strong
exogeneity of the focal regressor; under exogeneity the two agree
asymptotically, so nothing is lost by keeping the default.

nobootstrap computes the statistic but no critical value. The statistic
then has no reference distribution and should not be interpreted. Useful when
sweeping specifications.

nofreeze re-selects the stationary first-stage support on every bootstrap
path as well. By default that support is frozen while the level supports
are re-selected on every path, which is what the paper specifies: freezing the
level supports instead "leaves a size distortion that grows with T".

etafile(filename) reads the ±1 wild weights from a plain text file,
B rows by n columns, whitespace or comma separated, instead of drawing them.
This exists so a bootstrap can be reproduced across languages: feed the same
weights to this command and to the reference implementation and the two
p-values agree exactly. n here is the design length, e(N).

Reporting

showfirst lists the variables the first stage selected in each
projection, and the realised plug-in penalty.

graph plots the bootstrap null distribution with the observed statistic
and the critical value marked. graphblocks plots the cross-fitting
structure. name() and graphopts() control the result.

Interpreting the output

Header. The three lines on the left name the orthogonalised objects, in
the style ddml uses: the differenced outcome after its stationary
projection, and each lagged level after its projection on the control levels.
On the right are the design length (shorter than your sample by p+1), the
number of controls and how many are being treated as I(1).

DML-Bounds F. The F form of the Wald test that both residualised level
terms are zero. Read it only against the bootstrap critical value in the
next column. It is not a classical bounds statistic and the tabulated values do
not apply to it.

bootstrap p-value. The share of bootstrap draws at or above the observed
statistic. Small p means reject the null of no level relationship, that is,
evidence of a long-run relationship conditional on your controls.

alpha. The speed of adjustment, minus the coefficient on the residualised
lagged dependent variable. A well-behaved error-correction relation has alpha
positive in this parameterisation. Values near zero mean the level terms carry
almost no information and theta is then poorly identified, which shows up as a
very large se(theta).

theta. The long-run coefficient, the ratio of the two level coefficients,
with a delta-method standard error and a Wald confidence interval. When the test
does not reject, theta is describing a relation the data do not support: read it
as descriptive, not as an estimate of an equilibrium that has been established.

selected. With showfirst, the number of control levels the m_Z
projection retained is the empirical counterpart of the effective integrated
count. Zero means nothing was absorbed and orthogonalisation did nothing; a
large number means heavy absorption, and the long-run coefficient should be
checked for instability with
estat penalty.

Remarks and practical guidance

The estimand is conditional, and that is the central caveat. ardldml
does not ask whether depvar and focalvar cointegrate unconditionally.
It asks whether they cointegrate given controls(). If a control is
itself part of the equilibrium system, partialling it out absorbs the very trend
the test is meant to detect: k-tilde is driven toward zero and the relation is
destroyed rather than cleaned. A non-rejection then reflects over-absorption,
not the absence of a long-run relationship. This cannot be diagnosed from a
single fit. Always run
estat absorption.

Do not compare the statistic with 4.94 and 5.73. The generated-regressor
remainder is O(s log d / sqrt(T)), which is order one at the sample sizes
applied work uses: at s = 3, d = 40, T = 200 it is roughly 0.78. The first-order
asymptotics have not taken hold, so the tabulated value is a not-yet-reached
asymptotic reference rather than a valid finite-sample one. In the paper's Monte
Carlo, testing against 5.73 reaches a rejection rate of 0.737 under a null with
integrated nuisance. That is why every critical value here is computed.

Choosing K and h. K trades bias against the size of each training set;
5 or 6 is usual. h should cover the memory of the process. The buffer discards
at most h(K-1) observations from each training set, which is second-order
asymptotically but very much first-order at applied sample sizes. Check the cost
with estat blocks before committing.

Sample size. The design loses p+1 observations, and the bootstrap statistic
is computed on a path that is p+1 shorter again, exactly as in the reference
implementation. With p = 4 you need a comfortable margin above 50 observations
for the first stage to have degrees of freedom left after selection.

The integrated block should stay small. The paper's validity theory holds
the integrated nuisance block fixed-dimensional; selection over a growing
integrated block is a sparse-cointegration problem it explicitly leaves open,
and it can absorb trends spuriously. Prefer a short, economically motivated
integrated() list.

Relation to other Stata commands. ardl (Kripfganz and Schneider)
implements the classical bounds test with response-surface critical values and
is the right tool when the conditioning set is small; ardldml reduces to
that setting only in spirit, not in arithmetic. ddml (Ahrens, Hansen,
Schaffer and Wiemann) shares the orthogonalisation and cross-fitting logic but
targets a low-dimensional causal parameter under random or cluster-randomised
folds with asymptotic-normal inference; it cannot produce this test. Use
estat classical to put the classical benchmark
next to the DML-Bounds result on the same sample.

Reproducibility. Results are deterministic given the data and options; the
only randomness is the bootstrap, controlled by seed(). The first-stage
LASSO is solved to a tight tolerance by default so the selected support does not
depend on the machine.

Examples

Setup: exchange-rate pass-through to U.S. prices, late Great Moderation

. net get ardldml

. use ardldml_passthrough, clear

. foreach v in cpi neer m2 ip oil {

. replace `v' = ln(`v')

. }

. tsset mdate

. keep if inrange(mdate, tm(1999m1), tm(2007m12))

The paper's reference specification

. ardldml cpi neer, controls(m2 ffr ip unrate oil gs10 baa)

integrated(m2 ip oil gs10 baa ffr) lags(4) blocks(5) buffer(6)

breps(999) seed(20260625) showfirst

The diagnostic you must not skip

. estat absorption, drop(m2 oil)

Robustness across the penalty and the projection

. estat penalty

The classical benchmark on the same sample

. estat classical

See ardldml examples for a fuller worked
session.

Stored results

ardldml stores the following in e():

Scalars

e(N) design length (sample minus p+1)

e(N_raw) observations in the estimation sample

e(N_boot) design length of each bootstrap path

e(F) DML-Bounds statistic

e(crit) bootstrap critical value

e(p) bootstrap p-value

e(B) bootstrap replications used

e(n_fail) bootstrap draws that failed and were dropped

e(corr_ev) correlation of the conditional and marginal residuals

e(alpha) speed of adjustment

e(theta) long-run coefficient

e(theta_se) its delta-method standard error

e(df_rest) number of restrictions (2)

e(nsel_dy) variables selected in the differenced-target projection

e(nsel_z) control levels selected in the m_Z projection

e(lambda) realised plug-in penalty

e(k_ctrl), e(k_int) controls, and how many treated as I(1)

e(lags), e(blocks), e(buffer), e(case) the specification

e(cpen), e(ltol), e(dlags), e(st3cons) tuning

e(level) confidence level

Macros

e(cmd) ardldml

e(cmdline) command as typed

e(depvar), e(focal) the two tested variables

e(controls), e(integrated) the conditioning set

e(mzproj), e(penrule), e(bscheme) first stage and bootstrap

e(sel_dy), e(sel_z) names of the selected variables

e(tvar) time variable

e(predict) ardldml_p

Matrices

e(b), e(V) the two level coefficients and their covariance

e(draws) the usable bootstrap statistics

e(blocks_tab) the cross-fitting partition

Functions

e(sample) marks the estimation sample

References

Villena, M. 2026. Testing cointegration with many persistent controls. SSRN
working paper 6472826.
https://doi.org/10.2139/ssrn.6472826

Pesaran, M. H., Y. Shin, and R. J. Smith. 2001. Bounds testing approaches to the
analysis of level relationships. Journal of Applied Econometrics 16:
289-326.

Chernozhukov, V., D. Chetverikov, M. Demirer, E. Duflo, C. Hansen, W. Newey, and
J. Robins. 2018. Double/debiased machine learning for treatment and structural
parameters. Econometrics Journal 21: C1-C68.

Kripfganz, S., and D. C. Schneider. 2020. Response surface regressions for
critical value bounds and approximate p-values in equilibrium correction models.
Oxford Bulletin of Economics and Statistics 82: 1456-1481.

Zou, H. 2006. The adaptive lasso and its oracle properties. Journal of the
American Statistical Association 101: 1418-1429.

McCracken, M. W., and S. Ng. 2016. FRED-MD: A monthly database for
macroeconomic research. Journal of Business and Economic Statistics 34:
574-589.

Author

Dr Merwan Roudane

merwanroudane920@gmail.com

https://github.com/merwanroudane

The method is Villena's; this is an independent implementation of it, validated
against an independent Python implementation of the same paper. Please cite the
paper alongside the software.