Estimation interface

Package plm provides various functions for panel data estimation, among them:

• plm: estimation of the basic panel models and instrumental variable panel models, i.e., between and first-difference models and within and random effect models. Models are estimated internally using the lm function on transformed data,
• pvcm: estimation of models with variable coefficients,
• pgmm: estimation of generalized method of moments models,
• pggls: estimation of general feasible generalized least squares models,
• pmg: estimators for mean groups (MG), demeaned MG (DMG) and common correlated effects MG (CCEMG) for heterogeneous panel models,
• pcce: estimators for common correlated effects mean groups (CCEMG) and pooled (CCEP) for panel data with common factors,
• pldv: panel estimators for limited dependent variables.

The interface of these functions is consistent with the lm() function. Namely, their first two arguments are formula and data (which should be a data.frame and is mandatory). Three additional arguments are common to these functions:

• index: this argument enables the estimation functions to identify the structure of the data, i.e., the individual and the time period for each observation,
• effect: the kind of effects to include in the model, i.e., individual effects, time effects or both³,
• model: the kind of model to be estimated, most of the time a model with fixed effects or a model with random effects.

The results of these four functions are stored in an object which class has the same name of the function. They all inherit from class panelmodel. A panelmodel object contains: coefficients, residuals, fitted.values, vcov, df.residual and call and functions that extract these elements are provided.

Testing interface

The diagnostic testing interface provides both formula and panelmodel methods for most functions, with some exceptions. The user may thus choose whether to employ results stored in a previously estimated panelmodel object or to re-estimate it for the sake of testing.

Although the first strategy is the most efficient one, diagnostic testing on panel models mostly employs OLS residuals from pooling model objects, whose estimation is computationally inexpensive. Therefore most examples in the following are based on formula methods, which are perhaps the cleanest for illustrative purposes.

The (quasi–)demeaning framework

The estimation methods for the basic models in panel data econometrics, the pooled OLS, random effects and fixed effects (or within) models, can all be described inside the OLS estimation framework. In fact, while pooled OLS simply pools data, the standard way of estimating fixed effects models with, say, group (time) effects entails transforming the data by subtracting the average over time (group) to every variable, which is usually termed time-demeaning. In the random effects case, the various feasible GLS estimators which have been put forth to tackle the issue of serial correlation induced by the group-invariant random effect have been proven to be equivalent (as far as estimation of \(\beta\)s is concerned) to OLS on partially demeaned data, where partial demeaning is defined as:

