Estimating heterogeneous panels with csdm

Scope

csdm estimates heterogeneous panel models in which unobserved common factors may generate dependence across units. The package implements:

mean group (MG) estimation (Pesaran & Smith, 1995);
static common correlated effects (CCE) (Pesaran, 2006);
dynamic CCE (DCCE) (Chudik & Pesaran, 2015); and
a cross-sectionally augmented ARDL fit with implied adjustment and long-run parameters.

These estimators share a mean-group structure: a separate regression is fitted for every eligible unit and the unit-level coefficients are averaged. CCE-based models augment those regressions with cross-sectional averages that proxy the latent common-factor space.

The package adopts parts of the model structure used by Stata’s xtdcce2 (Ditzen, 2018), but it does not claim complete command or option parity. Pooled restrictions, estimation weights, CS-DL, CS-ECM, alternative fit-level covariance estimators, and prediction on new data are currently unavailable.

Data and sample

The bundled PWT_60_07 data contain 93 countries observed annually from 1960 through 2007. The variables follow the original xtdcce2 example:

log_rgdpo: log real output;
log_hc: log human capital;
log_ck: log physical capital; and
log_ngd: log population growth plus a 5% break-even investment rate.

The 93 values of log_ngd in 1960 are missing because its construction uses a growth rate. The worked example starts in 1970 and uses 15 countries to keep the vignette quick to build.

data(PWT_60_07, package = "csdm")

keep_ids <- unique(PWT_60_07$id)[1:15]
dat <- subset(PWT_60_07, id %in% keep_ids & year >= 1970)

dim(dat)
#> [1] 570   6
range(dat$year)
#> [1] 1970 2007

The examples below are levels regressions. Calling them growth regressions would require a differenced dependent variable, which is not part of the formula used here.

Estimators

Let $i=1,\ldots,N$ index units and $t=1,\ldots,T$ index time. A heterogeneous panel model is

$y_{it}=\alpha_i+\boldsymbol{\beta}_i'\mathbf{x}_{it}+u_{it}.$

Mean group

MG estimates the equation separately for each unit and averages the identified unit coefficients:

$\widehat{\boldsymbol{\beta}}_{MG} =\frac{1}{N_e}\sum_{i\in\mathcal E}\widehat{\boldsymbol{\beta}}_i,$

where $\mathcal E$ is the common set of eligible units and $N_e$ is its size. The reported covariance is the cross-unit sample covariance of the coefficient vectors divided by $N_e$ . Inference uses a large- $N$ normal approximation.

MG permits slope heterogeneity but does not model common-factor dependence.

Common correlated effects

CCE augments each unit equation with cross-sectional averages:

$y_{it}=\alpha_i+\boldsymbol{\beta}_i'\mathbf{x}_{it} +\boldsymbol{\gamma}_i'\bar{\mathbf z}_t+e_{it},$

where $\bar{\mathbf z}_t$ commonly contains the averages of the dependent variable and regressors. Under the CCE assumptions, these averages span the relevant common-factor space and allow consistent estimation of the unit slopes (Pesaran, 2006). Adding averages does not by itself guarantee that the factor space is adequately represented; rank, dimensions, and the choice of averages still matter.

Dynamic CCE

DCCE adds dynamics and lags of the cross-sectional averages:

$y_{it}=\alpha_i +\sum_{p=1}^{P}\phi_{ip}y_{i,t-p} +\sum_{q=0}^{Q}\boldsymbol{\beta}_{iq}'\mathbf{x}_{i,t-q} +\sum_{s=0}^{S}\boldsymbol{\delta}_{is}'\bar{\mathbf z}_{t-s} +e_{it}.$

csdm_lr(type = "ardl", ylags = P, xdlags = Q) controls model lags, and csdm_csa(lags = S) controls lags of the averages. Dynamic mean-group estimates can have short- $T$ bias; CCE augmentation does not remove that bias.

CS-ARDL output

model = "cs_ardl" fits the same unit-level levels ARDL used by DCCE and then transforms its coefficients. For unit $i$ , define

