{smcl}
{* *! version 1.1.0  29aug2026}{...}
{vieweralsosee "xtcspqardl" "help xtcspqardl"}{...}
{vieweralsosee "xtcspqardl postestimation" "help xtcspqardl_postestimation"}{...}
{vieweralsosee "xtqsh" "help xtqsh"}{...}
{vieweralsosee "xthst" "help xthst"}{...}
{viewerjumpto "Model" "xtcspqardl_methods##model"}{...}
{viewerjumpto "Step to equation map" "xtcspqardl_methods##map"}{...}
{viewerjumpto "Inference" "xtcspqardl_methods##inference"}{...}
{viewerjumpto "Long-run effects" "xtcspqardl_methods##lr"}{...}
{viewerjumpto "Two-step form" "xtcspqardl_methods##twostep"}{...}
{viewerjumpto "Diagnostics" "xtcspqardl_methods##diag"}{...}
{viewerjumpto "Changes in 1.1.0" "xtcspqardl_methods##changes"}{...}

{title:Title}

{phang}
{bf:xtcspqardl methods} {hline 2} the estimator, step by step, against the
source papers


{marker model}{title:The model}

{pstd}
Harding, Lamarche and Pesaran (2018) consider, for i = 1..N and
t = 1..T, the conditional quantile function (their equation 2.1)

{p 8 8 2}
Q_yit(tau | y_i,t-1, x_it, theta_i(tau), f_t)
= a_i(tau) + lambda_i(tau) y_i,t-1 + x_it' beta_i(tau) + f_t' gamma_i(tau),

{pstd}
with the regressors generated by (their 2.2)

{p 8 8 2}
x_it = mu_i + Gamma_i' f_t + v_it.