\[\begin{equation} (\#eq:ldemmodel) y_{it} - \theta \bar{y}_i = ( X_{it} - \theta \bar{X}_{i} ) \beta + ( u_{it} - \theta \bar{u}_i ) \end{equation}\]

where \(\theta=1-[\sigma_u^2 / (\sigma_u^2 + T \sigma_e^2)]^{1/2}\), \(\bar{y}\) and \(\bar{X}\) denote time means of \(y\) and \(X\), and the disturbance \(v_{it} - \theta \bar{v}_i\) is homoskedastic and serially uncorrelated. Thus the feasible RE estimate for \(\beta\) may be obtained estimating \(\hat{\theta}\) and running an OLS regression on the transformed data with lm(). The other estimators can be computed as special cases: for \(\theta=1\) one gets the fixed effects estimator, for \(\theta=0\) the pooled OLS one.

Moreover, instrumental variable estimators of all these models may also be obtained using several calls to lm().

For this reason the three above estimators have been grouped inside the same function.

On the output side, a number of diagnostics and a very general coefficients’ covariance matrix estimator also benefits from this framework, as they can be readily calculated applying the standard OLS formulas to the demeaned data, which are contained inside plm objects. This will be the subject of subsection inference in the panel model.

The object oriented approach to general GLS computations

The covariance matrix of errors in general GLS models is too generic to fit the quasi-demeaning framework, so this method calls for a full-blown application of GLS as in @ref(eq:naive). On the other hand, this estimator relies heavily on \(n\)–asymptotics, making it theoretically most suitable for situations which forbid it computationally: e.g., “short” micropanels with thousands of individuals observed over few time periods.

R has general facilities for fast matrix computation based on object orientation: particular types of matrices (symmetric, sparse, dense etc.) are assigned the relevant class and the additional information on structure is used in the computations, sometimes with dramatic effects on performance (see Bates (2004)) and packages Matrix (see Bates and Maechler (2016)) and SparseM (see Koenker and Ng (2016)). Some optimized linear algebra routines are available in the R package bdsmatrix (see Therneau (2014)) which exploit the particular block-diagonal and symmetric structure of \(\hat{V}\) making it possible to implement a fast and reliable full-matrix solution to problems of any practically relevant size.

The \(\hat{V}\) matrix is constructed as an object of class bdsmatrix. The peculiar properties of this matrix class are used for efficiently storing the object in memory and then by ad-hoc versions of the solve and crossprod methods, dramatically reducing computing times and memory usage. The resulting matrix is then used “the naive way” as in @ref(eq:naive) to compute \(\hat{\beta}\), resulting in speed comparable to that of the demeaning solution.

Random effects estimators

As observed above, the random effect model is obtained as a linear estimation on quasi-demeaned data. The parameter of this transformation is obtained using preliminary estimations.

Four estimators of this parameter are available, depending on the value of the argument random.method:

• "swar": from Swamy and Arora (1972), the default value,
• "walhus": from Wallace and Hussain (1969),
• "amemiya": from Amemiya (1971),
• "nerlove": from Nerlove (1971).
• "ht": for Hausman-Taylor-type instrumental variable (IV) estimation, discussed later, see Section Instrumental variable estimator.

For example, to use the amemiya estimator:

grun.amem <- plm(inv~value+capital, data=Grunfeld,
                 model="random", random.method="amemiya")

The estimation of the variance of the error components are performed using the ercomp function, which has a method and an effect argument, and can be used by itself:

ercomp(inv~value+capital, data=Grunfeld, method = "amemiya", effect = "twoways")

##                   var std.dev share
## idiosyncratic 2644.13   51.42 0.256
## individual    7452.02   86.33 0.721
## time           243.78   15.61 0.024
## theta: 0.868 (id) 0.2787 (time) 0.2776 (total)

Introducing time or two-ways effects

The default behavior of plm is to introduce individual effects. Using the effect argument, one may also introduce:

• time effects (effect = "time"),
• individual and time effects (effect = "twoways").

For example, to estimate a two-ways effect model for the Grunfeld data:

grun.tways <- plm(inv~value+capital, data = Grunfeld, effect = "twoways",
                  model = "random", random.method = "amemiya")
summary(grun.tways)

## Twoways effects Random Effect Model 
##    (Amemiya's transformation)
## 
## Call:
## plm(formula = inv ~ value + capital, data = Grunfeld, effect = "twoways", 
##     model = "random", random.method = "amemiya")
## 
## Balanced Panel: n = 10, T = 20, N = 200
## 
## Effects:
##                   var std.dev share
## idiosyncratic 2644.13   51.42 0.256
## individual    7452.02   86.33 0.721
## time           243.78   15.61 0.024
## theta: 0.868 (id) 0.2787 (time) 0.2776 (total)
## 
## Residuals:
##      Min.   1st Qu.    Median   3rd Qu.      Max. 
## -176.9062  -18.0431    3.2697   17.1719  234.1735 
## 
## Coefficients:
##               Estimate Std. Error z-value             Pr(>|z|)    
## (Intercept) -63.767791  29.851537 -2.1362              0.03267 *  
## value         0.111386   0.010909 10.2102 < 0.0000000000000002 ***
## capital       0.323321   0.018772 17.2232 < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    2066800
## Residual Sum of Squares: 518200
## R-Squared:      0.74927
## Adj. R-Squared: 0.74673
## Chisq: 588.717 on 2 DF, p-value: < 0.000000000000000222

In the “effects” section of the printed summary of the result, the variance of the three elements of the error term and the three parameters used in the transformation are printed.

Unbalanced panels

Estimations by plm support unbalanced panel models.

The following example is using data used by Harrison and Rubinfeld (1978) to estimate an hedonic housing prices function. It is reproduced in Baltagi and Chang (1994), table 2 (and in Baltagi (2005), pp. 172/4; Baltagi (2013), pp. 195/7 tables 9.1/3).

data("Hedonic", package = "plm")
Hed <- plm(mv~crim+zn+indus+chas+nox+rm+age+dis+rad+tax+ptratio+blacks+lstat,
           data = Hedonic, model = "random", index = "townid")
summary(Hed)

## Oneway (individual) effect Random Effect Model 
##    (Swamy-Arora's transformation)
## 
## Call:
## plm(formula = mv ~ crim + zn + indus + chas + nox + rm + age + 
##     dis + rad + tax + ptratio + blacks + lstat, data = Hedonic, 
##     model = "random", index = "townid")
## 
## Unbalanced Panel: n = 92, T = 1-30, N = 506
## 
## Effects:
##                   var std.dev share
## idiosyncratic 0.01696 0.13025 0.562
## individual    0.01324 0.11505 0.438
## theta:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.2505  0.5483  0.6284  0.6141  0.7147  0.7976 
## 
## Residuals:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -0.62902 -0.06712 -0.00156 -0.00216  0.06858  0.54973 
## 
## Coefficients:
##                 Estimate   Std. Error  z-value              Pr(>|z|)    
## (Intercept)  9.685866695  0.197510264  49.0398 < 0.00000000000000022 ***
## crim        -0.007411967  0.001047812  -7.0738   0.00000000000150795 ***
## zn           0.000078877  0.000650012   0.1213             0.9034166    
## indus        0.001556340  0.004034911   0.3857             0.6997051    
## chasyes     -0.004424737  0.029211764  -0.1515             0.8796041    
## nox         -0.005842506  0.001245183  -4.6921   0.00000270431168602 ***
## rm           0.009055167  0.001188629   7.6182   0.00000000000002573 ***
## age         -0.000857873  0.000467933  -1.8333             0.0667541 .  
## dis         -0.144418433  0.044093739  -3.2753             0.0010557 ** 
## rad          0.095983935  0.026610945   3.6069             0.0003098 ***
## tax         -0.000377396  0.000176926  -2.1331             0.0329190 *  
## ptratio     -0.029475776  0.009069842  -3.2499             0.0011546 ** 
## blacks       0.562775469  0.101973789   5.5188   0.00000003412743874 ***
## lstat       -0.291074917  0.023927306 -12.1650 < 0.00000000000000022 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    987.94
## Residual Sum of Squares: 8.9988
## R-Squared:      0.99091
## Adj. R-Squared: 0.99067
## Chisq: 1199.5 on 13 DF, p-value: < 0.000000000000000222

Measures for the unbalancedness of a panel data set or the data used in estimated models are provided by function punbalancedness. It gives the measures \(\gamma\) and \(\nu\) from Ahrens and Pincus (1981) where for both 1 represents balanced data and the more unbalanced the data the lower the value.

punbalancedness(Hed)

##     gamma        nu 
## 0.4715336 0.5188292

Instrumental variable estimators

All of the models presented above may be estimated using instrumental variables. The instruments are specified at the end of the formula after a | sign (pipe).

The instrumental variables estimator used is indicated with the inst.method argument:

• "bvk", from Balestra and Varadharajan–Krishnakumar (1987), the default value,
• "baltagi", from Baltagi (1981),
• "am", from Amemiya and MaCurdy (1986),
• "bms", from Breusch, Mizon, and Schmidt (1989).

An illustration is in the following example from Baltagi (2005), p. 120; Baltagi (2013), p. 137; Baltagi (2021), p. 165, table 7.3 (“G2SLS”).

data("Crime", package = "plm")
cr <- plm(lcrmrte ~ lprbarr + lpolpc + lprbconv + lprbpris + lavgsen +
          ldensity + lwcon + lwtuc + lwtrd + lwfir + lwser + lwmfg + lwfed +
          lwsta + lwloc + lpctymle + lpctmin + region + smsa + factor(year)
          | . - lprbarr - lpolpc + ltaxpc + lmix,
          data = Crime, model = "random")
summary(cr)

## Oneway (individual) effect Random Effect Model 
##    (Swamy-Arora's transformation)
## Instrumental variable estimation
##    (Balestra-Varadharajan-Krishnakumar's transformation)
## 
## Call:
## plm(formula = lcrmrte ~ lprbarr + lpolpc + lprbconv + lprbpris + 
##     lavgsen + ldensity + lwcon + lwtuc + lwtrd + lwfir + lwser + 
##     lwmfg + lwfed + lwsta + lwloc + lpctymle + lpctmin + region + 
##     smsa + factor(year) | . - lprbarr - lpolpc + ltaxpc + lmix, 
##     data = Crime, model = "random")
## 
## Balanced Panel: n = 90, T = 7, N = 630
## 
## Effects:
##                   var std.dev share
## idiosyncratic 0.02227 0.14924 0.326
## individual    0.04604 0.21456 0.674
## theta: 0.7457
## 
## Residuals:
##       Min.    1st Qu.     Median    3rd Qu.       Max. 
## -0.7485357 -0.0709883  0.0040648  0.0784455  0.4756273 
## 
## Coefficients:
##                  Estimate Std. Error z-value      Pr(>|z|)    
## (Intercept)    -0.4538501  1.7029831 -0.2665      0.789852    
## lprbarr        -0.4141383  0.2210496 -1.8735      0.060998 .  
## lpolpc          0.5049461  0.2277778  2.2168      0.026634 *  
## lprbconv       -0.3432506  0.1324648 -2.5913      0.009563 ** 
## lprbpris       -0.1900467  0.0733392 -2.5913      0.009560 ** 
## lavgsen        -0.0064389  0.0289407 -0.2225      0.823935    
## ldensity        0.4343449  0.0711496  6.1047 0.00000000103 ***
## lwcon          -0.0042958  0.0414226 -0.1037      0.917403    
## lwtuc           0.0444589  0.0215448  2.0636      0.039060 *  
## lwtrd          -0.0085579  0.0419829 -0.2038      0.838476    
## lwfir          -0.0040305  0.0294569 -0.1368      0.891166    
## lwser           0.0105602  0.0215823  0.4893      0.624630    
## lwmfg          -0.2018020  0.0839373 -2.4042      0.016208 *  
## lwfed          -0.2134579  0.2151046 -0.9923      0.321029    
## lwsta          -0.0601232  0.1203149 -0.4997      0.617275    
## lwloc           0.1835363  0.1396775  1.3140      0.188846    
## lpctymle       -0.1458703  0.2268086 -0.6431      0.520131    
## lpctmin         0.1948763  0.0459385  4.2421 0.00002214292 ***
## regionwest     -0.2281821  0.1010260 -2.2586      0.023905 *  
## regioncentral  -0.1987703  0.0607475 -3.2721      0.001068 ** 
## smsayes        -0.2595451  0.1499718 -1.7306      0.083518 .  
## factor(year)82  0.0132147  0.0299924  0.4406      0.659500    
## factor(year)83 -0.0847693  0.0320010 -2.6490      0.008074 ** 
## factor(year)84 -0.1062027  0.0387893 -2.7379      0.006183 ** 
## factor(year)85 -0.0977457  0.0511681 -1.9103      0.056097 .  
## factor(year)86 -0.0719451  0.0605819 -1.1876      0.235004    
## factor(year)87 -0.0396595  0.0758531 -0.5228      0.601081    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    30.169
## Residual Sum of Squares: 12.419
## R-Squared:      0.5923
## Adj. R-Squared: 0.57472
## Chisq: 542.478 on 26 DF, p-value: < 0.000000000000000222

The Hausman-Taylor model (see Hausman and Taylor (1981)) may be estimated with the plm⁶ function by setting parameters random.method = "ht" and inst.method = "baltagi" like in the example below. The following replicates Baltagi (2005), pp. 129/30; Baltagi (2013), pp. 145/6, tables 7.4/5; Baltagi (2021), pp. 174/5 tables 7.5/6:

ht <- plm(lwage ~ wks + south + smsa + married + exp + I(exp ^ 2) + 
              bluecol + ind + union + sex + black + ed |
              bluecol + south + smsa + ind + sex + black |
              wks + married + union + exp + I(exp ^ 2), 
          data = Wages, index = 595,
          model = "random", random.method = "ht", inst.method = "baltagi")
summary(ht)

## Oneway (individual) effect Random Effect Model 
##    (Hausman-Taylor's transformation)
## Instrumental variable estimation
##    (Baltagi's transformation)
## 
## Call:
## plm(formula = lwage ~ wks + south + smsa + married + exp + I(exp^2) + 
##     bluecol + ind + union + sex + black + ed | bluecol + south + 
##     smsa + ind + sex + black | wks + married + union + exp + 
##     I(exp^2), data = Wages, model = "random", random.method = "ht", 
##     inst.method = "baltagi", index = 595)
## 
## Balanced Panel: n = 595, T = 7, N = 4165
## 
## Effects:
##                   var std.dev share
## idiosyncratic 0.02304 0.15180 0.025
## individual    0.88699 0.94180 0.975
## theta: 0.9392
## 
## Residuals:
##       Min.    1st Qu.     Median    3rd Qu.       Max. 
## -12.643736  -0.466002   0.043285   0.524739  13.340263 
## 
## Coefficients:
##                 Estimate   Std. Error z-value              Pr(>|z|)    
## (Intercept)  2.912726279  0.283652215 10.2687 < 0.00000000000000022 ***
## wks          0.000837403  0.000599732  1.3963               0.16263    
## southyes     0.007439837  0.031955005  0.2328               0.81590    
## smsayes     -0.041833367  0.018958129 -2.2066               0.02734 *  
## marriedyes  -0.029850749  0.018979963 -1.5728               0.11578    
## exp          0.113132791  0.002470954 45.7851 < 0.00000000000000022 ***
## I(exp^2)    -0.000418865  0.000054598 -7.6718   0.00000000000001696 ***
## bluecolyes  -0.020704707  0.013780948 -1.5024               0.13299    
## ind          0.013603930  0.015237366  0.8928               0.37196    
## unionyes     0.032771447  0.014908437  2.1982               0.02794 *  
## sexfemale   -0.130923610  0.126658988 -1.0337               0.30129    
## blackyes    -0.285747871  0.155701854 -1.8352               0.06647 .  
## ed           0.137943957  0.021248489  6.4919   0.00000000008473689 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    243.04
## Residual Sum of Squares: 4163.6
## R-Squared:      0.60945
## Adj. R-Squared: 0.60833
## Chisq: 6891.87 on 12 DF, p-value: < 0.000000000000000222

Unobserved effects test

The unobserved effects test à la Wooldridge (see Wooldridge (2002) 10.4.4), is a semiparametric test for the null hypothesis that \(\sigma^2_{\mu}=0\), i.e. that there are no unobserved effects in the residuals. Given that under the null the covariance matrix of the residuals for each individual is diagonal, the test statistic is based on the average of elements in the upper (or lower) triangle of its estimate, diagonal excluded: \(n^{-1/2} \sum_{i=1}^n \sum_{t=1}^{T-1} \sum_{s=t+1}^T \hat{u}_{it} \hat{u}_{is}\) (where \(\hat{u}\) are the pooled OLS residuals), which must be “statistically close” to zero under the null, scaled by its standard deviation: \[W=\frac{ \sum_{i=1}^n \sum_{t=1}^{T-1} \sum_{s=t+1}^T \hat{u}_{it} \hat{u}_{is} }{ [{ \sum_{i=1}^n ( \sum_{t=1}^{T-1} \sum_{s=t+1}^T \hat{u}_{it} \hat{u}_{is} } )^2 ]^{1/2} }\]

This test is (\(n\)-) asymptotically distributed as a standard normal regardless of the distribution of the errors. It does also not rely on homoskedasticity.

It has power both against the standard random effects specification, where the unobserved effects are constant within every group, as well as against any kind of serial correlation. As such, it “nests” both random effects and serial correlation tests, trading some power against more specific alternatives in exchange for robustness.

While not rejecting the null favours the use of pooled OLS, rejection may follow from serial correlation of different kinds, and in particular, quoting Wooldridge (2002), “should not be interpreted as implying that the random effects error structure must be true”.

Below, the test is applied to the data and model in Munnell (1990):

pwtest(log(gsp)~log(pcap)+log(pc)+log(emp)+unemp, data=Produc)

## 
##  Wooldridge's test for unobserved individual effects
## 
## data:  formula
## z = 3.9383, p-value = 0.00008207
## alternative hypothesis: unobserved effect

Locally robust tests for serial correlation or random effects

The presence of random effects may affect tests for residual serial correlation, and the opposite. One solution is to use a joint test, which has power against both alternatives. A joint LM test for random effects and serial correlation under normality and homoskedasticity of the idiosyncratic errors has been derived by Baltagi and Li (1991) and Baltagi and Li (1995) and is implemented as an option in pbsytest:

pbsytest(log(gsp)~log(pcap)+log(pc)+log(emp)+unemp, data=Produc, test="j")

## 
##  Baltagi and Li AR-RE joint test - balanced panel
## 
## data:  formula
## chisq = 4187.6, df = 2, p-value < 0.00000000000000022
## alternative hypothesis: AR(1) errors or random effects

Rejection of the joint test, though, gives no information on the direction of the departure from the null hypothesis, i.e.: is rejection due to the presence of serial correlation, of random effects or of both?

Bera, Sosa–Escudero, and Yoon (2001) (hereafter BSY) derive locally robust tests both for individual random effects and for first-order serial correlation in residuals as “corrected” versions of the standard LM test (see plmtest). While still dependent on normality and homoskedasticity, these are robust to local departures from the hypotheses of, respectively, no serial correlation or no random effects. The authors observe that, although suboptimal, these tests may help detecting the right direction of the departure from the null, thus complementing the use of joint tests. Moreover, being based on pooled OLS residuals, the BSY tests are computationally far less demanding than likelihood-based conditional tests.

On the other hand, the statistical properties of these “locally corrected” tests are inferior to those of the non-corrected counterparts when the latter are correctly specified. If there is no serial correlation, then the optimal test for random effects is the likelihood-based LM test of Breusch and Godfrey (with refinements by Honda, see plmtest), while if there are no random effects the optimal test for serial correlation is, again, Breusch-Godfrey’s test¹². If the presence of a random effect is taken for granted, then the optimal test for serial correlation is the likelihood-based conditional LM test of Baltagi and Li (1995) (see pbltest).

The serial correlation version is the default:

pbsytest(log(gsp)~log(pcap)+log(pc)+log(emp)+unemp, data=Produc)

## 
##  Bera, Sosa-Escudero and Yoon locally robust test - balanced panel
## 
## data:  formula
## chisq = 52.636, df = 1, p-value = 0.0000000000004015
## alternative hypothesis: AR(1) errors sub random effects

The BSY test for random effects is implemented in the one-sided version¹³, which takes heed that the variance of the random effect must be non-negative:

pbsytest(log(gsp)~log(pcap)+log(pc)+log(emp)+unemp, data=Produc, test="re")

## 
##  Bera, Sosa-Escudero and Yoon locally robust test (one-sided) -
##  balanced panel
## 
## data:  formula
## z = 57.914, p-value < 0.00000000000000022
## alternative hypothesis: random effects sub AR(1) errors

Conditional LM test for AR(1) or MA(1) errors under random effects

Baltagi and Li (1991) and Baltagi and Li (1995) derive a Lagrange multiplier test for serial correlation in the idiosyncratic component of the errors under (normal, heteroskedastic) random effects. Under the null of serially uncorrelated errors, the test turns out to be identical for both the alternative of AR(1) and MA(1) processes. One- and two-sided versions are provided, the one-sided having power against positive serial correlation only. The two-sided is the default, while for the other one must specify the alternative option to "onesided":

pbltest(log(gsp)~log(pcap)+log(pc)+log(emp)+unemp, 
        data=Produc, alternative="onesided")

## 
##  Baltagi and Li one-sided LM test
## 
## data:  log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp
## z = 21.69, p-value < 0.00000000000000022
## alternative hypothesis: AR(1)/MA(1) errors in RE panel model

As usual, the LM test statistic is based on residuals from the maximum likelihood estimate of the restricted model (random effects with serially uncorrelated errors). In this case, though, the restricted model cannot be estimated by OLS anymore, therefore the testing function depends on lme() in the nlme package for estimation of a random effects model by maximum likelihood. For this reason, the test is applicable only to balanced panels.

No test has been implemented to date for the symmetric hypothesis of no random effects in a model with errors following an AR(1) process, but an asymptotically equivalent likelihood ratio test is available in the nlme package (see Section plm versus nlme and lme4).

General serial correlation tests

A general testing procedure for serial correlation in fixed effects (FE), random effects (RE) and pooled-OLS panel models alike can be based on considerations in Wooldridge (2002), 10.7.2.

Recall that plm model objects are the result of OLS estimation performed on “demeaned” data, where, in the case of individual effects (else symmetric), this means time-demeaning for the FE (within) model, quasi-time-demeaning for the RE (random) model and original data, with no demeaning at all, for the pooled OLS (pooling) model (see Section software approach).

For the random effects model, Wooldridge (2002) observes that under the null of homoskedasticity and no serial correlation in the idiosyncratic errors, the residuals from the quasi-demeaned regression must be spherical as well. Else, as the individual effects are wiped out in the demeaning, any remaining serial correlation must be due to the idiosyncratic component. Hence, a simple way of testing for serial correlation is to apply a standard serial correlation test to the quasi-demeaned model. The same applies in a pooled model, w.r.t. the original data.

The FE case needs some qualification. It is well-known that if the original model’s errors are uncorrelated then FE residuals are negatively serially correlated, with \(cor(\hat{u}_{it}, \hat{u}_{is})=-1/(T-1)\) for each \(t,s\) (see Wooldridge (2002) 10.5.4). This correlation clearly dies out as T increases, so this kind of AR test is applicable to within model objects only for T “sufficiently large”¹⁴. On the converse, in short panels the test gets severely biased towards rejection (or, as the induced correlation is negative, towards acceptance in the case of the one-sided DW test with alternative="greater"). See below for a serial correlation test applicable to “short” FE panel models.

plm objects retain the “demeaned” data, so the procedure is straightforward for them. The wrapper functions pbgtest and pdwtest re-estimate the relevant quasi-demeaned model by OLS and apply, respectively, standard Breusch-Godfrey and Durbin-Watson tests from package lmtest:

pbgtest(grun.fe, order = 2)

## 
##  Breusch-Godfrey/Wooldridge test for serial correlation in panel models
## 
## data:  inv ~ value + capital
## chisq = 42.587, df = 2, p-value = 0.0000000005655
## alternative hypothesis: serial correlation in idiosyncratic errors

The tests share the features of their OLS counterparts, in particular the pbgtest allows testing for higher-order serial correlation, which might turn useful, e.g., on quarterly data. Analogously, from the point of view of software, as the functions are simple wrappers towards bgtest and dwtest, all arguments from the latter two apply and may be passed on through the ellipsis (the ... argument).

Wooldridge’s test for serial correlation in “short” FE panels

For the reasons reported above, under the null of no serial correlation in the errors, the residuals of a FE model must be negatively serially correlated, with \(cor(\hat{\epsilon}_{it}, \hat{\epsilon}_{is})=-1/(T-1)\) for each \(t,s\). Wooldridge suggests basing a test for this null hypothesis on a pooled regression of FE residuals on themselves, lagged one period: \[\hat{\epsilon}_{i,t}=\alpha + \delta \hat{\epsilon}_{i,t-1} + \eta_{i,t}\] Rejecting the restriction \(\delta = -1/(T-1)\) makes us conclude against the original null of no serial correlation.

The building blocks available in plm make it easy to construct a function carrying out this procedure: first the FE model is estimated and the residuals retrieved, then they are lagged and a pooling AR(1) model is estimated. The test statistic is obtained by applying the above restriction on \(\delta\) and supplying a heteroskedasticity- and autocorrelation-consistent covariance matrix (vcovHC with the appropriate options, in particular method="arellano")¹⁵.

pwartest(log(emp) ~ log(wage) + log(capital), data=EmplUK)

## 
##  Wooldridge's test for serial correlation in FE panels
## 
## data:  plm.model
## F = 312.3, df1 = 1, df2 = 889, p-value < 0.00000000000000022
## alternative hypothesis: serial correlation

The test is applicable to any FE panel model, and in particular to “short” panels with small \(T\) and large \(n\).

Wooldridge’s first-difference-based test

In the context of the first difference model, Wooldridge (2002), 10.6.3 proposes a serial correlation test that can also be seen as a specification test to choose the most efficient estimator between fixed effects (within) and first difference (fd).

The starting point is the observation that if the idiosyncratic errors of the original model \(u_{it}\) are uncorrelated, the errors of the (first) differenced model¹⁶ \(e_{it} \equiv u_{it}-u_{i,t-1}\) will be correlated, with \(cor(e_{it}, e_{i,t-1})=-0.5\), while any time-invariant effect, “fixed” or “random”, is wiped out in the differencing. So a serial correlation test for models with individual effects of any kind can be based on estimating the model \[\hat{u}_{i,t}= \delta \hat{u}_{i,t-1} + \eta_{i,t}\] and testing the restriction \(\delta = -0.5\), corresponding to the null of no serial correlation. Drukker (2003) provides Monte Carlo evidence of the good empirical properties of the test.

On the other extreme (see Wooldridge (2002) 10.6.1), if the differenced errors \(e_{it}\) are uncorrelated, as by definition \(u_{it} = u_{i,t-1} + e_{it}\), then \(u_{it}\) is a random walk. In this latter case, the most efficient estimator is the first difference (fd) one; in the former case, it is the fixed effects one (within).

The function pwfdtest allows testing either hypothesis: the default behaviour h0="fd" is to test for serial correlation in first-differenced errors:

pwfdtest(log(emp) ~ log(wage) + log(capital), data=EmplUK)

## 
##  Wooldridge's first-difference test for serial correlation in panels
## 
## data:  plm.model
## F = 1.5251, df1 = 1, df2 = 749, p-value = 0.2172
## alternative hypothesis: serial correlation in differenced errors

while specifying h0="fe" the null hypothesis becomes no serial correlation in original errors, which is similar to the pwartest.

pwfdtest(log(emp) ~ log(wage) + log(capital), data=EmplUK, h0="fe")

## 
##  Wooldridge's first-difference test for serial correlation in panels
## 
## data:  plm.model
## F = 131.55, df1 = 1, df2 = 749, p-value < 0.00000000000000022
## alternative hypothesis: serial correlation in original errors

Not rejecting one of the two is evidence in favour of using the estimator corresponding to h0. Should the truth lie in the middle (both rejected), whichever estimator is chosen will have serially correlated errors: therefore it will be advisable to use the autocorrelation-robust covariance estimators from the subsection robust covariance matrix estimation in inference.

CD and LM-type tests for global cross-sectional dependence

The function pcdtest implements a family of XSD tests which can be applied in different settings, ranging from those where \(T\) grows large with \(n\) fixed to “short” panels with a big \(n\) dimension and a few time periods. All are based on (transformations of–) the product-moment correlation coefficient of a model’s residuals, defined as

\[ \hat{\rho}_{ij}=\frac{\sum_{t=1}^T \hat{u}_{it} \hat{u}_{jt}}{(\sum_{t=1}^T \hat{u}^2_{it})^{1/2} (\sum_{t=1}^T \hat{u}^2_{jt})^{1/2} } \]

i.e., as averages over the time dimension of pairwise correlation coefficients for each pair of cross-sectional units.

The Breusch-Pagan (Breusch and Pagan 1980) LM test, based on the squares of \(\rho_{ij}\), is valid for \(T \rightarrow \infty\) with \(n\) fixed; defined as

\[LM=\sum_{i=1}^{n-1} \sum_{j=i+1}^{n} T_{ij} \hat{\rho}_{ij}^2\]

where in the case of an unbalanced panel only pairwise complete observations are considered, and \(T_{ij}=min(T_i,T_j)\) with \(T_i\) being the number of observations for individual \(i\); else, if the panel is balanced, \(T_{ij}=T\) for each \(i,j\). The test is distributed as \(\chi^2_{n(n-1)/2}\). It is inappropriate whenever the \(n\) dimension is “large”. A scaled version, applicable also if \(T \rightarrow \infty\) and then \(n \rightarrow \infty\) (as in some pooled time series contexts), is defined as

\[SCLM=\sqrt{\frac{1}{n(n-1)}} ( \sum_{i=1}^{n-1} \sum_{j=i+1}^{n} T_{ij} \hat{\rho}_{ij}^2 -1 )\]

and distributed as a standard normal (see Pesaran (2004)).

A bias-corrected scaled version, \(BCSCLM\), for the fixed effect model with individual effects only is also available which is simply the \(SCLM\) with a term correcting for the bias (Baltagi, Feng, and Kao (2012))¹⁸. This statistic is also asymptotically distributed as standard normal. \[BCSCLM=\sqrt{\frac{1}{n(n-1)}} ( \sum_{i=1}^{n-1} \sum_{j=i+1}^{n} T_{ij} \hat{\rho}_{ij}^2 -1)-\frac{n}{2(T-1)}\]

Pesaran’s (Pesaran (2004), Pesaran (2015)) \(CD\) test

\[CD=\sqrt{\frac{2}{n(n-1)}} ( \sum_{i=1}^{n-1} \sum_{j=i+1}^{n} \sqrt{T_{ij}} \hat{\rho}_{ij} )\]

based on \(\rho_{ij}\) without squaring (also distributed as a standard normal) is appropriate both in \(n\)– and in \(T\)–asymptotic settings. It has remarkable properties in samples of any practically relevant size and is robust to a variety of settings. The only big drawback is that the test loses power against the alternative of cross-sectional dependence if the latter is due to a factor structure with factor loadings averaging zero, that is, some units react positively to common shocks, others negatively.

The default version of the test is "cd" yielding Pesaran’s \(CD\) test. These tests are originally meant to use the residuals of separate estimation of one time-series regression for each cross-sectional unit, so this is the default behaviour of pcdtest.

pcdtest(inv~value+capital, data=Grunfeld)

## 
##  Pesaran CD test for cross-sectional dependence in panels
## 
## data:  inv ~ value + capital
## z = 5.3401, p-value = 0.00000009292
## alternative hypothesis: cross-sectional dependence

If a different model specification (within, random, …) is assumed consistent, one can resort to its residuals for testing¹⁹ by specifying the relevant model type. The main argument of this function may be either a model of class panelmodel or a formula and a data.frame; in the second case, unless model is set to NULL, all usual parameters relative to the estimation of a plm model may be passed on. The test is compatible with any consistent panelmodel for the data at hand, with any specification of effect. E.g., specifying effect = "time" or effect = "twoways" allows to test for residual cross-sectional dependence after the introduction of time fixed effects to account for common shocks.

pcdtest(inv~value+capital, data=Grunfeld, model="within")

## 
##  Pesaran CD test for cross-sectional dependence in panels
## 
## data:  inv ~ value + capital
## z = 4.6612, p-value = 0.000003144
## alternative hypothesis: cross-sectional dependence

If the time dimension is insufficient and model=NULL, the function defaults to estimation of a within model and issues a warning.

CD(p) test for local cross-sectional dependence

A local variant of the \(CD\) test, called \(CD(p)\) test (Pesaran 2004), takes into account an appropriate subset of neighbouring cross-sectional units to check the null of no XSD against the alternative of local XSD, i.e. dependence between neighbours only. To do so, the pairs of neighbouring units are selected by means of a binary proximity matrix like those used in spatial models. In the original paper, a regular ordering of observations is assumed, so that the \(m\)-th cross-sectional observation is a neighbour to the \((m-1)\)-th and to the \((m+1)\)-th. Extending the \(CD(p)\) test to irregular lattices, we employ the binary proximity matrix as a selector for discarding the correlation coefficients relative to pairs of observations that are not neighbours in computing the \(CD\) statistic. The test is then defined as

\[CD=\sqrt{\frac{1}{\sum_{i=1}^{n-1} \sum_{j=i+1}^{n} w(p)_{ij}}} ( \sum_{i=1}^{n-1} \sum_{j=i+1}^{n} [w(p)]_{ij} \sqrt{T_{ij}}\hat{\rho}_{ij} )\]

where \([w(p)]_{ij}\) is the \((i,j)\)-th element of the \(p\)-th order proximity matrix, so that if \(h,k\) are not neighbours, \([w(p)]_{hk}=0\) and \(\hat{\rho}_{hk}\) gets “killed”; this is easily seen to reduce to formula (14) in Pesaran (Pesaran 2004) for the special case considered in that paper. The same can be applied to the \(LM\), \(SCLM\), and \(BCSCLM\) tests.

Therefore, the local version of either test can be computed supplying an \(n \times n\) matrix (of any kind coercible to logical), providing information on whether any pair of observations are neighbours or not, to the w argument. If w is supplied, only neighbouring pairs will be used in computing the test; else, w will default to NULL and all observations will be used. The matrix needs not really be binary, so commonly used “row-standardized” matrices can be employed as well: it is enough that neighbouring pairs correspond to nonzero elements in w.²⁰

Overview of functions for panel unit root testing

Below, first an overview is provided which tests are implemented per functions. A theoretical treatment is given for a few of those tests later on. The package plm offers several panel unit root tests contained in three functions:

• purtest (Levin-Lin-Chu test, IPS test, several Fisher-type tests, Hadri’s test),
• cipstest (cross-sectionally augmented IPS test), and
• phansi (Simes’ test).

While purtest implements various tests which can be selected via its test argument, cipstest and phansi are functions for a specific test each.

Function purtest offers the following tests by setting argument test to:

• "levinlin" (default), for the Levin-Lin-Chu test (Levin, Lin, and Chu (2002)), see below for a theoretical exposition (Levin-Lin-Chu test)),
• "ips", for Im-Pesaran-Shin (IPS) test by Im, Pesaran, and Shin (2003), see below for a theoretical exposition (Im-Pesaran-Shin test)),
• "madwu", is the inverse \(\chi^2\) test by Maddala and Wu (1999), also called P test by Choi (2001),
• "Pm", is the modified P test proposed by Choi (2001) for large N,
• "invnormal", is the inverse normal test (Choi (2001)),
• "logit", is the logit test (Choi (2001)),
• "hadri", for Hadri’s test (Hadri (2000)).