$d_i=1-\sum_{p=1}^{P}\phi_{ip},\qquad \varphi_i=-d_i,\qquad \boldsymbol{\theta}_i= \frac{\sum_{q=0}^{Q}\boldsymbol{\beta}_{iq}}{d_i}.$

The package reports $\varphi_i$ as the adjustment coefficient and $\boldsymbol{\theta}_i$ as the long-run ratio before computing their mean-group summaries. Ratios are undefined when $d_i$ is numerically zero. The implementation also stores AR-root diagnostics, but it does not automatically discard unstable units.

These transformations do not fit a separate error-correction model and do not establish cointegration. The levels component is the fitted levels ARDL coefficient vector, rather than the complete transformed short-run ECM vector.

Specifying and fitting models

The formula contains the contemporaneous economic regressors. Two specification objects add CCE terms and dynamics:

csdm_csa(vars, lags) selects variables whose cross-sectional averages enter the model and their maximum lags;
csdm_lr(type, ylags, xdlags) selects lags of the dependent variable and regressors.

With vars = "_all", averages are constructed from the evaluated response and economic model-matrix columns, excluding intercepts and unit trends. An explicit character vector instead refers to numeric columns in data. vars = "_none" turns off CCE augmentation where the chosen estimator permits it.

form <- log_rgdpo ~ log_hc + log_ck + log_ngd
csa_vars <- c("log_rgdpo", "log_hc", "log_ck", "log_ngd")

static_csa <- csdm_csa(vars = csa_vars)
dynamic_csa <- csdm_csa(vars = csa_vars, lags = 3)
ardl_1_0 <- csdm_lr(type = "ardl", ylags = 1, xdlags = 0)

MG

fit_mg <- csdm(
  form, data = dat, id = "id", time = "year", model = "mg"
)
summary(fit_mg)
#> csdm summary: Mean Group Model (MG)
#> Formula: log_rgdpo ~ log_hc + log_ck + log_ngd
#> N: 15, T: 38
#> Number of obs: 570
#> R-squared (mg): 0.9449
#> CD = 3.2379, p = 0.0012
#> (For additional CD diagnostics, use cd_test())
#> 
#> Mean Group:
#>              Coef. Std. Err.      z  P>|z| Signif. CI 2.5% CI 97.5%
#> (Intercept) 6.5609    0.8648 7.5867 0.0000     ***  4.8659   8.2558
#> log_hc      0.1725    0.8530 0.2022 0.8398         -1.4993   1.8443
#> log_ck      0.3056    0.1359 2.2485 0.0245       *  0.0392   0.5720
#> log_ngd     0.7800    0.3777 2.0652 0.0389       *  0.0398   1.5203
#> 
#> Mean Group Variables: log_hc, log_ck, log_ngd
#> Cross Sectional Averaged Variables: none (lags=0)
#> 
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

CCE

fit_cce <- csdm(
  form, data = dat, id = "id", time = "year", model = "cce",
  csa = static_csa
)
summary(fit_cce)
#> csdm summary: Static Common Correlated Error Model (CCE)
#> Formula: log_rgdpo ~ log_hc + log_ck + log_ngd
#> N: 15, T: 38
#> Number of obs: 570
#> R-squared (mg): 0.9777
#> CD = -2.6758, p = 0.0075
#> (For additional CD diagnostics, use cd_test())
#> 
#> Mean Group:
#>               Coef. Std. Err.       z  P>|z| Signif. CI 2.5% CI 97.5%
#> (Intercept)  1.9003    2.1195  0.8965 0.3700         -2.2540   6.0545
#> log_hc      -1.4921    1.0152 -1.4697 0.1416         -3.4819   0.4977
#> log_ck       0.1367    0.0956  1.4298 0.1528         -0.0507   0.3240
#> log_ngd      0.8075    0.2972  2.7175 0.0066      **  0.2251   1.3899
#> 
#> Mean Group Variables: log_hc, log_ck, log_ngd
#> Cross Sectional Averaged Variables: log_rgdpo, log_hc, log_ck, log_ngd (lags=0)
#> 
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

