--- title: "Westerlund: Panel Cointegration Tests Based on Westerlund (2007)" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Westerlund: Panel Cointegration Tests Based on Westerlund (2007)} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} bibliography: references.bib --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` # 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_bounds_2001; @strom_autoregressive_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_general_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_panel_2004]. In cases of cross-sectional dependence, tests based on structural dynamics like the @westerlund_testing_2007 test have to be used, which is now only available on Stata [@persyn_error-correctionbased_2008]. This vignette introduces an R implementation of the @westerlund_testing_2007 tests, \CRANpkg{Westerlund}, which can also be accessed at [https://github.com/bosco-hung/WesterlundTest](https://github.com/bosco-hung/WesterlundTest).^[A Python version is also available at PyPI ([https://pypi.org/project/Westerlund/](https://pypi.org/project/Westerlund/)); the R and Python versions are maintained with aligned user-facing functionality, although this vignette focuses on R.] The implementation follows the logic of the Stata command `xtwest` developed by @persyn_error-correctionbased_2008, including calculation of the four primary test statistics ($G_\tau$, $G_\alpha$, $P_\tau$, $P_\alpha$), asymptotic standardization using moments from @westerlund_testing_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_testing_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})$.^[ $-\alpha_i\beta'_i$ is sometimes written as $\lambda'_i$.] 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: ```r 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_automatic_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_testing_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: ```r 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: ```r print(res) summary(res) ``` 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$: 1. **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. 2. **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. 3. **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. 4. **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$. 5. **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 \CRANpkg{ggplot2} 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 \CRANpkg{Westerlund} 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_general_2021] in `pcdtest` of the \CRANpkg{plm} package [@croissant_plm_2006]. Regarding the order of integration, the \CRANpkg{plm} package provides several first-generation tests, including the Levin-Lin-Chu (LLC) [@levin_unit_2002], Im-Pesaran-Shin (IPS) [@im_testing_2003], and Hadri stationarity tests [@hadri_testing_2000]. In case there is cross-sectional dependence, the \CRANpkg{plm} package also offers the second-generation Cross-sectionally Augmented IPS (CIPS) test proposed by @pesaran_simple_2007 which is robust to cross-sectional dependence. Finally, researchers can use \CRANpkg{PooledMeanGroup} or \CRANpkg{ARDL} to run mean group (MG) or panel mean group (PMG) ARDL models [@zientara_pooledmeangroup_2017; @natsiopoulos_ardl_2020]. The \CRANpkg{Westerlund} 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. ```{r, fig.width=6, fig.height=4} 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) # 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) # 5. Visualize the Bootstrap Results p <- plot( res_boot, conf_level = 0.05, show_robust_p = TRUE ) print(p) ``` # Points to Note The R and Python implementations are designed to reproduce the same core `xtwest` estimation and bootstrap logic, including restricted-model lag/lead selection, time-cluster resampling, unit-specific panel lengths, and the original `currlead` compatibility behavior. Exact numerical identity across software should nevertheless not be expected. Non-bootstrap results should agree up to numerical precision when the same specification and data are used. Bootstrap realizations and resulting p-values can differ because the random-number streams and low-level numerical implementations are software-specific. In addition, this package intentionally applies the finite-sample correction $(r+1)/(B+1)$ rather than the original `xtwest` reporting rule $r/B$. ## Randomization The random-number generators used by R, Python, and Stata do not generally produce identical streams. The R implementation uses `Mersenne-Twister` through `RNGkind`; supplying `seed` to `westerlund_test()` makes repeated R runs reproducible. The Python implementation likewise accepts an explicit seed. However, using the same integer seed in R, Python, and Stata does not imply replication-by-replication identical bootstrap draws, so cross-software validation should focus on the implemented algorithm and the resulting bootstrap distribution rather than matching each random draw. ## OLS Equations The implementations rely on different numerical libraries for solving the OLS systems. R uses `stats::lm`, which is based on QR decomposition with pivoting. The Python implementation requests QR decomposition for the primary ECM and pooled regressions where parity is most important, although the surrounding model-selection and numerical-library machinery still differs. Near-collinearity and floating-point conditioning can therefore generate very small cross-software differences even when the model specification is identical. ## Floating Point Drifts Since the levels are built via integration and recursion, tiny differences are magnified. # Conclusion The \CRANpkg{Westerlund} package offers an R-native implementation of the Westerlund panel cointegration tests for users who prefer an R workflow. The core estimation and bootstrap procedures are designed to track the original `xtwest` logic closely, while the package adds reusable result objects, summary and print methods, individual-ECM output, reproducible bootstrap control, and configurable visualization. Meanwhile, the current \CRANpkg{Westerlund} implementation faces the same limitation of the `xtwest` function in that they can only address cases with at most six regressors due to their reliance on a hard-coded table. Future work could explore allowing users to simulate the moment to enhance the precision of estimations and allow the use of more regressors. However, as the Westerlund test has a high parameter density which makes it power-hungry and macro panel data often suffers from a small $N$ problem, use cases involving more than six regressors are probably less likely. Other pathways of future work would be to provide users with options to choose whether to apply finite sample correction, estimate the runtime, etc.