The tests in purtest are often called first generation panel unit root tests as they do assume absence of cross-sectional correlation; all these, except Hadri’s test, are based on the estimation of augmented Dickey-Fuller (ADF) regressions for each time series. A statistic is then computed using the t-statistics associated with the lagged variable. I a different manner, the Hadri residual-based LM statistic is the cross-sectional average of individual KPSS statistics (Kwiatkowski et al. (1992)), standardized by their asymptotic mean and standard deviation. Among the tests in purtest, "madwu", "Pm", "invormal", and "logit" are Fisher-type tests.²¹

purtest returns an object of class "purtest" which contains details about the test performed, among them details about the individual regressions/statistics for the test. Associated summary and print.summary methods can be used to extract/display the additional information.

Function cipstest implements Pesaran’s (Pesaran (2007)) cross-sectionally augmented version of the Im-Pesaran-Shin panel unit root test and is a so-called second-generation panel unit root test.

Function phansi implements the idea of Hanck (2013) to apply Simes’ testing approach for intersection of individual hypothesis tests to panel unit root testing, see below for a more thorough treatment of Simes’ approach for intersecting hypotheses.

Preliminary results

We consider the following model:

\[ y_{it} = \delta y_{it-1} + \sum_{L=1}^{p_i} \theta_i \Delta y_{it-L}+\alpha_{mi} d_{mt}+\epsilon_{it} \]