DCCE

fit_dcce <- csdm(
  form, data = dat, id = "id", time = "year", model = "dcce",
  csa = dynamic_csa,
  lr = ardl_1_0
)
summary(fit_dcce)
#> csdm summary: Dynamic Common Correlated Error Model (DCCE)
#> Formula: log_rgdpo ~ log_hc + log_ck + log_ngd
#> N: 15, T: 38
#> Number of obs: 525
#> R-squared (mg): 0.9861
#> CD = -2.392, p = 0.0168
#> (For additional CD diagnostics, use cd_test())
#> 
#> Mean Group:
#>                  Coef. Std. Err.       z  P>|z| Signif. CI 2.5% CI 97.5%
#> (Intercept)     9.2888    8.0386  1.1555 0.2479         -6.4666  25.0441
#> log_hc         -1.9558    1.5659 -1.2490 0.2117         -5.0249   1.1132
#> log_ck          0.6666    0.2424  2.7499 0.0060      **  0.1915   1.1417
#> log_ngd        -0.3178    1.1271 -0.2819 0.7780         -2.5268   1.8913
#> lag1_log_rgdpo -0.0745    0.0555 -1.3424 0.1795         -0.1834   0.0343
#> 
#> Mean Group Variables: log_hc, log_ck, log_ngd, lag1_log_rgdpo
#> Cross Sectional Averaged Variables: log_rgdpo, log_hc, log_ck, log_ngd (lags=3)
#> 
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

CS-ARDL

fit_cs_ardl <- csdm(
  form, data = dat, id = "id", time = "year", model = "cs_ardl",
  csa = dynamic_csa,
  lr = ardl_1_0
)

summary(fit_cs_ardl)
#> csdm summary: Cross-Sectional ARDL (CS-ARDL)
#> Formula: log_rgdpo ~ log_hc + log_ck + log_ngd
#> N: 15, T: 38
#> Number of obs: 525
#> R-squared (mg): 0.9861
#> 
#> CD = -2.392, p = 0.0168
#> (For additional CD diagnostics, use cd_test())
#> 
#> Levels ARDL Est.
#>                  Coef. Std. Err.       z  P>|z| Signif. CI 2.5% CI 97.5%
#> (Intercept)     9.2888    8.0386  1.1555 0.2479         -6.4666  25.0441
#> log_hc         -1.9558    1.5659 -1.2490 0.2117         -5.0249   1.1132
#> log_ck          0.6666    0.2424  2.7499 0.0060      **  0.1915   1.1417
#> log_ngd        -0.3178    1.1271 -0.2819 0.7780         -2.5268   1.8913
#> lag1_log_rgdpo -0.0745    0.0555 -1.3424 0.1795         -0.1834   0.0343
#> 
#> Adjust. Term
#>                Coef. Std. Err.        z P>|z| Signif. CI 2.5% CI 97.5%
#> lr_log_rgdpo -1.0745    0.0555 -19.3513     0     *** -1.1834  -0.9657
#> 
#> Long Run Est.
#>              Coef. Std. Err.       z  P>|z| Signif. CI 2.5% CI 97.5% n_used
#> lr_log_hc  -1.8378    1.4743 -1.2465 0.2126         -4.7274   1.0518     15
#> lr_log_ck   0.6106    0.2076  2.9409 0.0033      **  0.2037   1.0175     15
#> lr_log_ngd -0.4098    1.2364 -0.3314 0.7403         -2.8330   2.0134     15
#> 
#> Mean Group Variables: lag1_log_rgdpo, log_hc, log_ck, log_ngd
#> Cross Sectional Averaged Variables: log_rgdpo, log_hc, log_ck, log_ngd (lags=3)
#> Long Run Variables: log_hc, log_ck, log_ngd
#> Cointegration variable(s): log_rgdpo
#> 
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coef(fit_cs_ardl, component = "long_run")
#>  lr_log_hc  lr_log_ck lr_log_ngd 
#> -1.8377835  0.6105883 -0.4097798

