*! twostep_nardl.sthlp version 3.0.0 09mar2026

Title

twostep_nardl -- Two-step estimation of the Nonlinear
Autoregressive Distributed Lag (NARDL) model with comprehensive
post-estimation analysis

Syntax

twostep_nardl
depvar
indepvars
[if] [in]
[, options]

Model
decompose(varlist) variables to decompose into positive/negative partial sums (required)

lags(numlist) lag structure for each variable; default is lags(4)

maxlags(numlist) maximum lag for automatic selection

aic use AIC for lag selection

bic use BIC for lag selection (default when maxlags() specified)

onestep use one-step OLS estimation (alternative to two-step)

step1(method) first-step estimator: fmols (default), ols, tols, fmtols

threshold(numlist) threshold for decomposition; default is 0

bwidth(#) HAC bandwidth; default is floor(T^(1/4))

Deterministics
noconstant suppress constant (Case 1)

trendvar(varname) include a time trend (Case 4/5)

restricted restrict deterministic terms to LR relationship

exog(varlist) exogenous regressors

Reporting
noctable suppress coefficient table

noheader suppress header

nowaldtest suppress Wald test display

dots display lag selection progress

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

tsset must be set before using twostep_nardl.

Time-series operators are allowed for depvar and indepvars.

Data must be time-series (not panel).

Description

twostep_nardl implements the two-step estimation framework for the
Nonlinear Autoregressive Distributed Lag (NARDL) model developed by
Cho, Greenwood-Nimmo, and Shin (2019).

The NARDL model (Shin, Yu, and Greenwood-Nimmo, 2014)
is an asymmetric extension of the ARDL model of
Pesaran, Shin, and Smith (2001). It decomposes
changes in independent variables into positive and negative partial sums,
allowing separate long-run and short-run coefficients for increases vs. decreases.

The NARDL Error Correction Model:

Delta_y_t = rho * y_ + theta+ * x+_ + theta- * x-_
+ sum(phi_j * Delta_y_) + sum(pi+_j * Delta_x+_)
+ sum(pi-_j * Delta_x-_) + u_t

where x+_t = sum(max(Delta_x_s, 0)) and x-_t = sum(min(Delta_x_s, 0)) are
partial sum decompositions, rho is the speed of adjustment, and the long-run
multipliers are beta+ = -theta+/rho and beta- = -theta-/rho.

Linear (non-decomposed) variables:

Independent variables listed in indepvars but not in decompose()
are treated as linear variables. They enter the long-run equation in levels
and the short-run equation in first differences, without partial sum
decomposition. This allows mixing asymmetric and symmetric regressors in a
single model:

Why two-step? The single-step OLS estimator of the NARDL model suffers
from an asymptotic singularity problem that impedes the development of
standard asymptotic theory. The two-step procedure addresses this:

Step 1 (Long-run): Estimates parameters of the reparameterized long-run
cointegrating relationship using FM-OLS
(Phillips and Hansen, 1990) for k=1 or FM-TOLS
for k>1. These estimators are robust to endogeneity and serial correlation
and have mixed normal limiting distributions.

Step 2 (Short-run): Constructs the error correction term (ECT) from
Step 1 and estimates the ECM by OLS. Due to super-consistency of the
Step 1 estimator, treating the ECT as known does not affect the
asymptotic distribution of the Step 2 estimator.

The command also supports the standard one-step OLS approach
(Shin, Yu, and Greenwood-Nimmo, 2014) via the onestep option.
In one-step, the NARDL model is estimated by unrestricted OLS and long-run
coefficients are derived using the delta method via nlcom.

Options

Model

lags(numlist) specifies the lag structure. A single value applies to
all variables. Multiple values must match the number of variables (depvar
first, then each indepvar). Values of . trigger automatic lag
selection for that variable.

maxlags(numlist) sets the maximum lag for automatic selection. When
specified without lags(), all variables undergo automatic lag selection.

method(string) specifies the estimation method.
twostep (default) uses the two-step procedure.
onestep uses single-equation OLS with delta-method long-run coefficients.

step1(method) specifies the first-step estimator (two-step only).
fmols — Fully-Modified OLS (default for k=1).
ols — standard OLS (for comparison).
fmtols — Fully-Modified Transformed OLS (default for k>1).
tols — Transformed OLS.
FM-OLS and FM-TOLS are recommended as they correct for endogeneity and
serial correlation.

decompose(varlist) specifies which independent variables are decomposed
into positive and negative partial sums. This option is required.
Variables listed in indepvars but not in decompose() are treated as
linear variables — they enter the long-run equation in levels and the
short-run equation in first differences, without decomposition.

threshold(numlist) specifies the threshold value for decomposition.
Default is 0 (positive = increases, negative = decreases). Non-zero thresholds
create dead bands: only changes exceeding the threshold are classified.

bwidth(#) specifies the bandwidth for the Newey-West HAC
estimator used in FM-OLS/FM-TOLS.
Default is floor(T^{1/4}).

Deterministics

noconstant suppresses the constant (Case 1). Not recommended unless
theory dictates.

trendvar(varname) includes a time trend. Without restricted,
this gives Case 5 (unrestricted trend). With restricted, Case 4
(trend restricted to LR relationship).

restricted restricts the constant (Case 2) or trend (Case 4) to the
long-run relationship.

exog(varlist) includes additional exogenous regressors in the
short-run dynamic equation.

Deterministic cases (PSS framework)

The five cases follow
Pesaran, Shin, and Smith (2001):

Case 1: No intercept, no trend (noconstant)

Case 2: Restricted intercept, no trend (restricted)

Case 3: Unrestricted intercept, no trend (default)

Case 4: Unrestricted intercept, restricted trend (trendvar() restricted)

Case 5: Unrestricted intercept, unrestricted trend (trendvar())

Post-estimation commands

The following commands are available after twostep_nardl:

Command Description

estat ectest PSS bounds test for cointegration

estat waldtest Wald tests for LR and SR symmetry

estat diagnostics residual diagnostics panel

estat multiplier cumulative dynamic multipliers

estat halflife half-life and persistence profile

estat asymadj asymmetric adjustment speed

estat irf impulse response functions

estat ecmtable paper-style parameter estimates table (CGS 2019)

predict fitted values, residuals, ECT

estat ecmtable

estat ecmtable

Displays the estimated ECM parameters in the publication-style format of
Table 9 in Cho, Greenwood-Nimmo, and Shin (2019).
The output contains two panels:

Panel A: Long-Run Estimates — Long-run multipliers (beta+, beta-)
with standard errors and t-statistics from the FM-OLS (two-step) or
delta-method (one-step) VCE.

Panel B: Short-Run Dynamic Estimates — ECM coefficient plus all
lagged differences with standard errors and t-statistics. For two-step,
the first row is ECM(t-1); for one-step, the level variables (y_,
x+_, x-_) are shown.

A diagnostics footer reports Adjusted R-squared, RMSE, F-statistic,
number of observations, and a normality test p-value.

estat ectest

estat ectest [, siglevels(numlist) asymptotic]

Performs the Pesaran, Shin, and Smith (2001)
bounds test for the existence of a level (cointegrating) relationship.
Reports both the F-statistic (PSS) for joint significance of the level
variables and the t-statistic (BDM) on the speed of adjustment coefficient.

If the ardl package is installed, finite-sample critical values and
approximate p-values from
Kripfganz and Schneider (2020) are displayed
along with a decision matrix: .a = no rejection, . = inconclusive,
.r = rejection.

Decision rule: Reject H0 (no cointegration) if both F and t exceed
their I(1) critical values. Do not reject if either statistic is closer
to zero than the I(0) critical value. Otherwise, the result is inconclusive.

Options:

siglevels(numlist) significance levels for critical values
(default: 10 5 1).

asymptotic use asymptotic instead of finite-sample critical values.

estat waldtest

estat waldtest [, lrsymmetry srsymmetry impact
pairwise all]

Tests for asymmetry in the long-run and short-run coefficients:

Long-run symmetry: H0: beta+ = beta-, tested using delta-method
VCE for two-step or nlcom VCE for one-step.

Short-run symmetry (additive): H0: sum(pi+_j) = sum(pi-_j),
tested using robust (HC) standard errors.

Short-run symmetry (impact): H0: pi+_0 = pi-_0, the
contemporaneous impact.

estat diagnostics

estat diagnostics [, graph]

Residual diagnostic panel:

Serial correlation: Breusch-Godfrey LM test.

Heteroskedasticity: White's test.

Normality: Skewness-kurtosis test.

Functional form: Ramsey RESET test.

With graph, a CUSUM stability plot is produced.

Note: Rejection of diagnostic tests indicates potential model
misspecification. Non-rejection supports the validity of OLS inference.

estat multiplier

estat multiplier [, horizon(#) graph notable
saving(filename) poscolor(color) negcolor(color)
asymcolor(color) scheme(scheme) title(string)]

Computes cumulative dynamic multipliers showing how the dependent variable
adjusts over time to a unit positive or negative shock. The cumulative
multipliers converge to the long-run coefficients beta+ and beta-.

The difference between positive and negative multiplier paths reveals the
adjustment asymmetry pattern. With graph, a publication-quality
plot is produced showing positive, negative, and asymmetry paths.

horizon(#) number of periods; default is 40.

estat halflife

estat halflife [, horizon(#) graph]

Reports:

A. ECM Half-Life: Based on the speed of adjustment coefficient
rho. The half-life is -ln(2)/ln(1+rho), the number of periods to correct
50% of a disequilibrium shock.

B. Persistence Profile
(Pesaran and Shin, 1996):
Shows the full convergence path from 100% disequilibrium to 0%.
More informative than the simple half-life because it captures
non-monotonic convergence patterns from higher-order AR dynamics.

Warning: If rho >= 0, the ECM is NOT convergent and no half-life exists.

estat asymadj

estat asymadj [, horizon(#) graph]

Compares the adjustment speed for positive vs. negative shocks. Reports
impact multipliers, effective long-run values, half-lives for each direction,
and whether adjustment is faster for increases or decreases.

estat irf

estat irf [, horizon(#) graph]

Impulse response functions for positive and negative shocks, showing both
period-by-period responses and cumulative effects.

predict

predict [type] newvar [, xb residuals ecterm]

Options:

xb fitted values from the ECM (default).

residuals residuals from the ECM.

ecterm error correction term.

Interpretation guide

Speed of adjustment (rho)

The ECM coefficient rho (coefficient on L.ect in ADJ equation) should be:

• Negative and significant: indicates valid error correction
(convergence to long-run equilibrium).

• Between -1 and 0: monotonic convergence (most common).

• Less than -1: oscillatory convergence (overshooting).

• Non-negative: no error correction — the model is misspecified or
there is no cointegrating relationship.

Long-run coefficients (beta+, beta-)

These measure the equilibrium effect of a permanent unit increase (beta+)
or decrease (beta-) in x on y. If beta+ != beta- (Wald test rejects), the
long-run relationship is asymmetric. Example: if beta+ = 0.8 and beta- = 0.3,
a permanent unit increase in x raises y by 0.8, but a unit decrease only
lowers y by 0.3.

Short-run coefficients

The pi+ and pi- coefficients capture the immediate and lagged dynamic
responses to positive and negative shocks. Impact asymmetry (pi+_0 != pi-_0)
means that the initial response to increases differs from decreases.
Additive asymmetry (sum pi+ != sum pi-) captures overall short-run asymmetry.

Bounds test interpretation

The PSS bounds test has three possible outcomes at each significance level:

• Reject (.r): Both F and t statistics exceed I(1) bounds.
Evidence of cointegration regardless of integration order.

• No rejection (.a): Either statistic falls below I(0) bounds.
No evidence of cointegration.

• Inconclusive (.): Statistics between I(0) and I(1) bounds.
Further testing (e.g., unit root tests) needed.

Dynamic multipliers

Cumulative dynamic multipliers show the time path of adjustment to a unit
shock. They converge to the long-run coefficients. The difference between
positive and negative paths shows how much and how long asymmetry persists.

Persistence profile

PP(h) starts at 1.0 (full disequilibrium) and converges to 0 (equilibrium).
The half-life from the persistence profile may differ from the simple
ECM half-life when there are higher-order AR dynamics.

Recommended workflow

1. Verify variables are I(0) or I(1) (unit root tests).

2. Estimate with twostep_nardl.

3. Run estat ectest: bounds test for cointegration.

4. Check rho is negative and significant.

5. Run estat diagnostics: verify no serial correlation, normality, etc.

6. Run estat waldtest: test for LR and SR asymmetry.

7. Run estat multiplier, graph: visualize adjustment paths.

8. Run estat halflife, graph: assess convergence speed.

Examples

Example 1: Two-step NARDL with FMOLS (recommended)

. webuse lutkepohl2, clear

. tsset qtr

. twostep_nardl ln_consump ln_inc, decompose(ln_inc) lags(4 4) step1(fmols)

Example 2: One-step NARDL (SYG 2014 approach)

. twostep_nardl ln_consump ln_inc, decompose(ln_inc) lags(4 4) onestep

Example 3: Automatic lag selection (BIC)

. twostep_nardl ln_consump ln_inc, decompose(ln_inc) maxlags(6) bic dots

Example 4: Mixed model — asymmetric income, linear investment

. twostep_nardl ln_consump ln_inc ln_inv, decompose(ln_inc) lags(4 4 4) step1(fmols)

. * ln_inc is decomposed; ln_inv enters as a linear variable

Example 5: Bounds test for cointegration

. estat ectest

Example 6: Full post-estimation analysis

. estat waldtest, all

. estat diagnostics, graph

. estat multiplier, horizon(30) graph

. estat halflife, horizon(30) graph

. estat asymadj, horizon(30) graph

. estat irf, horizon(20) graph

Example 7: Predictions

. predict double yhat, xb

. predict double resid, residuals

. predict double ect, ecterm

Example 8: Custom threshold (asymmetric dead band)

. twostep_nardl y x, decompose(x) lags(4) threshold(0.05)

Example 9: Case 5 with trend

. twostep_nardl ln_consump ln_inc, decompose(ln_inc) lags(4 4) trendvar(qtr)

Stored results

twostep_nardl stores the following in e():

Scalars

e(N) number of observations

e(df_m) model degrees of freedom

e(df_r) residual degrees of freedom

e(r2) R-squared (step 2 / ECM)

e(r2_a) adjusted R-squared

e(F) model F-statistic

e(rmse) root mean squared error

e(ll) log-likelihood

e(rss) residual sum of squares

e(tau2) long-run variance (step 1)

e(rho) speed of adjustment coefficient

e(t_bdm) BDM t-statistic on rho

e(F_pss) PSS F-statistic for bounds test

e(case) deterministic case (1-5)

e(k) number of decomposed (asymmetric) variables

e(k_lin) number of linear (non-decomposed) variables

e(p_lag) lag order for dependent variable (in differences)

e(q_lag) maximum lag order for x variables

e(W_lr) LR symmetry Wald statistic

e(p_lr) LR symmetry p-value

e(W_sr) SR additive symmetry Wald statistic

e(p_sr) SR additive symmetry p-value

e(W_impact) impact symmetry Wald statistic

e(p_impact) impact symmetry p-value

e(lagselect) 1 if automatic lag selection was used

Macros

e(cmd) twostep_nardl

e(cmdline) command as typed

e(depvar) dependent variable name

e(xvars) independent variable(s)

e(asymvars) decomposed (asymmetric) variable(s)

e(linvars) linear (non-decomposed) variable(s)

e(method) twostep or onestep

e(step1) first-step method (twostep only)

e(step1_label) display label for step 1

e(lagstructure) lag structure string

e(pos_vars) positive partial sum variable names

e(neg_vars) negative partial sum variable names

e(ect_var) ECT variable name

e(trendvar) trend variable (if used)

e(exogvars) exogenous variables (if used)

e(predict) twostep_nardl_p

e(estat_cmd) twostep_nardl_estat

Matrices

e(b) combined coefficient vector [ADJ | LR | SR]

e(V) combined variance-covariance matrix

e(b_lr) long-run coefficients for decomposed vars (beta+, beta-)

e(V_lr) long-run VCE for decomposed vars

e(b_lin) long-run coefficients for linear vars (if any)

e(V_lin) long-run VCE for linear vars (if any)

e(b_sr) short-run coefficients (from OLS)

e(V_sr) short-run VCE

e(lags) selected lag structure

e(maxlags) maximum lag limits

Functions

e(sample) estimation sample

Requirements

Stata version: 14 or later.

Required: Time-series data with tsset.

Recommended packages:

ardl — Required for estat ectest to display Kripfganz-Schneider
finite-sample critical values. Install with: ssc install ardl

Data requirements:

• Variables should be I(0) or I(1). The ARDL/NARDL bounds test is
invalid if any variable is I(2).

• Sufficient observations for the chosen lag structure.
As a rule of thumb, N > 30 + total number of parameters.

• No structural breaks in the sample period (or control for them
with exogenous dummies).

Warnings:

• The bounds test has low power in small samples.
Consider finite-sample critical values (ardlbounds).

• If e(rho) >= 0, there is no error correction and the model
should be reconsidered.

• If diagnostic tests reject (serial correlation, heteroskedasticity),
inference may be unreliable. Consider adjusting the lag structure.

• One-step delta-method standard errors for LR coefficients may be
inaccurate in small samples. The two-step FM-OLS approach provides
more reliable LR inference.

References

Cho, J.S., Greenwood-Nimmo, M. and Shin, Y. (2019). Two-step estimation of
the nonlinear autoregressive distributed lag model. Working papers
2019rwp-154, Yonsei University, Yonsei Economics Research Institute.

Shin, Y., Yu, B. and Greenwood-Nimmo, M. (2014). Modelling asymmetric
cointegration and dynamic multipliers in a nonlinear ARDL framework. In:
Sickles, R., Horrace, W. (eds) Festschrift in Honor of Peter Schmidt.
Springer, New York, 281-314.

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

Phillips, P.C.B. and Hansen, B.E. (1990). Statistical inference in
instrumental variables regression with I(1) processes. Review of
Economic Studies, 57(1), 99-125.

Pesaran, M.H. and Shin, Y. (1996). Cointegration and speed of convergence
to equilibrium. Journal of Econometrics, 71(1-2), 117-143.

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

Kripfganz, S. and Schneider, D.C. (2023). ardl: Estimating
autoregressive distributed lag and equilibrium correction models.
Stata Journal, 23(4), 983-1012.

Authors

Dr. Merwan Roudane

merwanroudane920@gmail.com

Independent Researcher

Please cite as:

Roudane, M. (2026). twostep_nardl: Two-step estimation of the nonlinear
ARDL model. Statistical Software Components, Boston College Department of
Economics.

Acknowledgments:

This package builds upon the ardl package by Kripfganz and Schneider
for bounds testing infrastructure, and implements the methodology of
Cho, Greenwood-Nimmo, and Shin (2019).

Also see

ardl — ARDL estimation and bounds testing

ardlbounds — Critical values for bounds tests

tsset — Declare time-series data

regress — OLS regression

nlcom — Nonlinear combinations