The unit root hypothesis is \(\rho = 1\). The model can be rewritten in difference:

\[ \Delta y_{it} = \rho y_{it-1} + \sum_{L=1}^{p_i} \theta_i \Delta y_{it-L}+\alpha_{mi} d_{mt}+\epsilon_{it} \]

So that the unit-root hypothesis is now \(\rho = 0\).

Some of the unit-root tests for panel data are based on preliminary results obtained by running the above Augmented Dickey-Fuller (ADF) regression.

First, we have to determine the optimal number of lags \(p_i\) for each time-series. Several possibilities are available. They all have in common that the maximum number of lags have to be chosen first. Then, \(p_i\) can be chosen by using:

• the Schwarz information criterion (SIC) (also known as Bayesian information criterion (BIC)),
• the Akaike information criterion (AIC),
• the Hall’s method, which consist in removing the higher lags while they are not significant.

The ADF regression is run on \(T-p_i-1\) observations for each individual, so that the total number of observations is \(n\times \tilde{T}\) where \(\tilde{T}=T-p_i-1\)

\(\bar{p}\) is the average number of lags. Call \(e_{i}\) the vector of residuals.

Estimate the variance of the \(\epsilon_i\) as:

\[ \hat{\sigma}_{\epsilon_i}^2 = \frac{\sum_{t=p_i+1}^{T} e_{it}^2}{df_i} \]