The four calls use the same levels formula, so their coefficient meanings are directly comparable only after accounting for their different dynamic and CCE terms. In particular, a contemporaneous coefficient in DCCE is not a long-run effect.

Samples, missing values, and identification

csdm() requires unique, nonmissing unit-time keys. By default, numeric time indexes advance in steps of one; set time_step when the intended grid differs. Lags respect gaps rather than treating the previous observed row as the previous period.

subset is evaluated before model construction. The supported missing-value policies are na.omit, na.exclude, and na.fail. Unit regressions must retain positive residual degrees of freedom, and every economic coefficient must be identified after projecting out the CCE terms. Units that fail these checks are listed with reasons in fit$units. At least two eligible units are required.

For csa = csdm_csa("_all"), the source sample for cross-sectional averages is the complete base formula sample after subsetting and before dynamic lag trimming. It therefore need not equal the final set of fitted observations. Set fullsample = TRUE to calculate each average from all finite observations of that variable in the selected sample. This is useful when the averaging variables have different missing-value patterns and mirrors the fullsample option used in xtdcce2.

Inference and R model methods

Standard methods expose the fitted model and its sample:

coef(fit_cce)
#> (Intercept)      log_hc      log_ck     log_ngd 
#>   1.9002762  -1.4921059   0.1366750   0.8075283
sqrt(diag(vcov(fit_cce)))
#> (Intercept)      log_hc      log_ck     log_ngd 
#>  2.11954323  1.01523715  0.09559111  0.29715332
nobs(fit_cce)
#> [1] 570
head(model.frame(fit_cce))
#>   log_rgdpo    log_hc   log_ck   log_ngd
#> 1  8.016315 0.7679234 11.79632 -2.727981
#> 2  8.035557 0.7790388 11.84772 -2.712420
#> 3  8.035343 0.7901542 11.90875 -2.702157
#> 4  8.060844 0.8012696 11.97760 -2.697793
#> 5  8.075463 0.8123851 12.03988 -2.701199
#> 6  8.034805 0.8235005 12.07300 -2.709801

head(residuals(fit_cce, format = "long"))
#>    id year .row    residual
#> 11  1 1970    1  0.01979243
#> 12  1 1971    2 -0.00607203
#> 13  1 1972    3 -0.03194126
#> 14  1 1973    4 -0.03727097
#> 15  1 1974    5  0.00548517
#> 16  1 1975    6  0.02278676
head(fitted(fit_cce, format = "vector"))
#>       11       12       13       14       15       16 
#> 7.996523 8.041629 8.067284 8.098115 8.069978 8.012019

For CS-ARDL, coef() and vcov() accept component = "levels", "adjustment", "long_run", or "all". Every component uses a common eligible-unit sample for its coefficients and covariance. Exact algebraic relationships can make the combined "all" covariance singular.

The package also supplies tidy(), glance(), and augment() methods:

broom::tidy(fit_cce, conf.int = TRUE)
broom::glance(fit_cce)
broom::augment(fit_cce)

modelsummary::modelsummary(
  list(MG = fit_mg, CCE = fit_cce, DCCE = fit_dcce),
  statistic = "std.error"
)

Residual cross-sectional dependence

Let $e_{it}$ denote the fitted residual and $T_{ij}$ the number of overlapping finite observations for units $i$ and $j$ . The classical statistic implemented by cd_test() is

$CD=\sqrt{\frac{2}{N(N-1)}} \sum_{i<j}\sqrt{T_{ij}}\,\widehat\rho_{ij}.$

It uses pairwise-complete correlations by default. Large absolute values are evidence against the null represented by the standard-normal reference approximation. Depending on the theoretical formulation, that null is stated as cross-sectional independence or sufficiently weak dependence (Pesaran, 2015, 2021).

For a balanced residual matrix, CDw draws one independent Rademacher weight $w_i\in\{-1,1\}$ per unit and computes (Juodis & Reese, 2021)

$CD_W= \left(\frac{1}{NT}\sum_{i,t}w_i^2\widetilde e_{it}^2\right)^{-1} \sqrt{\frac{2}{TN(N-1)}} \sum_t\sum_{i<j}w_i\widetilde e_{it}w_j\widetilde e_{jt},$