{pstd}
The r unobserved factors f_t are correlated with the regressors, so a
quantile regression that ignores them is inconsistent.  Averaging (2.1)
and (2.2) across units and solving for f_t gives their (2.14)-(2.15): the
factor space is spanned by current and lagged cross-sectional averages
zbar_t = (ybar_t, xbar_t')'.  Substituting and truncating the infinite
lag polynomial at pT gives the estimating equation (their 2.17)

{p 8 8 2}
y_it = a_i(tau) + lambda_i(tau) y_i,t-1 + x_it' beta_i(tau)
+ sum_{l=0..pT} zbar_t-l' delta_il(tau) + e_it(tau).


{marker map}{title:Step-to-equation map}

{pstd}
Each numbered step below is a block of code in the package.

{p2colset 5 12 14 2}{...}
{p2col :[1]}Build zbar_t = (ybar_t, xbar_t')' at each date over the units
in the estimation sample, and its lags 1..pT.  {it:HLP eq. (2.11)-(2.12)
and (2.15)}.  Implemented in {cmd:xtcspqardl.ado}; the averaging set is
the {bf:level} variables, which is what carries the factors.{p_end}

{p2col :[2]}For each unit i and each tau, minimise the check function over
the augmented regressor set,
pi_i(tau) = argmin sum_t rho_tau(y_it - X_it' pi_i).
{it:HLP eq. (2.19)-(2.20)}.  Implemented with {helpb qreg}, with the
quantile passed as a fraction so that tau is used exactly.{p_end}

{p2col :[3]}Extract vartheta_i(tau) = (lambda_i(tau), beta_i(tau)')', the
Hadamard selection in {it:HLP eq. (2.21)}, discarding the delta_il
nuisance coefficients from the parameter of interest (they are still
reported separately under {opt showcsa}).{p_end}

{p2col :[4]}Average: vartheta(tau) = (1/N) sum_i vartheta_i(tau).
{it:HLP eq. (2.21)}.{p_end}

{p2col :[5]}Variance: Vv = (1/(N-1)) sum_i (vartheta_i - vartheta)
(vartheta_i - vartheta)', so that sqrt(N)(vartheta - vartheta_0) is
asymptotically normal with that variance.  {it:HLP section 2.3 and
Theorem 4}.  The reported standard error is sqrt(Vv_jj/N).{p_end}

{p2col :[6]}Long run: theta_j(tau) = beta_j(tau)/(1 - lambda(tau)),
the plug-in of the mean-group estimates.  {it:HLP section 3, p.25}.{p_end}

{p2col :[7]}Pooling, only under {opt qccepmg}: W_i = Omega_i^-1 from the
unit-level quantile variance, vartheta_P = (sum W_i)^-1 sum W_i
vartheta_i, with the heterogeneity-robust variance
(1/N) Psi^-1 R Psi^-1, Psi = (1/N) sum W_i,
R = (1/(N-1)) sum W_i d_i d_i' W_i, d_i = vartheta_i - vartheta_MG.
{it:Pesaran (2006), eq. (67), transposed to the quantile setting}.{p_end}
{p2colreset}{...}

{pstd}
The whole vector of unit estimates is retained, so the mean-group
covariance spans coefficients {bf:and} quantiles jointly.  That is what
makes the long-run delta method and the inter-quantile tests correct.


{marker inference}{title:Inference}

{pstd}
Vv is the estimator HLP give for the heterogeneous case (their Theorem 4,
convergence at rate sqrt(N)).  Under slope homogeneity the relevant limit
is their Theorem 5, at rate sqrt(NT), with the sandwich V_psi built from
their (2.22)-(2.23); the pooled variance stored in
{cmd:e(V_pooled_hom)} under {opt qccepmg} plays that role here.  The
package reports the heterogeneous-case variance by default because it is
valid either way, only less efficient under homogeneity.

{pstd}
Whether homogeneity is tenable is reported directly, using the
quantile-regression Swamy tests of Galvao, Juhl, Montes-Rojas and Olmo
(2017) computed on the unit-level estimates:

{p 8 8 2}
S(tau) = sum_i (b_i - b_MD)' (V_i/T)^-1 (b_i - b_MD)
{c 174} chi2((n-1)k), T {c 174} inf, n fixed,{break}
D(tau) = sqrt(n) [ S/n - k ] / sqrt(2k) {c 174} N(0,1), (T,n) {c 174} inf,

{pstd}
where b_MD is the minimum-distance (inverse-variance weighted) estimator
and (V_i/T)^-1 is the inverse of the variance of b_i, taken from the
unit-level quantile regression.  With {opt unitvce(robust)} that variance
is the Powell (1986) kernel sandwich Galvao et al. use.

{pstd}
{bf:Finite-T behaviour.}  Both forms are asymptotic in T, and here they
are computed on unit regressions that also carry the lagged dependent
variable and the cross-sectional averages.  When T is short relative to
that parameter count, V_i is estimated too small, S is inflated and the
test over-rejects.  A homogeneous-slope simulation with N = 30, k = 3 and
cr_lags(2) gives D = 5.05 at T = 60, 1.84 at T = 150 and 0.92 at T = 400,
with S falling from 158 to 102 against an expected 87.  The convergence is
the expected one, but a rejection at short T is weak evidence.  Reducing
{opt cr_lags()} helps; for the homogeneity question on its own, the
static {helpb xtqsh} is the better instrument.  The mean-regression
analogue, which reports the Pesaran-Yamagata Delta and adjusted Delta, is
{helpb xthst} (Bersvendsen and Ditzen, SSC); it answers a different
question and its figures should not be read against these.

{pstd}
The Pesaran-Yamagata (2008) bias adjustment, sqrt(2k(T-k-1)/(T+1)), is
{bf:not} applied: it is derived from the finite-T moments of the
least-squares Swamy statistic under normal errors and has no counterpart
in the quantile case.  Versions of this package before 1.1.0 reported the
adjusted figure and attributed it to Pesaran and Yamagata; the statistic
being computed was in fact the Galvao et al. one, and it is now labelled
and standardized accordingly.  For the test as a standalone command, with
the Powell kernel variance and a HAC option for beta-mixing data, see
{helpb xtqsh} (Roudane, SSC).


{marker lr}{title:Long-run effects and the covariance term}

{pstd}
theta = beta/(1 - lambda) is a ratio of two estimates that come from the
{bf:same} units and are correlated.  The delta method therefore needs

{p 8 8 2}
Var(theta) = g_b^2 Var(beta) + g_l^2 Var(lambda)
+ 2 g_b g_l Cov(beta, lambda),

{pstd}
with g_b = 1/(1-lambda) and g_l = beta/(1-lambda)^2.  Dropping the third
term is not conservative in a predictable direction: Cov(beta, lambda) is
typically negative in dynamic panels, so omitting it {it:inflates} the
standard error.  In the shipped example it inflates it by about 40 per
cent.  The package computes the full quadratic form G V G'.

{pstd}
In the one-step CS-PQARDL the same machinery is reused by storing
lambda = 1 + phi internally, where phi is the speed of adjustment, so
that theta = -xi/phi = xi/(1-lambda).  The reparameterisation is affine,
so the covariance matrix and the Jacobian are unchanged.

{pstd}
The half-life solves lambda^h = 1/2 exactly,
h = ln(0.5)/ln(lambda), with
dh/dlambda = -ln(0.5)/(lambda (ln lambda)^2).  The continuous-time
approximation ln(2)/|phi| is {it:not} used; it differs materially at
moderate persistence (for phi = -0.5 the exact value is 1.00 and the
approximation is 1.39).


{marker twostep}{title:The two-step form}

{pstd}
{opt ecm} implements Ul-Durar et al. (2025) literally.

{pstd}
{bf:Step 1}, their equation (2): a single pooled panel quantile regression
of the level dependent variable on the level regressors and the
cross-sectional averages of every variable, so the reported sample is the
whole panel.  The variance comes from a block bootstrap that resamples
whole cross-sectional units, which is the panel-correct analogue of the
bootstrap the paper uses; {opt reps(0)} falls back on the asymptotic
Koenker sparsity variance.  Reported with the pseudo R1 and the joint Wald
statistic, i.e. their Table 7.

{pstd}
{bf:Step 2}, their equation (3): unit-by-unit quantile regression of the
differences on the differences and the lagged step-1 residual, averaged
over units.  This is their Table 8; the unit-level adjustment
coefficients are their Figure 3 and Table 9, available with
{opt unittable} and in the {cmd:unit} figure.

{pstd}
The two steps are separate regressions, so the joint covariance is block
diagonal by construction.  Cross-block hypothesis tests are not valid in
this mode.


{marker diag}{title:Diagnostics}

{phang}
{bf:Pseudo R1} (Koenker and Machado 1999) is 1 - Vhat/Vtilde, where Vhat
is the minimised check-function value of the fitted model and Vtilde that
of the intercept-only model.  The package pools the numerator and
denominator over the units used.

{phang}
{bf:CD} (Pesaran 2004) is
sqrt(2/(N(N-1))) sum_{i<j} sqrt(T_ij) rho_ij, standard normal under no
cross-sectional dependence.  It is reported twice: on the residuals of
the unit-level fit {it:without} the cross-sectional averages, and on the
residuals {it:with} them.  The fall in |CD| is the evidence that the
augmentation worked.  Once the averages are in, the augmented residuals
sum to about zero across units at each date, which pushes the statistic
negative by an amount of order sqrt(T)/(N-1); read the change rather than
the level.


{marker changes}{title:What changed in version 1.1.0}

{pstd}
These are corrections to the estimator, not cosmetic changes; results from
version 1.0.2 will differ.

{p2colset 5 8 10 2}{...}
{p2col :1.}{opt qccepmg} previously returned results {it:numerically
identical} to {opt qccemg}: the pooling was described in the code header
but never implemented.  It now performs the inverse-variance
combination.{p_end}
{p2col :2.}The long-run point estimate was the mean of the unit-level
ratios beta_i/(1-lambda_i), while its standard error was the delta method
for the plug-in beta/(1-lambda).  The two are different estimands, so the
estimate and its standard error did not belong together.  The plug-in
form of the paper is now used for both; it is also far more stable, since
a single unit with lambda_i near one no longer dominates.{p_end}
{p2col :3.}The delta-method variance omitted Cov(beta, lambda).{p_end}
{p2col :4.}The mean-group covariance was block diagonal across quantiles,
so every inter-quantile contrast and every test of quantile constancy used
Var(b1) + Var(b2) instead of Var(b1) + Var(b2) - 2 Cov(b1, b2).  In the
shipped example this made the contrast standard errors about 2.7 times too
large.{p_end}
{p2col :5.}The requested quantile was rounded to a whole percentage point
({cmd:quantile(round(tau*100))}), so {cmd:tau(0.125)} was silently
estimated at 0.13.  The quantile is now passed exactly.{p_end}
{p2col :6.}A regressor of interest dropped for collinearity in a unit was
recorded as a genuine coefficient of exactly zero and averaged into the
mean group.  Such unit-quantiles are now excluded and counted.{p_end}
{p2col :7.}The {opt ecm} path differenced the dependent variable a second
time and read the level block from the wrong variables, producing speeds
of adjustment outside (-2, 0).  It now implements the two-step procedure
of the source paper.{p_end}
{p2col :8.}For the CS-PQARDL the cross-sectional averages were built from
the {it:differenced} regressors rather than the levels named in
{opt lr()}, so the augmentation did not span the factor space of the
level relation.{p_end}
{p2col :9.}The half-life used the continuous-time approximation
ln(2)/|phi| and was averaged across units; it is now the exact discrete
value computed from the mean-group persistence, with a standard error.{p_end}
{p2col :10.}The options {opt srtable}, {opt level()}, {opt halflife},
{opt irf}, {opt lagsel}, {opt pmax}, {opt qmax}, {opt ec} and
{opt replace} were parsed but never used, so the short-run table was
computed and discarded and confidence intervals were fixed at 95 per cent.
The useful ones are implemented; the rest have been removed rather than
left as silent no-ops.{p_end}
{p2col :11.}No {cmd:e(b)} or {cmd:e(V)} was posted, so {helpb test},
{helpb lincom}, {helpb nlcom}, {cmd:coefplot} and {cmd:esttab} did not
work after the command.  They now do, with one equation per quantile.{p_end}
{p2col :12.}New output: pseudo R1, the joint Wald statistic, the Pesaran
CD test before and after augmentation, the Galvao-Juhl-Montes-Rojas-Olmo
slope-homogeneity tests (both the chi-squared Swamy form and the
standardized normal form), the unit-level adjustment table, and standard
errors on the half-life.{p_end}
{p2col :13.}The slope-homogeneity statistic is labelled and standardized
as Galvao et al. (2017), the quantile-regression Swamy test, rather than
as Pesaran-Yamagata (2008), whose finite-T bias adjustment applies to the
least-squares statistic and was being applied here in error.{p_end}
{p2colreset}{...}


{title:Author}

{pstd}
Dr Merwan Roudane{break}
merwanroudane920@gmail.com{break}
{browse "https://github.com/merwanroudane":github.com/merwanroudane}