Levin-Lin-Chu model

Then, as per Levin, Lin, and Chu (2002), compute artificial regressions of \(\Delta y_{it}\) and \(y_{it-1}\) on \(\Delta y_{it-L}\) and \(d_{mt}\) and get the two vectors of residuals \(z_{it}\) and \(v_{it}\).

Standardize these two residuals and run the pooled regression of \(z_{it}/\hat{\sigma}_i\) on \(v_{it}/\hat{\sigma}_i\) to get \(\hat{\rho}\), its standard deviation \(\hat{\sigma}({\hat{\rho}})\) and the t-statistic \(t_{\hat{\rho}}=\hat{\rho}/\hat{\sigma}({\hat{\rho}})\).

Compute the long run variance of \(y_i\) :

\[ \hat{\sigma}_{yi}^2 = \frac{1}{T-1}\sum_{t=2}^T \Delta y_{it}^2 + 2 \sum_{L=1}^{\bar{K}}w_{\bar{K}L}\left[\frac{1}{T-1}\sum_{t=2+L}^T \Delta y_{it} \Delta y_{it-L}\right] \]

Define \(\bar{s}_i\) as the ratio of the long and short term variance and \(\bar{s}\) the mean for all the individuals of the sample

\[ s_i = \frac{\hat{\sigma}_{yi}}{\hat{\sigma}_{\epsilon_i}} \]