where $\widetilde e_{it}$ is demeaned within unit. Set seed to reproduce the random weights without changing the caller’s random-number state.

CDw+ adds the power-enhancement screening term (Fan et al., 2015):

$CD_{W+}=CD_W+ \sum_{i<j}|\widehat\rho_{ij}|\, 1\left\{|\widehat\rho_{ij}|>2\sqrt{\log(N)/T}\right\}.$

The threshold is applied to the ordinary residual correlation, without a $\sqrt{T}$ multiplier. The enhancement is nonnegative, so CDw+ is not a second independent random-sign test.

CD* applies the Pesaran-Xie bias correction after removing n_pc principal components (Pesaran & Xie, 2021). Its approximation requires a nondegenerate bias-correction denominator; a numerically computable result alone does not establish that the asymptotic assumptions are suitable.

cd_test(fit_mg, type = "CD")
#> Cross-sectional dependence tests
#> N = 15, T = 38
#> 
#>    statistic p.value
#> CD     3.238   0.001
cd_test(fit_cce, type = "all", seed = 42)
#> Cross-sectional dependence tests
#> N = 15, T = 38
#> 
#>        statistic p.value
#> CD        -2.676   0.007
#> CDw       -0.714   0.475
#> CDw+       2.247   0.025
#> CDstar    -2.830   0.005

Periods with no finite residual for any retained unit are outside the effective residual sample and are removed automatically. Partially observed periods are handled as follows:

classical CD remains pairwise-complete under the default na.action = "pairwise";
CDw and CDw+ require a balanced matrix and otherwise error;
CD* requires a balanced matrix and otherwise returns NA with a warning.

To evaluate all diagnostics on one common time sample, use:

cd_test(
  fit_cce,
  type = "all",
  seed = 42,
  na.action = "drop.incomplete.times"
)

A rejection concerns the residual dependence targeted by the selected test; it does not by itself identify the source of misspecification. A non-rejection does not prove residual independence, particularly in small samples or weak-power settings.

Practical limits

Large- $N$ normal inference does not correct short- $T$ dynamic bias.
CCE validity requires enough informative averages and suitable factor and loading conditions.
Long-run ratios require a stable, nonzero AR denominator and an economically defensible long-run interpretation.
Reported mean-group fit statistics summarize unit regressions; they are not pooled-regression goodness-of-fit measures.
Saved fits from older releases should be refitted because corrected samples, transformations, and covariance calculations can change results.

References

Chudik, A., & Pesaran, M. H. (2015). Common correlated effects estimation of heterogeneous dynamic panel data models with weakly exogenous regressors. Journal of Econometrics, 188(2), 393–420.

Ditzen, J. (2018). Estimating dynamic common-correlated effects in STATA. The STATA Journal, 18(3), 585–617. https://doi.org/10.1177/1536867X1801800306

Fan, J., Liao, Y., & Yao, J. (2015). Power enhancement in high-dimensional cross-section tests. Econometrica, 83(4), 1497–1541.

Juodis, A., & Reese, S. (2021). The incidental parameters problem in testing for remaining cross-sectional correlation. Journal of Business and Economic Statistics, 40(3), 1191–1203.

Pesaran, M. H. (2006). Estimation and inference in large heterogeneous panels with multifactor error structure. Econometrica, 74(4), 967–1012.

Pesaran, M. H. (2015). Testing weak cross-sectional dependence in large panels. Econometric Reviews, 34(6-10), 1089–1117.

Pesaran, M. H. (2021). General diagnostic tests for cross-sectional dependence in panels. Empirical Economics, 60(1), 13–50.

Pesaran, M. H., & Smith, R. (1995). Estimating long-run relationships from dynamic heterogeneous panels. Journal of Econometrics, 68(1), 79–113.

Pesaran, M. H., & Xie, Y. (2021). A bias-corrected CD test for error cross-sectional dependence in panel models. Econometric Reviews, 41(6), 649–677.

Joao Claudio Macosso