Introduction
Econometric models, such as the Autoregressive Distributed Lag (ARDL)
models, are increasingly used to analyze the long-run equilibrium
relationship between a dependent variable \(y\) and a set of independent variables
\(x\) across a panel of \(N\) cross-sectional units over \(T\) time periods (Pesaran, Shin, and Smith 2001; Pesaran and Shin
1999). Before applying such models to the panel data, one has to
first determine if a long-run equilibrium relationship exists between
the variables using a cointegration test. However, panel data often
suffer from the problem of cross-sectional dependence (CD/CSD), i.e. the
units of observations are not independent of one another, and this is
typically caused by unobserved common factors that affect all units
(e.g. global financial crisis in the case of investment flows) (Pesaran 2021). Traditional residuals-based
tests like the Pedroni panel cointegration test enforce strict
exogeneity constraints and are not robust to cross-sectional dependence,
so they are vulnerable to over-rejection (Pedroni
2004). In cases of cross-sectional dependence, tests based on
structural dynamics like the Westerlund
(2007) test have to be used, which is now only available on Stata
(Persyn and Westerlund 2008).
This vignette introduces an R implementation of the Westerlund (2007) tests, , which can also be
accessed at https://github.com/bosco-hung/WesterlundTest. The
implementation follows the logic of the Stata command
xtwest developed by Persyn and
Westerlund (2008), including calculation of the four primary test
statistics (\(G_\tau\), \(G_\alpha\), \(P_\tau\), \(P_\alpha\)), asymptotic standardization
using moments from Westerlund (2007), and
a Stata-aligned residual bootstrap for handling cross-sectional
dependence. The package additionally provides reusable result objects,
individual-ECM output, summary methods, and configurable bootstrap
visualization.
This vignette starts by discussing the econometric model developed by
Westerlund (2007). Then, it discusses the
R implementations. It proceeds with some points to note and
implementation examples. It closes with some concluding remarks.
The econometric model
The Westerlund test is based on a structural equation known as the
Conditional Error Correction Model (ECM):
\[\begin{equation}
\Delta y_{it} = \underbrace{\delta'_i
d_t}_{\text{Deterministics}} + \underbrace{\alpha_i (y_{i,t-1} -
\beta'_i x_{i,t-1})}_{\text{Long-Run Equilibrium}} +
\underbrace{\sum_{j=1}^{p_i} \alpha_{ij} \Delta y_{i,t-j} +
\sum_{j=-q_i}^{p_i} \gamma_{ij} \Delta x_{i,t-j}}_{\text{Short-Run
Dynamics}} + e_{it}
(\#eq:ECM)
\end{equation}\]
where its core lies in the long-run equilibrium term \(\alpha_i (y_{i,t-1} - \beta'_i
x_{i,t-1})\). The expression \(y_{i,t-1} - \beta'_i x_{i,t-1}\)
represents the deviation from equilibrium in the previous time period.
If variables \(y\) and \(x\) are cointegrated, they move together in
the long run according to the relationship \(y_{it} = \beta_i x_{it}\). Therefore, the
residual of \(y_{i,t-1} - \beta'_i
x_{i,t-1}\) represents the “error” or “disequilibrium” at time
\(t-1\). The parameter \(\alpha_i\) determines how the system
responds to that disequilibrium. If \(\alpha_i
< 0\) (Error Correction), the system is stable. If \(y_{t-1}\) was “too high” relative to \(x_{t-1}\) (positive error), the negative
\(\alpha_i\) forces \(\Delta y_{it}\) to be negative. The
variable \(y\) falls to restore
equilibrium. On the other hand, if \(\alpha_i
= 0\) (No Error Correction), the change in \(y\) (\(\Delta
y_{it}\)) does not depend on the previous level. The variables
drift apart indiscriminately. This implies no cointegration.
The summation terms \(\sum_{j=1}^{p_i}
\phi_{ij} \Delta y_{i,t-j} + \sum_{j=-q_i}^{p_i} \gamma_{ij} \Delta
x_{i,t-j}\) represent the short-run shocks. These terms \(\sum_{j=1}^{p_i} \phi_{ij} \Delta
y_{i,t-j}\) capture the inertia in the system. By including
sufficient lags of \(\Delta y\) and
\(\Delta x\), we ensure that the
regression residual \(e_{it}\) is free
of serial correlation (white noise), so as to ensure the validity of the
OLS \(t\)-statistics. On the other
hand, the term \(\sum_{j=-q_i}^{0} \gamma_{ij}
\Delta x_{i,t-j}\) includes current and future differences of
\(x\) for correcting for endogeneity.
By projecting the error onto the leads and lags of \(\Delta x\), the regressors become strictly
exogenous with respect to the error term, allowing for valid asymptotic
inference on \(\alpha_i\).
The term \(\delta'_i d_t\)
accounts for trends in the data that are not related to the stochastic
relationship between \(x\) and \(y\). In Case 1 (\(d_t = 0\)), there is no intercept. The data
has a zero mean; In Case 2 (\(d_t =
\{1\}\)), there is a constant intercept. The data fluctuates
around a non-zero level; Finally, in Case 3 (\(d_t = \{1, t\}\)), there is both an
intercept and linear trend. The data drifts upward or downward over
time.
Hypothesis Testing Framework
The statistical test tests the value of \(\alpha_i\) derived from the equation
outlined in the following. Under the Null Hypothesis (\(H_0\)), there is no error correction,
i.e. there is no cointegration exists for any unit:
\[\begin{equation}
H_0: \alpha_i = 0 \quad \text{for all } i = 1 \dots N
(\#eq:H0)
\end{equation}\]
The Alternative Hypothesis (\(H_1\))
depends on which statistic is used (Group vs. Panel): In the case of
group mean statistics (\(G_\tau,
G_\alpha\)), they do not assume a common speed of adjustment.
They test whether cointegration exists for at least one unit in the
panel:
\[\begin{equation}
H_1^G: \alpha_i < 0 \quad \text{for at least one } i (\#eq:H1-G)
\end{equation}\]
In the case of panel mean statistics (\(P_\tau, P_\alpha\)), they “pool” the
information across the cross-sectional dimension and test whether
cointegration exists for the entire panel as a whole, assuming a
homogeneous speed of adjustment:
\[\begin{equation}
H_1^P: \alpha_i = \alpha < 0 \quad \text{for all } i (\#eq:H1-P)
\end{equation}\]
Code description
The code features a main driver: westerlund_test. This
primary interface connects the functions to perform input validation and
data sorting, as well as to handle the branching logic between
asymptotic and bootstrap inference, which are summarized in the
structure diagram below. It takes the following arguments:
westerlund_test(data,
yvar,
xvars,
idvar,
timevar,
constant = FALSE,
trend = FALSE,
lags = 1, # Single integer or c(min, max)
leads = NULL, # Single integer or c(min, max)
westerlund = FALSE,
aic = TRUE,
bootstrap = -1, # <= 0 means no bootstrap
indiv.ecm = FALSE,
lrwindow = 2,
seed = NULL,
verbose = FALSE)
where data is the panel data input, yvar is
a string specifying the dependent variable, and xvars is a
character vector of regressors entering the long-run relationship. Up to
six regressors are supported by the asymptotic moment tables, while
westerlund=TRUE restricts the model to at most one
regressor and requires at least a constant. The idvar and
timevar arguments identify the panel and time dimensions.
Inclusion of an intercept or linear trend is controlled by
constant and trend; setting
trend=TRUE requires constant=TRUE, and the
trend is constructed within each unit as the sequence \(1,\dots,T_i\), matching the Stata
implementation.
The short-run lag length \(p\) is
defined by lags, which defaults to 1, and the lead length
\(q\) is defined by leads,
which defaults to 0 when NULL. Either can be a fixed
integer or a two-element range. A range activates unit-specific
automatic selection using AIC or BIC according to aic; with
westerlund=TRUE, the Westerlund-specific information
criterion is used. The bootstrap argument controls the
number of bootstrap replications, indiv.ecm=TRUE retains
individual ECM regression tables, and lrwindow sets the
Bartlett-kernel window. The optional integer seed makes a
bootstrap run reproducible within R; when seed=NULL, the
current R RNG state is used. The verbose argument controls
detailed console output.
The implementation enforces strict time-series continuity by checking
for time gaps: \[\begin{equation}
\Delta t = t_{i,s} - t_{i,s-1}
(\#eq:delta)
\end{equation}\]
If \(\Delta t > 1\), a “hole” is
identified, and the function halts to prevent invalid differencing. The
data is sorted to ensure the vector \(Y\) follows a block-unit structure: \[\begin{equation}
Y = [y_{1,1}, y_{1,2}, \dots, y_{1,T}, \quad y_{2,1}, \dots, y_{2,T},
\quad \dots, \quad y_{N,T}]'
(\#eq:block)
\end{equation}\]
The implementation supports both balanced and unbalanced panels.
Panel lengths \(T_i\) may differ across
units, including in the bootstrap procedure. What is required is a
continuous time index within each unit after observations with missing
model variables are excluded. In applied work, a balanced panel can
still be advantageous because the Westerlund ECM is parameter intensive
and retaining more observations improves estimation precision,
particularly in small macro panels.
Econometric Estimation Logic
The WesterlundPlain function estimates ECMs and
aggregates them into panel statistics.
Unit-Level Estimation
For each unit \(i\), the code
estimates the unrestricted ECM explained in Equation @ref(eq:ECM). If
auto selection is enabled, the code minimizes the Akaike
information criterion (AIC) or the Bayesian information criterion (BIC).
If westerlund=TRUE, it uses the specific penalty:
\[\begin{equation}
AIC(p, q) = \ln\left(\frac{RSS_i}{T_i - p - q - 1}\right) +
\frac{2(p + q + \text{det} + 1)}{T_i - p_{max} - q_{max}} (\#eq:AIC)
\end{equation}\]
Long-Run Variance (HAC) Estimation
The helper calc_lrvar_bartlett computes the Newey-West
variance using a Bartlett kernel (Newey and West
1994). Following the Stata convention, autocovariances \(\hat{\gamma}_j\) are scaled by \(1/n\):
\[\begin{equation}
\hat{\omega}^2 = \hat{\gamma}_0 + 2 \sum_{j=1}^{M} \left(1 -
\frac{j}{M+1}\right) \hat{\gamma}_j, \quad \hat{\gamma}_j =
\frac{1}{n}\sum_{t=j+1}^{n} \hat{e}_t \hat{e}_{t-j} (\#eq:HAC)
\end{equation}\]
Panel Statistics Construction
Group Mean Statistics (\(G\))
The group mean statistics (\(G\))
are computed by averaging individual unit results. Following the
completion of the unit-level ECMs, the code aggregates the individual
speed-of-adjustment coefficients (\(\hat{\alpha}_i\)) and their standard
errors:
\[\begin{equation}
G_\tau = \frac{1}{N} \sum_{i=1}^N
\frac{\hat{\alpha}_i}{SE(\hat{\alpha}_i)} \quad \text{and} \quad
G_\alpha = \frac{1}{N} \sum_{i=1}^N \frac{T_i
\hat{\alpha}_i}{\hat{\alpha}_i(1)} (\#eq:G-tau)
\end{equation}\]
where \(\hat{\alpha}_i(1) =
\frac{\hat{\omega}_{ui}}{\hat{\omega}_{yi}}\) is the ratio of
long-run standard deviations derived from the Bartlett-kernel HAC
estimation described earlier.
Pooled Panel Statistics (\(P\))
For pooled statistics (\(P_\tau,
P_\alpha\)), the code partials out short-run dynamics and
deterministics from \(\Delta y_{it}\)
and \(y_{i,t-1}\) using the average lag
and lead orders calculated across all units. The pooled coefficient
\(\hat{\alpha}_{pooled}\) is then
estimated by aggregating these filtered residuals:
\[\begin{equation}
\hat{\alpha}_{pooled} = \left( \sum_{i=1}^N \sum_{t=1}^T
\tilde{y}_{i,t-1}^2 \right)^{-1} \sum_{i=1}^N \sum_{t=1}^T
\frac{1}{\hat{\alpha}_i(1)} \tilde{y}_{i,t-1} \Delta \tilde{y}_{it}
(\#eq:alpha)
\end{equation}\]
The final statistics are then constructed as:
\[\begin{equation}
P_\tau = \frac{\hat{\alpha}_{pooled}}{SE(\hat{\alpha}_{pooled})} \quad
\text{and} \quad P_\alpha = T \cdot \hat{\alpha}_{pooled} (\#eq:P-tau)
\end{equation}\]
where \(T\) represents the effective
sample size after accounting for the loss of observations due to the
chosen lag and lead lengths.
Standardization and Display
The function DisplayWesterlund converts raw statistics
\(S\) into Z-scores using asymptotic
moments \(\mu_S\) and \(\sigma_S\):
\[\begin{equation}
Z_S = \frac{\sqrt{N}(S - \mu_S)}{\sigma_S} (\#eq:Z-score)
\end{equation}\]
The moments are retrieved from hard-coded matrices indexed by
deterministic cases (1: None, 2: Constant, 3: Trend) and the number of
regressors \(K\) in the original Westerlund (2007) paper.
As the test is lower-tailed, the null of no cointegration is rejected
if the observed statistic is significantly negative.
Returned Results and Reporting
westerlund_test() returns an object of class
westerlund_test. In addition to the four raw statistics,
the object contains standardized Z-scores and asymptotic p-values,
bootstrap p-values and bootstrap distributions when requested,
unit-level estimates, individual ECM information, mean-group
coefficients and standard errors, model settings, and metadata. Commonly
used components include:
res$test_stats
res$z_scores
res$p_values
res$boot_pvals
res$unit_data
res$mean_group
res$bootstrap_distributions
The S3 methods provide compact and extended reporting:
When indiv.ecm = TRUE, unit-specific ECM coefficient
tables are retained in res$indiv_reg. Individual OLS
regression tables use residual degrees of freedom and \(t\) inference, while the mean-group
reporting tables use large-sample normal (\(z\)) inference.
Bootstrap Inference
As explained earlier, the Westerlund test can handle cross-sectional
dependence through a residual bootstrap. The implementation follows the
bootstrap logic of xtwest under the null hypothesis \(H_0:\alpha_i=0\):
Restricted Null-Model Estimation: For each unit
\(i\), the code estimates the short-run
model under the null, excluding the lagged level terms that generate
error correction. If lag or lead ranges are supplied, their optimal
values are selected from this restricted model. Residuals \(\hat e_{it}\) are centered within unit. The
differenced regressors \(\Delta
x_{it}\) are also centered within unit, with their means
calculated over all usable observations (touse) rather than
only the residual-estimation support.
Actual-Time Cluster Resampling: To preserve
contemporaneous cross-sectional dependence, bootstrap sampling is
performed on the actual time clusters rather than independently by unit
or by within-unit row position. The Stata sequence is emulated by
duplicating each panel cluster (expandcl 2), sampling the
eligible time clusters with replacement, generating the
newttt random ordering keys, aligning observations through
tussent, averaging these keys into newtt, and
then applying the common stable ordering within each unit. Because \(T_i\) is retained for each unit, unbalanced
panels remain unbalanced in the bootstrap rather than being forced to a
common minimum or maximum length.
Innovation Construction: The bootstrap
innovation \(u_{it}^*\) combines the
resampled residual with the lagged, contemporaneous, and led centered
regressor differences:
\[\begin{equation}
u_{it}^* = e_{it}^* + \sum_{j=-q}^{p} \hat{\gamma}_{ij} \Delta
x^*_{i,t-j}
(\#eq:innovation)
\end{equation}\]
The shiftNA helper pads out-of-range lags and leads with
NA. Boundary missings are then handled in the same
reconstruction sequence used by the Stata-aligned bootstrap before the
autoregressive recursion is initialized. For compatibility with the
original 2010 xtwest code, the bootstrap
lead-reconstruction loop also reproduces its currlead
behavior: when automatic lead selection differs across units, the final
processed unit’s selected lead controls the global lead loop, while
unavailable unit-specific coefficients are zero.
Recursive Simulation: The differenced dependent
variable is generated recursively from the unit-specific autoregressive
coefficients:
\[\begin{equation}
\Delta y^*_{it} = \sum_{j=1}^{p} \hat{\phi}_{ij} \Delta y^*_{i,t-j}
+ u_{it}^*
(\#eq:recursive)
\end{equation}\]
The initial lag periods are treated as a burn-in region, and each
unit is subsequently restricted using its own original \(T_i\).
Integration and Re-estimation: The simulated
differences are accumulated to recover bootstrap levels:
\[\begin{equation}
y^*_{it} = \sum_{s=1}^t \Delta y^*_{is}, \quad X^*_{it} =
\sum_{s=1}^t \Delta X^*_{is}
(\#eq:integration)
\end{equation}\]
The resulting bootstrap panel is passed back to
WesterlundPlain, which recalculates \(G_\tau\), \(G_\alpha\), \(P_\tau\), and \(P_\alpha\) for each replication.
The package intentionally reports the finite-sample-corrected
bootstrap p-value
\[\begin{equation}
p^* = \frac{\sum_{b=1}^B I(Stat^*_b \le Stat_{obs}) + 1}{B + 1}
(\#eq:robust-p)
\end{equation}\]
Visualization
The R implementation provides the plot.westerlund_test
S3 method, invoked with plot(result), to visualize
bootstrap inference. It uses to create a faceted \(2 \times 2\) grid for \(G_\tau\), \(G_\alpha\), \(P_\tau\), and \(P_\alpha\). The lower-tail critical
probability can be changed through conf_level; robust
bootstrap p-values can be shown or hidden with
show_robust_p; and plots can be saved directly with
save_path, dpi, and figsize.
Colors, line widths, density transparency, and grid display are also
configurable.
Mathematically, given \(B\)
bootstrap replications, the function estimates the empirical null
distribution using a kernel density estimator:
\[\begin{equation}
\hat{f}(s) = \frac{1}{B h} \sum_{b=1}^{B} K\left( \frac{s -
\text{Stat}^*_b}{h} \right)
(\#eq:kernel)
\end{equation}\]
where \(Stat^*_b\) represents the
simulated statistics for \(b = 1, \dots,
B\), \(K(\cdot)\) is the
Gaussian kernel, and \(h\) is the
bandwidth automatically determined by the geom_density
geometry.
The visualization logic is structured as follows. First, the area
under the kernel density curve represents the empirical null
distribution. Second, a solid vertical line marks the observed statistic
\(Stat_{obs}\). Third, a dashed
vertical line marks the empirical lower-tail critical value \(CV^*\) at the probability specified by
conf_level. Each facet reports the observed value and
critical value and, by default, the robust bootstrap p-value.
The null hypothesis is rejected at the \(\alpha\) level if the observed statistic
(solid line) lies to the left of the bootstrap critical value
(dashed line). If the density mass is located significantly to the right
of the observed statistic, the result is considered robust. Significant
overlap between the density and the observed statistic suggests that
asymptotic significance may be a false positive caused by
cross-sectional dependence.
Integration with the Other Packages
The package can be easily integrated into the workflow of the users
investigating the long-run relationship of variables in panel data.
Examinations of long-run dynamics often involve (1) testing the
existence of cross-sectional dependence, (2) testing the order of
integration, (3) testing the existence of cointegration, and finally (4)
running the error-correction model tests.
Regarding cross-sectional dependence, in the R ecosystem, the Pesaran
cross-sectional dependence test is available (Pesaran 2021) in pcdtest of the
package (Croissant, Millo, and Tappe
2006). Regarding the order of integration, the package provides
several first-generation tests, including the Levin-Lin-Chu (LLC) (Levin, Lin, and James Chu 2002),
Im-Pesaran-Shin (IPS) (Im, Pesaran, and Shin
2003), and Hadri stationarity tests (Hadri
2000). In case there is cross-sectional dependence, the package
also offers the second-generation Cross-sectionally Augmented IPS (CIPS)
test proposed by Pesaran (2007) which is
robust to cross-sectional dependence. Finally, researchers can use or to
run mean group (MG) or panel mean group (PMG) ARDL models (Zientara and Kujawski 2017; Natsiopoulos and Tzeremes
2020).
The package introduced by this vignette solves the road block
regarding the test of the existence of cointegration.
Simple Implementation Examples
This section provides a simple example of how to use
westerlund_test() for asymptotic and bootstrap inference,
inspect the returned result object with print() and
summary(), and visualize the bootstrap distribution with
the S3 plot() method.
This example uses a small synthetic panel dataset with 10
cross-sectional units and 30 time periods, where the dependent variable
y and the independent variable x1 are
generated as random normal variables. After data generation, this
example first runs the asymptotic version of the Westerlund test without
bootstrapping, and then runs the bootstrap version with a low number of
replications for demonstration purposes. The asymptotic test applies the
following settings: constant included, 1 lag, and no leads. The
bootstrap test applies the following settings: constant included,
automatic lag selection between 0 and 1 using AIC (by default), and 50
bootstrap replications.
library(Westerlund)
# 1. Generate a small synthetic panel dataset
set.seed(123)
N <- 10; T <- 30
df <- data.frame(
id = rep(1:N, each = T),
time = rep(1:T, N),
y = rnorm(N * T),
x1 = rnorm(N * T)
)
# 2. Run Asymptotic (Plain) Test
res_plain <- westerlund_test(
data = df, yvar = "y", xvars = "x1",
idvar = "id", timevar = "time",
constant = TRUE, lags = 1, leads = 0,
verbose = FALSE
)
print(res_plain)
#>
#> --- Westerlund (2007) Panel Cointegration Test ---
#> Statistic Value Z_score P_val_asymp
#> Gt -3.785 -7.064 8.074e-13
#> Ga -23.487 -9.495 1.105e-21
#> Pt -11.851 -7.316 1.279e-13
#> Pa -22.798 -13.226 3.116e-40
# 3. Run Bootstrap Test with automatic lag selection
# Note: bootstrap replications are kept low for example purposes
res_boot <- westerlund_test(
data = df, yvar = "y", xvars = "x1",
idvar = "id", timevar = "time",
constant = TRUE, lags = c(0, 1),
bootstrap = 50, seed = 123, verbose = FALSE
)
# 4. Inspect the full result summary
summary(res_boot)
#>
#> ======================================================================
#> Westerlund (2007) Panel Cointegration Test Summary
#> ======================================================================
#> Units: 10 | Average time periods: 30.00
#> Deterministic terms: constant
#> Lag/lead selection: AIC
#>
#> Statistic Value Z_score P_val_asymp P_val_boot
#> Gt -4.877 -10.91 5.139e-28 0.01961
#> Ga -25.015 -10.38 1.492e-25 0.01961
#> Pt -16.484 -11.98 2.382e-33 0.01961
#> Pa -24.070 -14.13 1.204e-45 0.01961
#>
#> Mean Group Estimates:
#> $mg_alpha
#> [1] -1.147496
#>
#> $se_mg_alpha
#> [1] 0.09689855
#>
#> $mg_betas
#> x1
#> -0.02627093
#>
#> $se_mg_betas
#> [1] 0.05726402
# 5. Visualize the Bootstrap Results
p <- plot(
res_boot,
conf_level = 0.05,
show_robust_p = TRUE
)
print(p)