\[ \bar{s} = \frac{\sum_{i=1}^n s_i}{n} \]

\[ t^*_{\rho}=\frac{t_{\rho}- n \bar{T} \bar{s} \hat{\sigma}_{\tilde{\epsilon}}^{-2} \hat{\sigma}({\hat{\rho}}) \mu^*_{m\tilde{T}}}{\sigma^*_{m\tilde{T}}} \]

follows a normal distribution under the null hypothesis of stationarity. \(\mu^*_{m\tilde{T}}\) and \(\sigma^*_{m\tilde{T}}\) are given in table 2 of the original paper and are also available in the package.

An example how the Levin-Lin-Chu test is performed with purtest using a lag of 2 and intercept and a time trend as exogenous variables in the ADF regressions is:

data("HousePricesUS", package = "pder")
lprice <- log(pdata.frame(HousePricesUS)$price)
(lev <- purtest(lprice, test = "levinlin", lags = 2, exo = "trend"))

## 
##  Levin-Lin-Chu Unit-Root Test (ex. var.: Individual Intercepts and
##  Trend)
## 
## data:  lprice
## z = -1.2573, p-value = 0.1043
## alternative hypothesis: stationarity

summary(lev) ### gives details

## Levin-Lin-Chu Unit-Root Test 
## Exogenous variables: Individual Intercepts and Trend 
## User-provided lags
## statistic: -1.257 
## p-value: 0.104 
## 
##    lags obs          rho        trho      p.trho     sigma2ST     sigma2LT
## 1     2  26 -0.092065357 -1.66309731 0.767613204 0.0003143120 0.0004013788
## 4     2  26 -0.124093984 -1.29563385 0.888755668 0.0010950144 0.0014736172
## 5     2  26 -0.104647566 -1.10814627 0.926357866 0.0007296044 0.0007451534
## 6     2  26 -0.219022744 -2.94312106 0.148774635 0.0007716609 0.0048254402
## 8     2  26 -0.052471794 -0.95375744 0.948405601 0.0006375257 0.0028152736
## 9     2  26 -0.181914333 -2.73331072 0.222919642 0.0021489671 0.0064455696
## 10    2  26 -0.232215125 -3.37321191 0.054989191 0.0005566400 0.0024147067
## 11    2  26 -0.356452679 -4.35943612 0.002479709 0.0008542529 0.0045574510
## 12    2  26  0.279936991  1.83482002 0.999998365 0.0004172617 0.0012951914
## 13    2  26 -0.062610441 -0.84216587 0.960499065 0.0003168316 0.0002981994
## 16    2  26 -0.159254884 -2.29683734 0.435109226 0.0007437190 0.0010203969
## 17    2  26 -0.237065476 -4.05050006 0.007367490 0.0005512405 0.0009463645
## 18    2  26 -0.140788644 -2.08598977 0.553093684 0.0005079423 0.0005697978
## 19    2  26 -0.099218199 -1.83853581 0.686000079 0.0008756343 0.0020399374
## 20    2  26 -0.046049208 -0.85174237 0.959567857 0.0003914722 0.0011128571
## 21    2  26 -0.102633777 -1.81503721 0.697696805 0.0004063481 0.0002345858
## 22    2  26 -0.115700485 -1.72146553 0.741996511 0.0010094113 0.0047169602
## 23    2  26 -0.218251170 -2.90863990 0.159598845 0.0010530905 0.0031254005
## 24    2  26 -0.293126134 -3.65827755 0.025157448 0.0004297907 0.0012860757
## 25    2  26 -0.107476475 -2.33427946 0.414643019 0.0008718077 0.0083545213
## 26    2  26 -0.135655633 -2.06664416 0.563904448 0.0009443422 0.0007705242
## 27    2  26 -0.005168776 -0.06565125 0.995432637 0.0007059964 0.0017093982
## 28    2  26 -0.101736562 -1.02991147 0.938382427 0.0008552994 0.0003212357
## 29    2  26 -0.106917037 -1.44344289 0.848353136 0.0004659842 0.0004528307
## 30    2  26 -0.143955051 -1.60594256 0.791068910 0.0016589513 0.0022706981
## 31    2  26 -0.093688191 -1.92279670 0.642427798 0.0004025885 0.0009260538
## 32    2  26 -0.313691108 -2.30732500 0.429359754 0.0009640889 0.0019245305
## 33    2  26 -0.151599029 -2.54586869 0.305755540 0.0011446680 0.0072973098
## 34    2  26 -0.113830637 -2.06152082 0.566761534 0.0008514932 0.0055913854
## 35    2  26 -0.220363663 -1.72391205 0.740887358 0.0005138068 0.0007015982
## 36    2  26 -0.211779244 -3.98522621 0.009148144 0.0011004960 0.0062947120
## 37    2  26 -0.161851244 -2.19906397 0.489557509 0.0002334460 0.0001298656
## 38    2  26 -0.222507555 -1.45762738 0.843900682 0.0048242337 0.0019584921
## 39    2  26 -0.119405321 -2.41405422 0.372087172 0.0004737829 0.0009459741
## 40    2  26 -0.066956522 -0.94176615 0.949844273 0.0011477969 0.0044950987
## 41    2  26 -0.107420235 -2.09998836 0.545259606 0.0009669881 0.0033414294
## 42    2  26 -0.211640785 -3.51839705 0.037371037 0.0004302930 0.0020971822
## 44    2  26 -0.160491538 -2.32116399 0.421790841 0.0016866384 0.0067053791
## 45    2  26  0.013957358  0.21048073 0.998138553 0.0002474183 0.0001310960
## 46    2  26 -0.125206819 -1.36523187 0.871040437 0.0014782610 0.0005893232
## 47    2  26 -0.146576570 -2.35613125 0.402831448 0.0002851628 0.0001796349
## 48    2  26 -0.106184312 -1.41243370 0.857717409 0.0006722417 0.0029218865
## 49    2  26 -0.110328029 -2.31986075 0.422503260 0.0007413810 0.0032547907
## 50    2  26 -0.336849990 -3.07534065 0.112195130 0.0015064070 0.0017150678
## 51    2  26 -0.219041498 -2.26882562 0.450564573 0.0004175437 0.0010455238
## 53    2  26 -0.249921002 -2.67545341 0.246874228 0.0008780514 0.0016600448
## 54    2  26 -0.092856496 -1.36610183 0.870804641 0.0011141161 0.0011656778
## 55    2  26 -0.119379994 -2.19438511 0.492186948 0.0008204933 0.0007526503
## 56    2  26 -0.094196887 -1.53868461 0.816472629 0.0016308380 0.0062631397

Im-Pesaran-Shin (IPS) test

This test by Im, Pesaran, and Shin (2003) does not require that \(\rho\) is the same for all the individuals. The null hypothesis is still that all the series have an unit root, but the alternative is that some may have a unit root and others have different values of \(\rho_i <0\).

The test is based on the average of the student statistic of the \(\rho\) obtained for each individual:

\[ \bar{t}=\frac{1}{n}\sum_{i=1}^n t_{\rho i} \]

The statistic is then:

\[ z = \frac{\sqrt{n}\left(\bar{t}- E(\bar{t})\right)}{\sqrt{V(\bar{t})}} \]

\(\mu^*_{m\tilde{T}}\) and \(\sigma^*_{m\tilde{T}}\) are given in table 2 of the original paper and are also available in the package.

An example of the IPS test with purtest with the same settings as in the previously performed Levin-Lin-Chu test is:

purtest(lprice, test = "ips", lags = 2, exo = "trend")

## 
##  Im-Pesaran-Shin Unit-Root Test (ex. var.: Individual Intercepts and
##  Trend)
## 
## data:  lprice
## Wtbar = 0.76622, p-value = 0.7782
## alternative hypothesis: stationarity

Simes’ approach: intersecting hypotheses

A different approach to panel unit root testing can be drawn from the general Simes’ test for intersection of individual hypothesis tests (Simes 1986). Hanck (2013) suggests to apply the approach for panel unit root testing: The tests works by combining p-values from single hypothesis tests (individual unit root tests) with a global (intersected) hypothesis and controls for the multiplicity in testing. Thus, it works “on top” of any panel unit root test which yield a p-value for each individual series. Unlike most other panel unit root tests, this approach allows to discriminate between individuals for which the individual H0 (unit root present for individual series) is rejected/is not rejected and requires a pre-specified significance level. Further, the test is robust versus general patterns of cross-sectional dependence.

The function phansi for this test takes as main input object either a numeric containing p-values of individual tests or a "purtest" object as produced by function purtest which holds a suitable pre-computed panel unit root test (one that produces p-values per individual series). The significance level is set by argument alpha (default 5 %). The function’s return value is a list with detailed evaluation of the applied Simes test. The associated print method gives a verbal evaluation.

The following examples shows both accepted ways of input, the first example replicates Hanck (2013), table 11 (left side), who applied some panel unit root test for a Purchasing Power Parity analysis per country (individual H0 hypotheses per series) to get the individual p-values and then used Simes’ approach for testing the global (intersecting) hypothesis for the whole panel.

### input is numeric (p-values), replicates Hanck (2013), Table 11 (left side)
pvals <- c(0.0001,0.0001,0.0001,0.0001,0.0001,0.0001,0.0050,0.0050,0.0050,
           0.0050,0.0175,0.0175,0.0200,0.0250,0.0400,0.0500,0.0575,0.2375,0.2475)
countries <- c("Argentina","Sweden","Norway","Mexico","Italy","Finland","France",
              "Germany","Belgium","U.K.","Brazil","Australia","Netherlands",
              "Portugal","Canada", "Spain","Denmark","Switzerland","Japan")
names(pvals) <- countries
h <- phansi(pvals)
print(h)

## 
##         Simes Test as Panel Unit Root Test (Hanck (2013))
## 
## H0: All individual series have a unit root
## HA: Stationarity for at least some individuals
## 
## Alpha: 0.05
## Number of individuals: 19
## 
## Evaluation:
##  H0 rejected (globally)
## 
##  Individual H0 rejected for 10 individual(s) (integer id(s)):
##   1, 2, 3, 4, 5, 6, 7, 8, 9, 10

h$rejected # logical indicating the individuals with rejected individual H0

##   Argentina      Sweden      Norway      Mexico       Italy     Finland 
##        TRUE        TRUE        TRUE        TRUE        TRUE        TRUE 
##      France     Germany     Belgium        U.K.      Brazil   Australia 
##        TRUE        TRUE        TRUE        TRUE       FALSE       FALSE 
## Netherlands    Portugal      Canada       Spain     Denmark Switzerland 
##       FALSE       FALSE       FALSE       FALSE       FALSE       FALSE 
##       Japan 
##       FALSE

### input is a (suitable) purtest object / different example
y <- data.frame(split(Grunfeld$inv, Grunfeld$firm))
obj <- purtest(y, pmax = 4, exo = "intercept", test = "madwu")
phansi(obj, alpha = 0.06) # test with significance level set to 6 %

The Laird-Ware representation for mixed models

A general representation for the linear mixed effects model is given in Laird and Ware (1982).

\[ \begin{array}{rcl} y_{it} & = & \beta_1 x_{1ij} + \dots + \beta_p x_{pij} \\ & & b_1 z_{1ij} + \dots + b_p z_{pij} + \epsilon_{ij} \\ b_{ik} & \sim & N(0,\psi^2_k), \phantom{p} Cov(b_k,b_{k'}) = \psi_{kk'} \\ \epsilon_{ij} & \sim & N(0,\sigma^2 \lambda_{ijj}), \phantom{p} Cov(\epsilon_{ij},\epsilon_{ij'}) = \sigma^2 \lambda_{ijj'} \\ \end{array} \]

where the \(x_1, \dots x_p\) are the fixed effects regressors and the \(z_1, \dots z_p\) are the random effects regressors, assumed to be normally distributed across groups. The covariance of the random effects coefficients \(\psi_{kk'}\) is assumed constant across groups and the covariances between the errors in group \(i\), \(\sigma^2 \lambda_{ijj'}\), are described by the term \(\lambda_{ijj'}\) representing the correlation structure of the errors within each group (e.g., serial correlation over time) scaled by the common error variance \(\sigma^2\).

Pooling and Within

The pooling specification in plm is equivalent to a classical linear model (i.e., no random effects regressor and spherical errors: \(b_{iq}=0 \phantom{p} \forall i,q, \phantom{p} \lambda_{ijj}=\sigma^2\) for \(j=j'\), \(0\) else). The within one is the same with the regressors’ set augmented by \(n-1\) group dummies. There is no point in using nlme as parameters can be estimated by OLS which is also ML.

Random effects

In the Laird and Ware notation, the RE specification is a model with only one random effects regressor: the intercept. Formally, \(z_{1ij}=1 \phantom{p}\forall i,j, \phantom{p} z_{qij}=0 \phantom{p} \forall i, \forall j, \forall q \neq 1\) \(\lambda_{ij}=1\) for \(i=j\), \(0\) else). The composite error is therefore \(u_{ij}=1b_{i1} + \epsilon_{ij}\). Below we report coefficients of Grunfeld’s model estimated by GLS and then by ML:

library(nlme)
reGLS <- plm(inv~value+capital, data=Grunfeld, model="random")

reML <- lme(inv~value+capital, data=Grunfeld, random=~1|firm)

coef(reGLS)

## (Intercept)       value     capital 
## -57.8344149   0.1097812   0.3081130

summary(reML)$coefficients$fixed

## (Intercept)       value     capital 
## -57.8644245   0.1097897   0.3081881

Variable coefficients, “random”

Swamy’s variable coefficients model (Swamy 1970) has coefficients varying randomly (and independently of each other) around a set of fixed values, so the equivalent specification is \(z_{q}=x_{q} \phantom{p} \forall q\), i.e. the fixed effects and the random effects regressors are the same, and \(\psi_{kk'}=\sigma_\mu^2 I_N\), and \(\lambda_{ijj}=1\), \(\lambda_{ijj'}=0\) for \(j \neq j'\), that’s to say they are not correlated.

Estimation of a mixed model with random coefficients on all regressors is rather demanding from the computational side. Some models from our examples fail to converge. The below example is estimated on the Grunfeld data and model with time effects.

vcm <- pvcm(inv~value+capital, data=Grunfeld, model="random", effect="time")

vcmML <- lme(inv~value+capital, data=Grunfeld, random=~value+capital|year)

coef(vcm)

## (Intercept)       value     capital 
## -18.5538638   0.1239595   0.1114579

summary(vcmML)$coefficients$fixed

## (Intercept)       value     capital 
## -26.3558395   0.1241982   0.1381782

Variable coefficients, “within”

This specification actually entails separate estimation of \(T\) different standard linear models, one for each group in the data, so the estimation approach is the same: OLS. In nlme this is done by creating an lmList object, so that the two models below are equivalent (output suppressed):

vcmf <- pvcm(inv~value+capital, data=Grunfeld, model="within", effect="time")

vcmfML <- lmList(inv~value+capital|year, data=Grunfeld)

General FGLS

The general, or unrestricted, feasible GLS (FGLS), pggls in the plm nomenclature, is equivalent to a model with no random effects regressors (\(b_{iq}=0 \phantom{p} \forall i,q\)) and an error covariance structure which is unrestricted within groups apart from the usual requirements. The function for estimating such models with correlation in the errors but no random effects is gls().

This very general serial correlation and heteroskedasticity structure is not estimable for the original Grunfeld data, which have more time periods than firms, therefore we restrict them to firms 4 to 6.

sGrunfeld <- Grunfeld[Grunfeld$firm %in% 4:6, ]

ggls <- pggls(inv~value+capital, data=sGrunfeld, model="pooling")

gglsML <- gls(inv~value+capital, data=sGrunfeld,
              correlation=corSymm(form=~1|year))

coef(ggls)

## (Intercept)       value     capital 
##  1.19679342  0.10555908  0.06600166

summary(gglsML)$coefficients

## (Intercept)       value     capital 
##  -2.4156266   0.1163550   0.0735837

The within case is analogous, with the regressor set augmented by \(n-1\) group dummies.

AR(1) pooling or random effects panel

Linear models with groupwise structures of time-dependence³¹ may be fitted by gls(), specifying the correlation structure in the correlation option³²:

Grunfeld$year <- as.numeric(as.character(Grunfeld$year))
lmAR1ML <- gls(inv~value+capital,data=Grunfeld,
               correlation=corAR1(0,form=~year|firm))

and analogously the random effects panel with, e.g., AR(1) errors (see Baltagi (2005); Baltagi (2013); Baltagi (2021), ch. 5), which is a very common specification in econometrics, may be fit by lme specifying an additional random intercept:

reAR1ML <- lme(inv~value+capital, data=Grunfeld,random=~1|firm,
               correlation=corAR1(0,form=~year|firm))

The regressors’ coefficients and the error’s serial correlation coefficient may be retrieved this way:

summary(reAR1ML)$coefficients$fixed

##  (Intercept)        value      capital 
## -40.27650822   0.09336672   0.31323330

coef(reAR1ML$modelStruct$corStruct, unconstrained=FALSE)

##      Phi 
## 0.823845

Significance statistics for the regressors’ coefficients are to be found in the usual summary object, while to get the significance test of the serial correlation coefficient one can do a likelihood ratio test as shown in the following.

An LR test for serial correlation and one for random effects

A likelihood ratio test for serial correlation in the idiosyncratic residuals can be done as a nested models test, by anova(), comparing the model with spherical idiosyncratic residuals with the more general alternative featuring AR(1) residuals. The test takes the form of a zero restriction test on the autoregressive parameter.

This can be done on pooled or random effects models alike. First we report the simpler case.

We already estimated the pooling AR(1) model above. The GLS model without correlation in the residuals is the same as OLS, and one could well use lm() for the restricted model. Here we estimate it by gls().

lmML <- gls(inv~value+capital, data=Grunfeld)
anova(lmML, lmAR1ML)

##         Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## lmML        1  4 2400.217 2413.350 -1196.109                        
## lmAR1ML     2  5 2094.936 2111.352 -1042.468 1 vs 2 307.2813  <.0001

The AR(1) test on the random effects model is to be done in much the same way, using the random effects model objects estimated above:

anova(reML, reAR1ML)

##         Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## reML        1  5 2205.851 2222.267 -1097.926                        
## reAR1ML     2  6 2094.802 2114.501 -1041.401 1 vs 2 113.0496  <.0001

A likelihood ratio test for random effects compares the specifications with and without random effects and spherical idiosyncratic errors:

anova(lmML, reML)

##      Model df      AIC      BIC    logLik   Test L.Ratio p-value
## lmML     1  4 2400.217 2413.350 -1196.109                       
## reML     2  5 2205.851 2222.267 -1097.926 1 vs 2 196.366  <.0001

The random effects, AR(1) errors model in turn nests the AR(1) pooling model, therefore a likelihood ratio test for random effects sub AR(1) errors may be carried out, again, by comparing the two autoregressive specifications:

anova(lmAR1ML, reAR1ML)

##         Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## lmAR1ML     1  5 2094.936 2111.352 -1042.468                        
## reAR1ML     2  6 2094.802 2114.501 -1041.401 1 vs 2 2.134349   0.144

whence we see that the Grunfeld model specification doesn’t seem to need any random effects once we control for serial correlation in the data.