CS计算机代考程序代写 Hive flex GMM scheme data structure JSS

JSS
Journal of Statistical Software
July 2008, Volume 27, Issue 2. http://www.jstatsoft.org/
Panel Data Econometrics in R: The plm Package Yves Croissant Giovanni Millo
Universit ́e Lumi`ere Lyon 2 University of Trieste and Generali SpA
Abstract
Panel data econometrics is obviously one of the main fields in the profession, but most of the models used are difficult to estimate with R. plm is a package for R which intends to make the estimation of linear panel models straightforward. plm provides functions to estimate a wide variety of models and to make (robust) inference.
Keywords: panel data, covariance matrix estimators, generalized method of moments, R.
1. Introduction
Panel data econometrics is a continuously developing field. The increasing availability of data observed on cross-sections of units (like households, firms, countries etc.) and over time has given rise to a number of estimation approaches exploiting this double dimensionality to cope with some of the typical problems associated with economic data, first of all that of unobserved heterogeneity.
Timewise observation of data from different observational units has long been common in other fields of statistics (where they are often termed longitudinal data). In the panel data field as well as in others, the econometric approach is nevertheless peculiar with respect to experimental contexts, as it is emphasizing model specification and testing and tackling a number of issues arising from the particular statistical problems associated with economic data.
Thus, while a very comprehensive software framework for (among many other features) max- imum likelihood estimation of linear regression models for longitudinal data, packages nlme (Pinheiro, Bates, DebRoy, and Sarkar 2007) and lme4 (Bates 2007), is available in the R (R Development Core Team 2008) environment and can be used, e.g., for estimation of random effects panel models, its use is not intuitive for a practicing econometrician, and maximum likelihood estimation is only one of the possible approaches to panel data econometrics. More- over, economic panel datasets often happen to be unbalanced (i.e., they have a different number

2 Panel Data Econometrics in R: The plm Package
of observations between groups), which case needs some adaptation to the methods and is not compatible with those in nlme. Hence the need for a package doing panel data “from the econometrician’s viewpoint” and featuring at a minimum the basic techniques econometri- cians themselves are used to: random and fixed effects estimation of static linear panel data models, variable coefficients models, generalized method of moments estimation of dynamic models; and the basic toolbox of specification and misspecification diagnostics.
Furthermore, we felt there was the need for automation of some basic data management tasks as lagging, summing and, more in general, applying (in the R sense) functions to the data, which, although conceptually simple, become cumbersome and error-prone on two- dimensional data, especially in the case of unbalanced panels.
The result of our work is bundled in the plm add-on package, available from the Comprehensive R Archive Network at http://CRAN.R-project.org/package=plm.
The paper is organized as follows: Section 2 presents a very short overview of the typical model taxonomy1. Section 3 discusses the software approach used in the package. The next three sections present the functionalities of the package in more detail: data management (Section 4), estimation (Section 5) and testing (Section 6), giving a short description and illustrating them with examples. Section 7 compares the approach in plm to that of nlme and lme4, highlighting the features of the latter two that an econometrician might find most useful. Section 8 concludes the paper.
2. The linear panel model
The basic linear panel models used in econometrics can be described through suitable restric- tions of the following general model:
y =α +β⊤x +u (1) it it it it it
where i = 1,…n is the individual (group, country, …) index, t = 1,…,T is the time index and uit a random disturbance term of mean 0.
Of course the latter is not estimable with N = n × T data points. A number of assumptions are usually made about the parameters, the errors and the exogeneity of the regressors, giving rise to a taxonomy of feasible models for panel data.
The most common one is parameter homogeneity, which means that αit = α for all i, t and βit = β for all i, t. The resulting model
yit = α + β⊤xit + uit (2) is a standard linear model pooling all the data across i and t.
To model individual heterogeneity, one often assumes that the error term has two separate components, one of which is specific to the individual and does not change over time2. This is called the unobserved effects model:
1Comprehensive treatments are to be found in many econometrics textbooks, e.g. Baltagi (2001) or Wooldridge (2002): the reader is referred to these, especially to the first nine chapters of Baltagi (2001).
2For the sake of exposition we are considering only the individual effects case here. There may also be time effects, which is a symmetric case, or both of them, so that the error has three components: uit = μi +λt +εit.

Journal of Statistical Software 3
yit =α+β⊤xit +μi +εit (3)
The appropriate estimation method for this model depends on the properties of the two error components. The idiosyncratic error εit is usually assumed well-behaved and independent from both the regressors xit and the individual error component μi. The individual component may be in turn either independent from the regressors or correlated.
If it is correlated, the ordinary least squares (OLS) estimator for β would be inconsistent, so it is customary to treat the μi as a further set of n parameters to be estimated, as if in the general model αit = αi for all t. This is called the fixed effects (also known as within or least squares dummy variables) model, usually estimated by OLS on transformed data, and gives consistent estimates for β.
If the individual-specific component μi is uncorrelated with the regressors, a situation which is usually termed random effects, the overall error uit also is, so the OLS estimator is consistent. Nevertheless, the common error component over individuals induces correlation across the composite error terms, making OLS estimation inefficient, so one has to resort to some form of feasible generalized least squares (GLS) estimators. This is based on the estimation of the variance of the two error components, for which there are a number of different procedures available.
If the individual component is missing altogether, pooled OLS is the most efficient estimator for β. This set of assumptions is usually labelled pooling model, although this actually refers to the errors’ properties and the appropriate estimation method rather than the model itself. If one relaxes the usual hypotheses of well-behaved, white noise errors and allows for the id- iosyncratic error εit to be arbitrarily heteroskedastic and serially correlated over time, a more general kind of feasible GLS is needed, called the unrestricted or general GLS. This specifi- cation can also be augmented with individual-specific error components possibly correlated with the regressors, in which case it is termed fixed effects GLS.
Another way of estimating unobserved effects models through removing time-invariant indi- vidual components is by first-differencing the data: lagging the model and subtracting, the time-invariant components (the intercept and the individual error component) are eliminated, and the model
∆yit = β⊤∆xit + ∆uit (4)
(where ∆yit = yit − yi,t−1, ∆xit = xit − xi,t−1 and, from (3), ∆uit = uit − ui,t−1 = ∆εit for t = 2, …, T ) can be consistently estimated by pooled OLS. This is called the first-difference, or FD estimator. Its relative efficiency, and so reasons for choosing it against other consis- tent alternatives, depends on the properties of the error term. The FD estimator is usually preferred if the errors uit are strongly persistent in time, because then the ∆uit will tend to be serially uncorrelated.
Lastly, the between model, which is computed on time (group) averages of the data, discards all the information due to intragroup variability but is consistent in some settings (e.g., non- stationarity) where the others are not, and is often preferred to estimate long-run relationships.
Variable coefficients models relax the assumption that βit = β for all i, t. Fixed coefficients models allow the coefficients to vary along one dimension, like βit = βi for all t. Random

4 Panel Data Econometrics in R: The plm Package
coefficients models instead assume that coefficients vary randomly around a common average,
as βit = β + ηi for all t, where ηi is a group- (time-) specific effect with mean zero.
The hypotheses on parameters and error terms (and hence the choice of the most appropriate
estimator) are usually tested by means of:
􏰀 pooling tests to check poolability, i.e., the hypothesis that the same coefficients apply across all individuals,
􏰀 if the homogeneity assumption over the coefficients is established, the next step is to establish the presence of unobserved effects, comparing the null of spherical residuals with the alternative of group (time) specific effects in the error term,
􏰀 the choice between fixed and random effects specifications is based on Hausman-type tests, comparing the two estimators under the null of no significant difference: if this is not rejected, the more efficient random effects estimator is chosen,
􏰀 even after this step, departures of the error structure from sphericity can further affect inference, so that either screening tests or robust diagnostics are needed.
Dynamic models and in general lack of strict exogeneity of the regressors, pose further prob- lems to estimation which are usually dealt with in the generalized method of moments (GMM) framework.
These were, in our opinion, the basic requirements of a panel data econometrics package for the R language and environment. Some, as often happens with R, were already fulfilled by packages developed for other branches of computational statistics, while others (like the fixed effects or the between estimators) were straightforward to compute after transforming the data, but in every case there were either language inconsistencies w.r.t. the standard econometric toolbox or subtleties to be dealt with (like, for example, appropriate computation of standard errors for the demeaned model, a common pitfall), so we felt there was need for an “all in one” econometrics-oriented package allowing to make specification searches, estimation and inference in a natural way.
3. Software approach
3.1. Data structure
Panel data have a special structure: each row of the data corresponds to a specific individual and time period. In plm the data argument is an ordinary data.frame, but an argument called index has to be added to indicate the structure of the data. This can be:
􏰀 NULL (the default value), it is then assumed that the first two columns contain the individual and the time index and that observations are ordered by individual and by time period,
􏰀 a character string, which should be the name of the individual index,
􏰀 a character vector of length two containing the names of the individual and the time index,

Journal of Statistical Software 5 􏰀 an integer which is the number of individuals (only in case of a balanced panel with
observations ordered by individual).
The plm.data function is then called, which returns a data.frame with the two first columns containing the individual and the time indexes.
3.2. Interface
Estimation interface
plm provides four functions for estimation:
􏰀 plm: estimation of the basic panel models, i.e., within, between and random effect
models. Models are estimated using the lm function to 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.
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 both3,
􏰀 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
3Although in most models the individual and time effects cases are symmetric, there are exceptions: es- timating the fd model on time effects is meaningless because cross-sections do not generally have a natural ordering, so here the effect will always be set to “individual”.

6 Panel Data Econometrics in R: The plm Package
inexpensive. Therefore most examples in the following are based on formula methods, which
are perhaps the cleanest for illustrative purposes.
3.3. Computational approach to estimation
The feasible GLS methods needed for efficient estimation of unobserved effects models have a simple closed-form solution: once the variance components have been estimated and hence the covariance matrix of errors Vˆ , model parameters can be estimated as
βˆ = (X⊤Vˆ −1X)−1(X⊤Vˆ −1y) (5)
Nevertheless, in practice plain computation of βˆ has long been an intractable problem even for moderate-sized datasets because of the need to invert the N × N Vˆ matrix. With the advances in computer power, this is no more so, and it is possible to program the “naive” estimator (5) in R with standard matrix algebra operators and have it working seamlessly for the standard “guinea pigs”, e.g., the Grunfeld data. Estimation with a couple of thousands of data points also becomes feasible on a modern machine, although excruciatingly slow and definitely not suitable for everyday econometric practice. Memory limits would also be very near because of the storage needs related to the huge Vˆ matrix. An established solution exists for the random effects model which reduces the problem to an ordinary least squares computation.
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 βs is concerned) to OLS on partially demeaned data, where partial demeaning is defined as:
yit − θy ̄i = (Xit − θX ̄i)β + (uit − θu ̄i) (6)
where θ = 1−[σu2/(σu2 +Tσe2)]1/2, y ̄ and X ̄ denote time means of y and X, and the disturbance vit − θv ̄i is homoskedastic and serially uncorrelated. Thus the feasible RE estimate for β may be obtained estimating θˆ and running an OLS regression on the transformed data with lm(). The other estimators can be computed as special cases: for θ = 1 one gets the fixed effects estimator, for θ = 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

Journal of Statistical Software 7 standard OLS formulas to the demeaned data, which are contained inside plm objects. This
will be the subject of Section 3.4.
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 (5). 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 2007) and SparseM (see Koenker and Ng 2007). Some optimized linear algebra routines are available in the R package kinship (see Atkinson and Therneau 2007) which exploit the particular block- diagonal and symmetric structure of Vˆ making it possible to implement a fast and reliable full-matrix solution to problems of any practically relevant size.
The 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 (5) to compute βˆ, resulting in speed comparable to that of the demeaning solution.
3.4. Inference in the panel model
General frameworks for restrictions and linear hypotheses testing are available in the R en- vironment4. These are based on the Wald test, constructed as βˆ⊤Vˆ−1βˆ, where βˆ and Vˆ are consistent estimates of β and V (β), The Wald test may be used for zero-restriction (i.e., signifi- cance) testing and, more generally, for linear hypotheses in the form (Rβˆ−r)⊤[RVˆ R⊤]−1(Rβˆ− r)5. To be applicable, the test functions require extractor methods for coefficients’ and covari- ance matrix estimates to be defined for the model object to be tested. Model objects in plm all have coef() and vcov() methods and are therefore compatible with the above functions.
In the same framework, robust inference is accomplished substituting (“plugging in”) a robust estimate of the coefficient covariance matrix into the Wald statistic formula. In the panel context, the estimator of choice is the White system estimator. This called for a flexible method for computing robust coefficient covariance matrices `a la White for plm objects.
A general White system estimator for panel data is:
n
VˆR(β) = (X⊤X)−1 􏰆 Xi⊤EiXi(X⊤X)−1 (7)
i=1
where Ei is a function of the residuals eˆit, t = 1, . . . T chosen according to the relevant
heteroskedasticity and correlation structure. Moreover, it turns out that the White covariance
4See packages lmtest (Zeileis and Hothorn 2002) and car (Fox 2007).
5Moreover, coeftest() provides a compact way of looking at coefficient estimates and significance diag- nostics.

8 Panel Data Econometrics in R: The plm Package
matrix calculated on the demeaned model’s regressors and residuals (both part of plm objects) is a consistent estimator of the relevant model’s parameters’ covariance matrix, thus the method is readily applicable to models estimated by random or fixed effects, first difference or pooled OLS methods. Different pre-weighting schemes taken from package sandwich (Zeileis 2004) are also implemented to improve small-sample performance. Robust estimators with any combination of covariance structures and weighting schemes can be passed on to the testing functions.
4. Managing data and formulas
The package is now illustrated by application to some well-known examples. It is loaded using
R> library(“plm”)
The four datasets used are EmplUK data available in plm which was used by Arellano and Bond (1991), the Grunfeld data available in package Ecdat (Croissant 2006) which is used in several econometric books, the Produc data used by Munnell (1990) and the Wages used by Cornwell and Rupert (1988).
R> data(“EmplUK”, package = “plm”)
R> data(“Produc”, package = “Ecdat”)
R> data(“Grunfeld”, package = “Ecdat”)
R> data(“Wages”, package = “Ecdat”)
4.1. Data structure
As observed above, the current version of plm is capable of working with a regular data.frame without any further transformation, provided that the individual and time indexes are in the first two columns, as in all the example datasets but Wages. If this were not the case, an index optional argument would have to be passed on to the estimating and testing functions.
As an example, as the Wages dataset lacks indexes altogether, but we know it is a balanced panel already sorted in (id, time) order, these can be explicitly reconstructed and put into the first two columns by plm.data specifying the number of individuals:
R> Wages <- plm.data(Wages, index = 595) once for all. As an alternative, one could keep specifying index = 595 into every command. Datasets that have indexes but are, or are suspected to be, in the wrong order can be set in the same way, specifying the names of the id and time indexes instead of the number of individuals. For the sake of illustration, here is an example on the EmplUK dataset (although this is in the correct order already): R> EmplUK <- plm.data(EmplUK, index = c("firm", "year")) d=00 1−1...00 . . . . .. . . ....... 0 0 0 0 ... 1 −1 Journal of Statistical Software 9 4.2. Data transformation Panel data estimation requires to apply different transformations to raw series. If x is a series of length nT (where n is the number of individuals and T is the number of time periods), the transformed series x ̃ is obtained as x ̃ = Mx where M is a transformation matrix. Denoting j a vector of one of length T and In the identity matrix of dimension n, we get: 􏰀 the between transformation: P = 1 In ⊗ jj⊤ returns a vector containing the individual T means. The Between and between functions performs this operation, the first one returning a vector of length nT , the second one a vector of length n, 􏰀 the within transformation: Q = InT − P returns a vector containing the values in deviation from the individual means. The within function performs this operation. 􏰀 the first difference transformation D = In ⊗ d where 1−10 0...00 0 1 −1 0 ...0 0  is of dimension (T − 1, T ). Note that R’s diff() and lag() functions do not compute correctly these transformations for panel data because they are unable to identify when there is a change in individual in the data. They have been therefore rewritten in order to handle correctly panel data. The Between, between, within, lag and diff functions can be used directly in formulas to estimate models in all the estimation functions previously described. 4.3. Dynamic formula Using the lag and diff functions to write a formula can be very cumbersome for dynamic models. A dynformula function is provided to ease the writing of such formulas. Suppose for example than one wants to estimate a model where employment depends on its own first two lags and on the second and third lags of wages and capital, all the variables except capital being in logs, and capital being differenced. The formula is then: R> log(emp) ~ lag(log(emp), 1) + lag(log(emp), 2) + lag(log(wage), 2) +
+ lag(log(wage), 3) + diff(capital, 2) + diff(capital, 3)
which can be easily obtained with dynformula:
R> dynformula(emp ~ wage + capital, log = list(capital = FALSE, TRUE),
+ lag = list(emp = 2, c(2, 3)), diff = list(FALSE, capital = TRUE))
log(emp) ~ lag(log(emp), 1) + lag(log(emp), 2) + lag(log(wage),
2) + lag(log(wage), 3) + diff(capital, 2) + diff(capital,
3)

10 Panel Data Econometrics in R: The plm Package The arguments lag, diff and log are lists which can be:
􏰀 unnamed, in this case, the length of the list must equal the number of variables,
􏰀 named, in this case the missing variables get the default values (no lags, no logs, no
differences),
􏰀 partially named, in this case the unnamed element is the user defined default value.
The elements of diff and log lists should be logical values. For lag, it is either a number or a vector of two numbers: 3 means current value and first three lags, c(1,2) means first two lags.
5. Model estimation 5.1. Estimation of the basic models with plm
Several models can be estimated with plm by filling the model argument:
􏰀 the fixed effects model (within), 􏰀 the pooling model (pooling),
􏰀 the first-difference model (fd),
􏰀 the between model (between),
􏰀 the error components model (random).
The basic use of plm is to indicate the model formula, the data and the model to be estimated.
For example, the fixed effects model and the random effects model are estimated using:
R> grun.fe <- plm(inv ~ value + capital, data = Grunfeld, model = "within") R> grun.re <- plm(inv ~ value + capital, data = Grunfeld, model = "random") R> summary(grun.re)
Oneway (individual) effect Random Effect Model (Swamy-Arora’s transformation) Call:
plm(formula = inv ~ value + capital, data = Grunfeld, model = “random”)
Balanced Panel: n=10, T=20, N=200
Effects:
idiosyncratic 2784.458 52.768 0.282
var std.dev share

Journal of Statistical Software 11 individual 7089.800 84.201 0.718
theta: 0.86122
Residuals :
Min. 1st Qu. Median 3rd Qu. Max.
-178.00 -19.70 4.69 19.50 253.00
Coefficients :
Estimate Std. Error t-value Pr(>|t|)
(intercept) -57.834415 28.898935 -2.0013 0.04536 *
value 0.109781 0.010493 10.4627 < 2e-16 *** capital 0.308113 0.017180 17.9339 < 2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Total Sum of Squares: 2381400 Residual Sum of Squares: 548900 Multiple R-Squared: 0.7695 F-statistic: 328.837 on 197 and 2 DF, p-value: 0.0030364 For a random model, the summary method gives information about the variance of the com- ponents of the errors. Fixed effects may be extracted easily using fixef: R> fixef(grun.fe)
123456 -11.552778 160.649753 -176.827902 30.934645 -55.872873 35.582644
7 8 9 10
-7.809534 1.198282 -28.478333 52.176096
The fixef function returns an object of class fixef. A summary method is provided, which prints the effects (in deviation from the overall intercept), their standard errors and the test of equality to the overall intercept.
R> summary(fixef(grun.fe))
Estimate Std. Error t-value Pr(>|t|)
1 -11.55278
2 160.64975
3 -176.82790
48.28725 -0.2393 0.8109
21.97037 7.3121 2.629e-13 ***
21.39348 -8.2655 2.220e-16 ***
7.67856 4.0287 5.608e-05 ***
7.83815 -7.1283 1.016e-12 ***
4.61238 7.7146 1.221e-14 ***
5.07140 -1.5399 0.1236
7.52232 0.1593 0.8734
5.19403 -5.4829 4.184e-08 ***
4 5 6 7 8 9
30.93464
-55.87287
35.58264
-7.80953
1.19828
-28.47833

12 Panel Data Econometrics in R: The plm Package
10 52.17610 0.80774 64.5954 < 2.2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 In case of a two-ways effect model, an additional argument effect is required to extract fixed effects: R> grun.twfe <- plm(inv ~ value + capital, data = Grunfeld, model = "within", + effect = "twoways") R> fixef(grun.twfe, effect = “time”)
1935 1936 1937 1938 1939 1940 1941 1942 1943 1944
72.746 101.607 122.481 77.555 80.526 113.265 139.719 122.665 117.790 120.925
1945 1946 1947 1948 1949 1950 1951 1952 1953 1954
124.159 161.359 147.135 153.948 139.244 151.061 199.583 224.033 275.583 273.781
5.2. More advanced use of plm
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 estima- tions. 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).
For example, to use the amemiya estimator:
R> emp.amem <- plm(log(emp) ~ log(wage) + log(capital), data = EmplUK, + model = "random", random.method = "amemiya") 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"). Journal of Statistical Software 13 For example, to estimate a two-ways effect model for the Grunfeld data: R> grun.tways <- plm(inv ~ value + capital, data = Grunfeld, + effect = "twoways", model = "random", random.method = "amemiya") R> 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.135
individual 8294.716
time 270.529
theta : 0.87475 (id) 0.29695 (time) 0.29595 (total)
Residuals :
Min. 1st Qu. Median 3rd Qu. Max.
-176.00 -18.00 3.02 18.00 233.00
Coefficients :
Estimate Std. Error t-value Pr(>|t|)
(intercept) -64.351811 31.183651 -2.0636 0.03905 *
value 0.111593 0.011028 10.1192 < 2e-16 *** capital 0.324625 0.018850 17.2214 < 2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Total Sum of Squares: 2038000 Residual Sum of Squares: 514120 Multiple R-Squared: 0.74774 F-statistic: 291.965 on 197 and 2 DF, p-value: 0.0034191 In the “effects” section of the result, the variance of the three elements of the error term and the three parameters used in the transformation are now printed. The two-ways effect model is for the moment only available for balanced panels. Instrumental variable estimators One or all of the models may be estimated using instrumental variables. The instruments are specified either as a one side formula in the argument instruments, or at the end of the formula after a | sign. The instrumental variables estimator used is indicated with the inst.method argument: 51.421 0.2359 91.075 0.7400 16.448 0.0241 14 Panel Data Econometrics in R: The plm Package 􏰀 bvk, from Balestra and Varadharajan-Krishnakumar (1987), the default value, 􏰀 baltagi, from Baltagi (1981). We illustrate instrumental variables estimation with the employment equation. We use the first lag of wage as instrument of wage. The four following commands are similar: R> emp.iv <- plm(emp ~ wage + capital | lag(wage, 1) + capital, + data = EmplUK, model = "random") R> emp.iv <- plm(emp ~ wage + capital | . - wage + lag(wage, 1), + data = EmplUK, model = "random") R> emp.iv <- plm(emp ~ wage + capital, instruments = ~lag(wage, 1) + + capital, data = EmplUK, model = "random") R> emp.iv <- plm(emp ~ wage + capital, instruments = ~ . - wage + + lag(wage, 1), data = EmplUK, model = "random") The Hausman-Taylor model (see Hausman and Taylor 1981) may be estimated with plm by equating the model argument to "ht". R> form <- lwage ~ wks + south + smsa + married + exp + I(exp^2) + + bluecol + ind + union + sex + black + ed | sex + black + + bluecol + south + smsa + ind R> ht <- plm(form, data = Wages, model = "ht") R> summary(ht)
Oneway (individual) effect Hausman-Taylor Model
Call:
plm(formula = lwage ~ wks + south + smsa + married + exp + I(exp^2) +
bluecol + ind + union + sex + black + ed, data = Wages, model = “ht”,
instruments = ~sex + black + bluecol + south + smsa + ind)
T.V. exo : bluecol,south,smsa,ind
T.V. endo : wks,married,exp,I(exp^2),union
T.I. exo : sex,black
T.I. endo : ed
Balanced Panel: n=595, T=7, N=4165
Effects:
var std.dev share
idiosyncratic 0.023044 0.151803 0.0253
individual 0.886993 0.941803 0.9747
theta: 0.93919
Residuals :
Min. 1st Qu. Median 3rd Qu. Max.

Journal of Statistical Software 15
-1.92000 -0.07070 0.00657 0.07970 2.03000
Coefficients :
Estimate Std. Error t-value Pr(>|t|)
(intercept) 2.7818e+00 3.0765e-01 9.0422 < 2.2e-16 *** wks 8.3740e-04 5.9973e-04 1.3963 0.16263 southyes 7.4398e-03 3.1955e-02 0.2328 0.81590 smsayes -4.1833e-02 1.8958e-02 -2.2066 0.02734 * marriedyes -2.9851e-02 1.8980e-02 -1.5728 0.11578 exp 1.1313e-01 2.4710e-03 45.7851 < 2.2e-16 *** I(exp^2) -4.1886e-04 5.4598e-05 -7.6718 1.696e-14 *** bluecolyes -2.0705e-02 1.3781e-02 -1.5024 0.13299 ind unionyes sexmale blackyes ed --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 1.3604e-02 1.5237e-02 0.8928 0.37196 3.2771e-02 1.4908e-02 2.1982 0.02794 * 1.3092e-01 1.2666e-01 1.0337 0.30129 -2.8575e-01 1.5570e-01 -1.8352 0.06647 . 1.3794e-01 2.1248e-02 6.4919 8.474e-11 *** Total Sum of Squares: 243.04 Residual Sum of Squares: 95.947 Multiple R-Squared: 0.60522 F-statistic: 489.524 on 4151 and 13 DF, p-value: 3.3651e-16 5.3. Variable coefficients model The pvcm function enables the estimation of variable coefficients models. Time or individual effects are introduced if effect is fixed to "time" or "individual" (the default value). Coefficients are assumed to be fixed if model = "within" or random if model = "random". In the first case, a different model is estimated for each individual (or time period). In the second case, the Swamy model (see Swamy 1970) model is estimated. It is a generalized least squares model which uses the results of the previous model. Denoting βˆi the vectors of coefficients obtained for each individual, we get: 􏰃n􏰄 βˆ = 􏰆 􏰁∆ˆ + σˆi2(Xi⊤Xi)−1􏰂−1 􏰁∆ˆ + σˆi2(Xi⊤Xi)−1􏰂−1 βˆi (8) i=1 where σˆi2 is the unbiased estimator of the variance of the errors for individual i obtained from the preliminary estimation and: 1 n􏰃 1n 􏰄􏰃 1n 􏰄⊤ 1n ∆ˆ = 􏰆 βˆi − 􏰆βˆi βˆi − 􏰆βˆi − 􏰆σˆi2(Xi⊤Xi)−1 (9) n − 1 i=1 n i=1 n i=1 n i=1 If this matrix is not positive-definite, the second term is dropped. With the Grunfeld data, we get: 16 Panel Data Econometrics in R: The plm Package R> grun.varw <- pvcm(inv ~ value + capital, data = Grunfeld, model = "within") R> grun.varr <- pvcm(inv ~ value + capital, data = Grunfeld, model = "random") R> summary(grun.varr)
Oneway (individual) effect Random coefficients model
Call:
pvcm(formula = inv ~ value + capital, data = Grunfeld, model = “random”)
Balanced Panel: n=10, T=20, N=200
Residuals:
Min. 1st Qu. Median Mean 3rd Qu. Max.
-211.500 -32.320 -4.283 9.048 12.710 579.200
Estimated mean of the coefficients:
Estimate Std. Error z-value Pr(>|z|)
(Intercept) -9.629285 17.035040 -0.5653 0.5718946
value 0.084587 0.019956 4.2387 2.248e-05 ***
capital 0.199418 0.052653 3.7874 0.0001522 ***
—
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ‘ 1
Estimated variance of the coefficients:
(Intercept) value capital
(Intercept) 2344.24402 -0.6852340 -4.0276612
value -0.68523 0.0031182 -0.0011847
capital -4.02766 -0.0011847 0.0244824
Total Sum of Squares: 9359900
Residual Sum of Squares: 2194300
Multiple R-Squared: 0.76557
5.4. Generalized method of moments estimator
The generalized method of moments is mainly used in panel data econometrics to estimate dynamic models (Arellano and Bond 1991; Holtz-Eakin, Newey, and Rosen 1988).
yit = ρyit−1 + β⊤xit + μi + εit (10) The model is first differenced to get rid of the individual effect:
∆yit = ρ∆yit−1 + β⊤∆xit + ∆εit (11)
Least squares are inconsistent because ∆εit is correlated with ∆yit−1. yit−2 is a valid, but weak instrument (see Anderson and Hsiao 1981). The GMM estimator uses the fact that the number of valid instruments is growing with t:

differences, one uses:
with:
􏰃n 􏰄−1
􏰆Wi⊤H(1)Wi (14)
i=1
 2 −1 0 … 0 −1 2 −1 … 0

H(1)=d⊤d= 0 −1 2 … 0 (15)
 . . . . .  …..
Journal of Statistical Software 17
􏰀 t = 3: y1,
􏰀 t=4: y1,y2,
􏰀 t=5: y1,y2,y3
For individual i, the matrix of instruments is then:
y1 0 0 0 0 0 … 0 0 0 0 xi3 
0y1y2000…0000xi4  Wi=000y1y2y3…0000xi5 (12)
 . . . . . . . . . .. . .  …………
0 0 0 0 ………y1 y2 …yt−2 xiT−2
The moment conditions are: 􏰅ni=1 Wi⊤ei(β) where ei(β) is the vector of residuals for individual
i. The GMM estimator minimize: 􏰃n􏰄􏰃n􏰄
􏰆 ei(β)⊤Wi A 􏰆 Wi⊤ei(β) (13) i=1 i=1
where A is the weighting matrix of the moments.
One-step estimators are computed using a known weighting matrix. For the model in first
A(1) =
0 0 0 −1 2
Two-steps estimators are obtained using H(2) = 􏰅n e(1)e(1)⊤ where e(1) are the residuals of
i i=1ii i
The GMM estimator is provided by the pgmm function. It’s main argument is a dynformula
the one step estimate.
which describes the variables of the model and the lag structure.
The effect argument is either NULL, “individual” (the default), or “twoways”. In the first case, the model is estimated in levels. In the second case, the model is estimated in first differences to get rid of the individuals effects. In the last case, the model is estimated in first differences and time dummies are included.
In a GMM estimation, there are “normal” instruments and “GMM” instruments. GMM in- struments are indicated with the gmm.inst argument (a one side formula) and the lags with

18 Panel Data Econometrics in R: The plm Package
the lag.gmm argument. By default, all the variables of the model that are not used as GMM
instruments are used as normal instruments, with the same lag structure.
The complete list of instruments can also be specified with the argument instruments which
should be a one side formula (or dynformula).
The model argument specifies whether a one-step or a two-steps model is required (“onestep”
or “twosteps”).
The following example is from Arellano and Bond (1991). Employment is explained by past values of employment (two lags), current and first lag of wages and output and current value of capital.
R> emp.gmm <- pgmm(dynformula(emp ~ wage + capital + output, lag = list(2, + 1, 0, 1), log = TRUE), EmplUK, effect = "twoways", model = "twosteps", + gmm.inst = ~log(emp), lag.gmm = list(c(2, 99))) R> summary(emp.gmm)
Twoways effects Two steps model
Call:
pgmm(formula = log(emp) ~ lag(log(emp), 1) + lag(log(emp), 2) +
log(wage) + lag(log(wage), 1) + log(capital) + log(output) +
lag(log(output), 1), data = EmplUK, effect = “twoways”,
model = “twosteps”, gmm.inst = ~log(emp), lag.gmm = list(c(2, 99)))
Unbalanced Panel: n=140, T=7-9, N=1031
Number of Observations Used: 611
Residuals
Min. 1st Qu. Median Mean 3rd Qu. Max.
-0.6191000 -0.0494800 -0.0004565 -0.0001841 0.0533500 0.6410000
Coefficients
Estimate Std. Error z-value Pr(>|z|)
lag(log(emp), 1) 0.474151
lag(log(emp), 2) -0.052967
log(wage) -0.513205
lag(log(wage), 1) 0.224640
log(capital) 0.292723
log(output) 0.609775
lag(log(output), 1) -0.446373
—
Signif. codes: 0 ‘***’ 0.001
0.085303 5.5584 2.722e-08 ***
0.027284 -1.9413 0.0522200 .
0.049345 -10.4003 < 2.2e-16 *** 0.080063 0.039463 0.108524 0.124815 -3.5763 0.0003485 *** 2.8058 0.0050192 ** 7.4177 1.191e-13 *** 5.6188 1.923e-08 *** '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Sargan Test: chisq(25) = 30.11247 (p.value=0.22011) Autocorrelation test (1): normal = -2.427829 (p.value=0.0075948) Autocorrelation test (2): normal = -0.3325401 (p.value=0.36974) Journal of Statistical Software 19 Wald test for coefficients: chisq(7) = 371.9877 (p.value=< 2.22e-16) Wald test for time dummies: chisq(6) = 26.90450 (p.value=0.00015090) 5.5. General FGLS models General FGLS estimators are based on a two-step estimation process: first an OLS model is estimated, then its residuals uˆit are used to estimate an error covariance matrix more general than the random effects one for use in a feasible-GLS analysis. Formally, the estimated error covariancematrixisVˆ=In⊗Ωˆ,with n uˆituˆ⊤ Ωˆ=􏰆 it i=1 n (see Wooldridge 2002, 10.4.3 and 10.5.5). This framework allows the error covariance structure inside every group (if effect = "individual") of observations to be fully unrestricted and is therefore robust against any type of intragroup heteroskedasticity and serial correlation. This structure, by converse, is assumed identical across groups and thus general FGLS is inefficient under group- wise heteroskedasticity. Cross-sectional correlation is excluded a priori. Moreover, the number of variance parameters to be estimated with N = n × T data points is T (T + 1)/2, which makes these estimators particularly suited for situations where n ≫ T , as e.g., in labour or household income surveys, while problematic for “long” panels, where Vˆ tends to become singular and standard errors therefore become biased downwards. In a pooled time series context (effect = "time"), symmetrically, this estimator is able to account for arbitrary cross-sectional correlation, provided that the latter is time-invariant (see Greene 2003, 13.9.1–2, p.321–322). In this case serial correlation has to be assumed away and the estimator is consistent with respect to the time dimension, keeping n fixed. The function pggls estimates general FGLS models, with either fixed of “random” effects6. The “random effect” general FGLS is estimated by: R> zz <- pggls(log(emp) ~ log(wage) + log(capital), data = EmplUK, + model = "random") R> summary(zz)
Oneway (individual) effect Random effects model
Call:
pggls(formula = log(emp) ~ log(wage) + log(capital), data = EmplUK,
model = “random”)
Unbalanced Panel: n=140, T=7-9, N=1031
Residuals
6The “random effect” is better termed “general FGLS” model, as in fact it does not have a proper random effects structure, but we keep this terminology for general language consistency.

20 Panel Data Econometrics in R: The plm Package Min. 1st Qu. Median Mean 3rd Qu. Max.
-1.8440000 -0.3908000 0.0388700 0.0005525 0.4153000 1.4920000
Coefficients
Estimate Std. Error z-value Pr(>|z|)
(intercept) 1.751448 0.180100 9.7249 < 2e-16 *** log(wage) -0.132986 0.054292 -2.4495 0.01431 * log(capital) 0.629621 0.018483 34.0648 < 2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Total Sum of Squares: 1853.6 Residual Sum of Squares: 386 Multiple R-squared: 0.79176 The fixed effects pggls (see Wooldridge 2002, p. 276) is based on the estimation of a within model in the first step; the rest follows as above. It is estimated by: R> zz <- pggls(log(emp) ~ log(wage) + log(capital), data = EmplUK, + model = "within") The pggls function is similar to plm in many respects. An exception is that the estimate of the group covariance matrix of errors (zz$sigma, a matrix, not shown) is reported in the model objects instead of the usual estimated variances of the two error components. 6. Tests As sketched in Section 2, specification testing in panel models involves essentially testing for poolability, for individual or time unobserved effects and for correlation between these latter and the regressors (Hausman-type tests). As for the other usual diagnostic checks, we provide a suite of serial correlation tests, while not touching on the issue of heteroskedasticity testing. Instead, we provide heteroskedasticity-robust covariance estimators, to be described in Section 6.6. 6.1. Tests of poolability pooltest tests the hypothesis that the same coefficients apply to each individual. It is a standard F test, based on the comparison of a model obtained for the full sample and a model based on the estimation of an equation for each individual. The first argument of pooltest is a plm object. The second argument is a pvcm object obtained with model=within. If the first argument is a pooling model, the test applies to all the coefficients (including the intercepts), if it is a within model, different intercepts are assumed. To test the hypothesis that all the coefficients in the Grunfeld example, excluding the inter- cepts, are equal, we use : R> znp = pvcm(inv ~ value + capital, data = Grunfeld, model = “within”)
R> zplm = plm(inv ~ value + capital, data = Grunfeld)
R> pooltest(zplm, znp)

Journal of Statistical Software 21
F statistic
data: inv ~ value + capital
F = 5.7805, df1 = 18, df2 = 170, p-value = 1.219e-10
alternative hypothesis: unstability
The same test can be computed using a formula as first argument of the pooltest function: R> pooltest(inv ~ value + capital, data = Grunfeld, model = “within”)
6.2. Tests for individual and time effects
plmtest implements Lagrange multiplier tests of individual or/and time effects based on the results of the pooling model. Its main argument is a plm object (the result of a pooling model) or a formula.
Two additional arguments can be added to indicate the kind of test to be computed. The argument type is one of:
􏰀 bp: Breusch and Pagan (1980),
􏰀 honda: Honda (1985), the default value,
􏰀 kw: King and Wu (1997),
􏰀 ghm: Gourieroux, Holly, and Monfort (1982).
The effects tested are indicated with the effect argument (one of individual, time or twoways).
To test the presence of individual and time effects in the Grunfeld example, using the Gourier- oux et al. (1982) test, we use:
R> g <- plm(inv ~ value + capital, data = Grunfeld, model = "pooling") R> plmtest(g, effect = “twoways”, type = “ghm”)
Lagrange Multiplier Test – two-ways effects (Gourieroux, Holly and
Monfort)
data: inv ~ value + capital
chisq = 798.1615, df = 2, p-value < 2.2e-16 alternative hypothesis: significant effects or R> plmtest(inv ~ value + capital, data = Grunfeld, effect = “twoways”,
+ type = “ghm”)

22 Panel Data Econometrics in R: The plm Package
pFtest computes F tests of effects based on the comparison of the within and the pooling models. Its main arguments are either two plm objects (the results of a pooling and a within model) or a formula.
R> gw <- plm(inv ~ value + capital, data = Grunfeld, effect = "twoways", + model = "within") R> gp <- plm(inv ~ value + capital, data = Grunfeld, model = "pooling") R> pFtest(gw, gp)
F test for effects
data: inv ~ value + capital
F = 17.4031, df1 = 28, df2 = 169, p-value < 2.2e-16 alternative hypothesis: significant effects R> pFtest(inv ~ value + capital, data = Grunfeld, effect = “twoways”)
6.3. Hausman test
phtest computes the Hausman test which is based on the comparison of two sets of estimates (see Hausman 1978). Its main arguments are two panelmodel objects or a formula. A classical application of the Hausman test for panel data is to compare the fixed and the random effects models:
R> gw <- plm(inv ~ value + capital, data = Grunfeld, model = "within") R> gr <- plm(inv ~ value + capital, data = Grunfeld, model = "random") R> phtest(gw, gr)
Hausman Test
data: inv ~ value + capital
chisq = 2.3304, df = 2, p-value = 0.3119
alternative hypothesis: one model is inconsistent
6.4. Tests of serial correlation
A model with individual effects has composite errors that are serially correlated by definition. The presence of the time-invariant error component7 gives rise to serial correlation which does not die out over time, thus standard tests applied on pooled data always end up rejecting the null of spherical residuals8. There may also be serial correlation of the “usual” kind in the idiosyncratic error terms, e.g., as an AR(1) process. By “testing for serial correlation” we mean testing for this latter kind of dependence.
7Here we treat fixed and random effects alike, as components of the error term, according with the modern approach in econometrics (see Wooldridge 2002).
8Neglecting time effects may also lead to serial correlation in residuals (as observed in Wooldridge 2002, 10.4.1).

Journal of Statistical Software 23
For these reasons, the subjects of testing for individual error components and for serially correlated idiosyncratic errors are closely related. In particular, simple (marginal) tests for one direction of departure from the hypothesis of spherical errors usually have power against the other one: in case it is present, they are substantially biased towards rejection. Joint tests are correctly sized and have power against both directions, but usually do not give any information about which one actually caused rejection. Conditional tests for serial correlation that take into account the error components are correctly sized under presence of both departures from sphericity and have power only against the alternative of interest. While most powerful if correctly specified, the latter, based on the likelihood framework, are crucially dependent on normality and homoskedasticity of the errors.
In plm we provide a number of joint, marginal and conditional ML-based tests, plus some semiparametric alternatives which are robust versus heteroskedasticity and free from distri- butional assumptions.
Unobserved effects test
The unobserved effects test `a la Wooldridge (see Wooldridge 2002, 10.4.4), is a semiparametric
test for the null hypothesis that σμ2 = 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 􏰅n 􏰅T −1 􏰅T uˆ uˆ (where uˆ are the i=1 t=1 s=t+1 it is
pooled OLS residuals), which must be “statistically close” to zero under the null, scaled by
its standard deviation:
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):
R> 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 = 8.207e-05
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.
􏰅n􏰅T−1􏰅T uˆuˆ i=1 t=1 s=t+1 it is
W=[􏰅n (􏰅T−1􏰅T uˆuˆ)2]1/2 i=1 t=1 s=t+1 it is

24 Panel Data Econometrics in R: The plm Package
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:
R> pbsytest(log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp, data = Produc,
+ test = “J”)
Baltagi and Li AR-RE joint test
data: log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp
chisq = 4187.597, df = 2, p-value < 2.2e-16 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) 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 test9. 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: R> pbsytest(log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp, data = Produc)
Bera, Sosa-Escudero and Yoon locally robust test
data: log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp
chisq = 52.6359, df = 1, p-value = 4.015e-13
alternative hypothesis: AR(1) errors sub random effects
The BSY test for random effects is implemented in the one-sided version10, which takes heed that the variance of the random effect must be non-negative:
9LM3 in Baltagi and Li (1995). 10Corresponding to RSOμ∗ in the original paper.

Journal of Statistical Software 25 R> pbsytest(log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp, data = Produc,
+ test = “RE”)
Bera, Sosa-Escudero and Yoon locally robust test
data: log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp
z = 57.9143, p-value < 2.2e-16 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) ran- dom 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: R> 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 < 2.2e-16 alternative hypothesis: AR(1)/MA(1) errors in RE panel models 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 any more, 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 7). 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 3). 26 Panel Data Econometrics in R: The plm Package For the random effects model, Wooldridge (2002) observes that under the null of homoskedas- ticity 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 com- ponent. 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(uˆit,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”11. 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 pack- age lmtest: R> grun.fe <- plm(inv ~ value + capital, data = Grunfeld, model = "within") R> pbgtest(grun.fe, order = 2)
Breusch-Godfrey/Wooldridge test for serial correlation in
panel models
data: inv ~ value + capital
chisq = 42.5867, df = 2, p-value = 5.655e-10
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 ‘. . . ’ operator.
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(εˆ , εˆ ) = −1/(T − 1)
11Baltagi and Li derive a basically analogous T-asymptotic test for first-order serial correlation in a FE panel model as a Breusch-Godfrey LM test on within residuals (see Baltagi and Li 1995, Paragraph 2.3 and Formula 12). They also observe that the test on within residuals can be used for testing on the RE model, as “the within transformation [time-demeaning, in our terminology] wipes out the individual effects, whether fixed or random”. Generalizing the Durbin-Watson test to FE models by applying it to fixed effects residuals is documented in Bhargava, Franzini, and Narendranathan (1982).
it is

Journal of Statistical Software 27 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: εˆ=α+δεˆ +η
i,t i,t−1 i,t
Rejecting the restriction δ = −1/(T − 1) makes us conclude against the original null of no
serial correlation.
The building blocks available in plm, together with the function linear.hypothesis() in package car, 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 linear.hypothesis() to the latter model to test the above restriction on δ, supplying a heteroskedasticity- and autocorrelation-consistent covariance matrix (pvcovHC with the appropriate options, in par- ticular method = “arellano”, see Section 6.6).
R> pwartest(log(emp) ~ log(wage) + log(capital), data = EmplUK)
Wooldridge’s test for serial correlation in FE panels
data: log(emp) ~ log(wage) + log(capital)
chisq = 312.2975, p-value < 2.2e-16 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 corre- lation 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 uit are uncorrelated, the errors of the (first) differenced model12 eit ≡ uit − ui,t−1 will be correlated, with cor(eit, ei,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 uˆi,t = δuˆi,t−1 + ηi,t and testing the restriction δ = −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 eit are uncorre- lated, as by definition uit = ui,t−1 + eit, then uit 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: 12Here, eit for notational simplicity (and as in Wooldridge): equivalent to ∆εit in the general notation of the paper. 28 Panel Data Econometrics in R: The plm Package R> pwfdtest(log(emp) ~ log(wage) + log(capital), data = EmplUK)
Wooldridge’s first-difference test for serial correlation in panels
data: log(emp) ~ log(wage) + log(capital)
Chisq = 0.9316, p-value = 0.3344
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.
R> pwfdtest(log(emp) ~ log(wage) + log(capital), data = EmplUK,
+ h0 = “fe”)
Wooldridge’s first-difference test for serial correlation in panels
data: log(emp) ~ log(wage) + log(capital)
Chisq = 131.011, p-value < 2.2e-16 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 Section 6.6 in inference. 6.5. Tests for cross-sectional dependence Next to the more familiar issue of serial correlation, over the last years a growing body of literature has been dealing with cross-sectional dependence (henceforth: XSD) in panels, which can arise, e.g., if individuals respond to common shocks (as in the literature on factor models) or if spatial diffusion processes are present, relating individuals in a way depending on a measure of distance (spatial models). The subject is huge, and here we touch only some general aspects of misspecification testing and valid inference. If XSD is present, the consequence is, at a minimum, inefficiency of the usual estimators and invalid inference when using the standard covariance matrix13.The plan is to have in plm both misspecification tests to detect XSD and robust covariance matrices to perform valid inference in its presence, like in the serial dependence case. For now, though, only misspecification tests are included. 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 13This is the case, e.g., if in an unobserved effects model when XSD is due to an unobservable factor structure, with factors that are uncorrelated with the regressors. In this case the within or random estimators are still consistent, although inefficient (see De Hoyos and Sarafidis 2006). Journal of Statistical Software 29 dimension and a few time periods. All are based on (transformations of) the product-moment correlation coefficient of a model’s residuals, defined as 􏰅Tt=1 uˆituˆjt ρˆij = (􏰅Tt=1 uˆ2it)1/2(􏰅Tt=1 uˆ2jt)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 ρij, is valid for T → ∞ with n fixed; defined as LM = n−1 n 􏰆􏰆 i=1 j=i+1 2 where in the case of an unbalanced panel only pairwise complete observations are considered, and Tij = min(Ti,Tj) with Ti being the number of observations for individual i; else, if the panel is balanced, Tij = T for each i, j. The test is distributed as χ2n(n−1)/2. It is inappropriate whenever the n dimension is “large”. A scaled version, applicable also if T → ∞ and then n → ∞ (as in some pooled time series contexts), is defined as  n−1 n 􏰇 􏰆􏰆 i=1 j=i+1 􏰈 SCLM= n(n−1) and distributed as a standard Normal. Pesaran’s (Pesaran 2004) CD test Tijρˆij 1 2 Tijρˆij  􏰈 2 n−1 n 􏰇  CD= 􏰆􏰆 Tijρˆij n(n − 1)  i=1 j=i+1 based on ρ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". 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. R> pcdtest(inv ~ value + capital, data = Grunfeld)
Pesaran CD test for cross-sectional dependence in panels
data: formula
z = 5.3401, p-value = 9.292e-08
alternative hypothesis: cross-sectional dependence

30 Panel Data Econometrics in R: The plm Package
If a different model specification (within, random, …) is assumed consistent, one can resort to its residuals for testing14 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.
R> pcdtest(inv ~ value + capital, data = Grunfeld, model = “within”)
Pesaran CD test for cross-sectional dependence in panels
data: formula
z = 4.6612, p-value = 3.144e-06
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
􏰈
1
􏰅n−1 􏰅n w(p)
i=1 j=i+1 ij

n−1n 􏰇 􏰆􏰆
 [w(p)]ij Tijρˆij i=1 j=i+1
CD =
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 ρˆ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 and SCLM tests.
Therefore, the local version of either test can be computed supplying an n × 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
14This is also the only solution when the time dimension’s length is insufficient for estimating the heteroge- neous model.

Journal of Statistical Software 31 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 w15. 6.6. Robust covariance matrix estimation
Robust estimators of the covariance matrix of coefficients are provided, mostly for use in Wald-type tests. pvcovHC estimates three “flavours” of White’s heteroskedasticity-consistent covariance matrix16 (known as the sandwich estimator). Interestingly, in the context of panel data the most general version also proves consistent versus serial correlation.
All types assume no correlation between errors of different groups while allowing for het- eroskedasticity across groups, so that the full covariance matrix of errors is V = In ⊗ Ωi; i = 1, .., n. As for the intragroup error covariance matrix of every single group of observations, “white1” allows for general heteroskedasticity but no serial correlation, i.e.,
σ2 … … 0  i1
 2 .
 0 σi2 . 
Ωi= . ..  (16)
 . . 0  0 … … σ2
while “white2” is “white1” restricted to a common variance inside every group, estimated as σi2 = 􏰅Tt=1 uˆ2it/T, so that Ωi = IT ⊗σi2 (see Greene (2003, 13.7.1–2) and Wooldridge (2002, 10.7.2); “arellano” (see ibid. and the original ref. Arellano 1987) allows a fully general structure w.r.t. heteroskedasticity and serial correlation:
 σ2 σi1,i2 … … σi1,iT  i1
 σi2,i1 σ2 .  i2
iT
….. (17)  . 
Ωi = 
 . σiT−1 σiT−1,iT
.2
σiT,i1 … … σiT,iT−1 σ2
iT
The latter is, as already observed, consistent w.r.t. timewise correlation of the errors, but on the converse, unlike the White 1 and 2 methods, it relies on large n asymptotics with small T.
The fixed effects case, as already observed in Section 6.4 on serial correlation, is complicated by the fact that the demeaning induces serial correlation in the errors. The original White estimator (white1) turns out to be inconsistent for fixed T as n grows, so in this case it is advisable to use the arellano version (see Stock and Watson 2006).
The errors may be weighted according to the schemes proposed by MacKinnon and White (1985) and Cribari-Neto (2004) to improve small-sample performance17.
15The very comprehensive package spdep for spatial dependence analysis (see Bivand 2008) contains features for creating, lagging and manipulating neighbour list objects of class nb, that can be readily converted to and from proximity matrices by means of the nb2mat function. Higher orders of the CD(p) test can be obtained lagging the corresponding nbs through nblag.
16See White (1980) and White (1984).
17The HC3 and HC4 weighting schemes are computationally expensive and may hit memory limits for nT in the thousands, where on the other hand it makes little sense to apply small sample corrections.

32 Panel Data Econometrics in R: The plm Package
The main use of pvcovHC is together with testing functions from the lmtest and car packages. These typically allow passing the vcov parameter either as a matrix or as a function (see Zeileis 2004). If one is happy with the defaults, it is easiest to pass the function itself:
R> library(“lmtest”)
R> re <- plm(inv ~ value + capital, data = Grunfeld, model = "random") R> coeftest(re, vcov = pvcovHC)
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
(intercept) -57.834415 23.449626 -2.4663 0.01451 *
value 0.109781 0.012984 8.4551 6.186e-15 *** capital 0.308113 0.051889 5.9379 1.284e-08 ***
—
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ‘ 1
else one may do the covariance computation inside the call to coeftest, thus passing on a matrix:
R> coeftest(re, vcov = pvcovHC(re, method = “white2”, type = “HC3”))
For some tests, e.g., for multiple model comparisons by waldtest, one should always provide a function18. In this case, optional parameters are provided as shown below (see also Zeileis 2004, p.12):
R> waldtest(re, update(re, . ~ . – capital),
+ vcov = function(x) pvcovHC(x, method = “white2”, type = “HC3”))
Wald test
Model 1: inv ~ value + capital
Model 2: inv ~ value
Res.Df Df Chisq Pr(>Chisq)
1 197
2 198 -1 87.828 < 2.2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Moreover, linear.hypothesis from package car may be used to test for linear restrictions: R> library(“car”)
R> linear.hypothesis(re, “2 * value = capital”, vcov = pvcovHC)
18Joint zero-restriction testing still allows providing the vcov of the unrestricted model as a matrix, see the documentation of package lmtest.

Journal of Statistical Software 33
Linear hypothesis test
Hypothesis:
2 value – capital = 0
Model 1: inv ~ value + capital
Model 2: restricted model
Note: Coefficient covariance matrix supplied.
Res.Df Df Chisq Pr(>Chisq)
1 197
2 198 -1 3.4783 0.06218 .
—
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ‘ 1
A specific pvcovHC method for pgmm objects is also provided which implements the robust covariance matrix proposed by Windmeijer (2005) for generalized method of moments esti- mators.
7. plm versus nlme/lme4
The models termed panel by the econometricians have counterparts in the statistics literature on mixed models (or hierarchical models, or models for longitudinal data), although there are both differences in jargon and more substantial distinctions. This language inconsistency between the two communities, together with the more complicated general structure of sta- tistical models for longitudinal data and the associated notation in the software, is likely to scare some practicing econometricians away from some potentially useful features of the R environment, so it may be useful to provide here a brief reconciliation between the typical panel data specifications used in econometrics and the general framework used in statistics for mixed models19.
R is particularly strong on mixed models’ estimation, thanks to the long-standing nlme pack- age (see Pinheiro et al. 2007) and the more recent lme4 package, based on S4 classes (see Bates 2007)20. In the following we will refer to the more established nlme to give some ex- amples of “econometric” panel models that can be estimated in a likelihood framework, also including some likelihood ratio tests. Some of them are not feasible in plm and make a useful complement to the econometric “toolbox” available in R.
7.1. Fundamental differences between the two approaches
Econometrics deal mostly with non-experimental data. Great emphasis is put on specifica- tion procedures and misspecification testing. Model specifications tend therefore to be very
19This discussion does not consider GMM models. One of the basic reasons for econometricians not to choose maximum likelihood methods in estimation is that the strict exogeneity of regressors assumption required for consistency of the ML models reported in the following is often inappropriate in economic settings.
20The standard reference on the subject of mixed models in S/R is Pinheiro and Bates (2000).

34 Panel Data Econometrics in R: The plm Package
simple, while great attention is put on the issues of endogeneity of the regressors, dependence structures in the errors and robustness of the estimators under deviations from normality. The preferred approach is often semi- or non-parametric, and heteroskedasticity-consistent techniques are becoming standard practice both in estimation and testing.
For all these reasons, although the maximum likelihood framework is important in testing21 and sometimes used in estimation as well, panel model estimation in econometrics is mostly accomplished in the generalized least squares framework based on Aitken’s Theorem and, when possible, in its special case OLS, which are free from distributional assumptions (al- though these kick in at the diagnostic testing stage). On the contrary, longitudinal data models in nlme and lme4 are estimated by (restricted or unrestricted) maximum likelihood. While under normality, homoskedasticity and no serial correlation of the errors OLS are also the maximum likelihood estimator, in all the other cases there are important differences.
The econometric GLS approach has closed-form analytical solutions computable by stan- dard linear algebra and, although the latter can sometimes get computationally heavy on the machine, the expressions for the estimators are usually rather simple. ML estimation of longitudinal models, on the contrary, is based on numerical optimization of nonlinear func- tions without closed-form solutions and is thus dependent on approximations and convergence criteria. For example, the “GLS” functionality in nlme is rather different from its “economet- ric” counterpart. “Feasible GLS” estimation in plm is based on a single two-step procedure, in which an inefficient but consistent estimation method (typically OLS) is employed first in order to get a consistent estimate of the errors’ covariance matrix, to be used in GLS at the second step; on the converse, “GLS” estimators in nlme are based on iteration until convergence of two-step optimization of the relevant likelihood.
7.2. Some false friends
The fixed/random effects terminology in econometrics is often recognized to be misleading, as both are treated as random variates in modern econometrics (see e.g. Wooldridge 2002, 10.2.1). It has been recognized since Mundlak’s classic paper (Mundlak 1978) that the fundamental issue is whether the unobserved effects are correlated with the regressors or not. In this last case, they can safely be left in the error term, and the serial correlation they induce is cared for by means of appropriate GLS transformations. On the contrary, in the case of correlation, “fixed effects” methods such as least squares dummy variables or time-demeaning are needed, which explicitly, although inconsistently22, estimate a group- (or time-) invariant additional parameter for each group (or time period).
Thus, from the point of view of model specification, having fixed effects in an econometric model has the meaning of allowing the intercept to vary with group, or time, or both, while the other parameters are generally still assumed to be homogeneous. Having random effects means having a group- (or time-, or both) specific component in the error term.
In the mixed models literature, on the contrary, fixed effect indicates a parameter that is
21Lagrange Multiplier tests based on the likelihood principle are suitable for testing against more general alternatives on the basis of a maintained model with spherical residuals and find therefore application in testing for departures from the classical hypotheses on the error term. The seminal reference is Breusch and Pagan (1980).
22For fixed effects estimation, as the sample grows (on the dimension on which the fixed effects are specified) so does the number of parameters to be estimated. Estimation of individual fixed effects is T- (but not n-) consistent, and the opposite.

Journal of Statistical Software 35 assumed constant, while random effects are parameters that vary randomly around zero ac-
cording to a joint multivariate Normal distribution.
So, the FE model in econometrics has no counterpart in the mixed models framework, unless reducing it to OLS on a specification with one dummy for each group (often termed least squares dummy variables, or LSDV model) which can trivially be estimated by OLS. The RE model is instead a special case of mixed model where only the intercept is specified as a random effect, while the “random” type variable coefficients model can be seen as one that has the same regressors in the fixed and random sets. The unrestricted generalized least squares can in turn be seen, in the nlme framework, as a standard linear model with a general error covariance structure within the groups and errors uncorrelated across groups.
7.3. A common taxonomy
To reconcile the two terminologies, in the following we report the specification of the panel models in plm according to the general expression of a mixed model in Laird-Ware form (see the web appendix to Fox 2002) and the nlme estimation commands for maximum likelihood estimation of an equivalent specification23.
The Laird-Ware representation for mixed models
A general representation for the linear mixed effects model is given in Laird and Ware (1982).
yit = β1x1ij +…+βpxpij b1z1ij +…+bpzpij +εij
bik ∼ N(0, ψk2), Cov(bk, bk′ ) = ψkk′
εij ∼ N(0,σ2λijj), Cov(εij,εij′) = σ2λijj′
where the x1, . . . xp are the fixed effects regressors and the z1, . . . zp are the random effects regressors, assumed to be normally distributed across groups. The covariance of the random effects coefficients ψkk′ is assumed constant across groups and the covariances between the errorsingroupi,σ2λijj′,aredescribedbythetermλijj′ representingthecorrelationstructure of the errors within each group (e.g., serial correlation over time) scaled by the common error variance σ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: biq = 0 ∀i,q, λijj = σ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.
23In doing so, we stress that “equivalence” concerns only the specification of the model, and neither the appropriateness nor the relative efficiency of the relevant estimation techniques, which will of course be depen- dent on the context. Unlike their mixed model counterparts, the specifications in plm are, strictly speaking, distribution-free. Nevertheless, for the sake of exposition, in the following we present them in the setting which ensures consistency and efficiency (e.g., we consider the hypothesis of spherical errors part of the specification of pooled OLS and so forth).

36 Panel Data Econometrics in R: The plm Package Random effects
In the Laird and Ware notation, the RE specification is a model with only one random effects regressor: the intercept. Formally, z1ij = 1 ∀i, j, zqij = 0 ∀i, ∀j, ∀q ̸= 1 λij = 1 for i = j, 0 else). The composite error is therefore uij = 1bi1 + εij. Below we report coefficients of Grunfeld’s model estimated by GLS and then by ML
R> data(Grunfeld, package = “Ecdat”)
R> reGLS <- plm(inv ~ value + capital, data = Grunfeld, model = "random") R> reML <- lme(inv ~ value + capital, data = Grunfeld, random = ~1 | firm) R> coef(reGLS)
(intercept) value capital
-57.8344149 0.1097812 0.3081130
R> summary(reML)$coef$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 zq = xq ∀q, i.e., the fixed effects and the random effects regressors are the same, and ψkk′ = σμ2IN, and λijj = 1, λijj′ = 0 for j ̸= 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.
R> vcm <- pvcm(inv ~ value + capital, data = Grunfeld, model = "random", + effect = "time") R> vcmML <- lme(inv ~ value + capital, data = Grunfeld, + random = ~ value + capital | year) R> coef(vcm)
(Intercept) value capital
-18.5538638 0.1239595 0.1114579
R> summary(vcmML)$coef$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

Journal of Statistical Software 37 is done by creating an lmList object, so that the two models below are equivalent (output
suppressed):
R> vcmf <- pvcm(inv ~ value + capital, data = Grunfeld, model = "within", + effect = "time") R> vcmfML <- lmList(inv ~ value + capital | year, data = Grunfeld) Unrestricted FGLS The general, or unrestricted, feasible GLS, pggls in the plm nomenclature, is equivalent to a model with no random effects regressors (biq = 0 ∀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. R> sGrunfeld <- Grunfeld[Grunfeld$firm %in% 4:6, ] R> ggls <- pggls(inv ~ value + capital, data = sGrunfeld, model = "random") R> gglsML <- gls(inv ~ value + capital, data = sGrunfeld, + correlation = corSymm(form = ~ 1 | year)) R> coef(ggls)
(intercept) value capital
1.19679342 0.10555908 0.06600166
R> summary(gglsML)$coef
(Intercept) value capital
-2.4156266 0.1163550 0.0735837
The within case is analogous, with the regressors’ set augmented by n − 1 group dummies. 7.4. Some useful “econometric” models in nlme
Finally, amongst the many possible specifications estimable with nlme, we report a couple cases that might be especially interesting to applied econometricians.
AR(1) pooling or random effects panel
Linear models with groupwise structures of time-dependence24 may be fitted by gls(), spec- ifying the correlation structure in the correlation option:
R> lmAR1ML <- gls(inv ~ value + capital, data = Grunfeld, + correlation = corAR1(0, form = ~ year | firm)) 24Take heed that here, in contrast to the usual meaning of serial correlation in time series, we always speak of serial correlation between the errors of each group. 38 Panel Data Econometrics in R: The plm Package and analogously the random effects panel with, e.g., AR(1) errors (see Baltagi 2001, Chap- ter 5), which is a very common specification in econometrics, may be fit by lme specifying an additional random intercept: R> 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: R> summary(reAR1ML)$coef$fixed
(Intercept) value capital
-40.27650822 0.09336672 0.31323330
R> 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().
R> lmML <- gls(inv ~ value + capital, data = Grunfeld) R> 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: R> anova(reML, reAR1ML)

Journal of Statistical Software 39
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: R> 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.3660 <.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: R> 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.
8. Conclusions
With plm we aim at providing a comprehensive package containing the standard function- alities that are needed for the management and the econometric analysis of panel data. In particular, we provide: functions for data transformation; estimators for pooled, random and fixed effects static panel models and variable coefficients models, general GLS for general covariance structures, and generalized method of moments estimators for dynamic panels; specification and diagnostic tests. Instrumental variables estimation is supported. Most es- timators allow working with unbalanced panels. While among the different approaches to longitudinal data analysis we take the perspective of the econometrician, the syntax is con- sistent with the basic linear modeling tools, like the lm function.
On the input side, formula and data arguments are used to specify the model to be estimated. Special functions are provided to make writing formulas easier, and the structure of the data is indicated with an index argument.
On the output side, the model objects (of the new class panelmodel) are compatible with the general restriction testing frameworks of packages lmtest and car. Specialized methods are also provided for the calculation of robust covariance matrices; heteroskedasticity- and correlation-consistent testing is accomplished by passing these on to testing functions, together with a panelmodel object.

40 Panel Data Econometrics in R: The plm Package
The main functionalities of the package have been illustrated here by applying them on some well-known datasets from the econometric literature. The similarities and differences with the maximum likelihood approach to longitudinal data have also been briefly discussed.
We plan to expand the methods in this paper to systems of equations and to the estimation of models with autoregressive errors. Addition of covariance estimators robust versus cross- sectional correlation are also in the offing. Lastly, conditional visualization features in the R environment seem to offer a promising toolbox for visual diagnostics, which is another subject for future work.
Acknowledgments
While retaining responsibility for any error, we thank Jeffrey Wooldridge, Achim Zeileis and three anonymous referees for useful comments. We also acknowledge kind editing assistance by Lisa Benedetti.
References
Amemiya T (1971). “The Estimation of the Variances in a Variance-Components Model.” International Economic Review, 12, 1–13.
Anderson TW, Hsiao C (1981). “Estimation of Dynamic Models with Error Components.” Journal of the American Statistical Association, 76, 598–606.
Arellano M (1987). “Computing Robust Standard Errors for Within Group Estimators.” Oxford Bulletin of Economics and Statistics, 49, 431–434.
Arellano M, Bond S (1991). “Some Tests of Specification for Panel Data: Monte Carlo Evidence and an Application to Employment Equations.” Review of Economic Studies, 58, 277–297.
Atkinson B, Therneau T (2007). kinship: Mixed-Effects Cox Models, Sparse Matri- ces, and Modeling Data from Large Pedigrees. R package version 1.1.0-18, URL http: //CRAN.R-project.org/package=kinship.
Balestra P, Varadharajan-Krishnakumar J (1987). “Full Information Estimations of a System of Simultaneous Equations with Error Components.” Econometric Theory, 3, 223–246.
Baltagi BH (1981). “Simultaneous Equations with Error Components.” Journal of Economet- rics, 17, 21–49.
Baltagi BH (2001). Econometric Analysis of Panel Data. 3rd edition. John Wiley & Sons, Hoboken, NJ.
Baltagi BH, Li Q (1991). “A Joint Test for Serial Correlation and Random Individual Effects.” Statistics and Probability Letters, 11, 277–280.
Baltagi BH, Li Q (1995). “Testing AR(1) Against MA(1) Disturbances in an Error Component Model.” Journal of Econometrics, 68, 133–151.

Journal of Statistical Software 41 Bates D (2004). “Least Squares Calculations in R.” R News, 4(1), 17–20. URL http:
//CRAN.R-project.org/doc/Rnews/.
Bates D (2007). lme4: Linear Mixed-Effects Models Using S4 Classes. R package ver-
sion 0.99875-9, URL http://CRAN.R-project.org/package=lme4.
Bates D, Maechler M (2007). matrix: A Matrix Package for R. R package version 0.99875-2,
URL http://CRAN.R-project.org/package=Matrix.
Bera AK, Sosa-Escudero W, Yoon M (2001). “Tests for the Error Component Model in the
Presence of Local Misspecification.” Journal of Econometrics, 101, 1–23.
Bhargava A, Franzini L, Narendranathan W (1982). “Serial Correlation and the Fixed Effects
Model.” Review of Economic Studies, 49, 533–554.
Bivand R (2008). spdep: Spatial Dependence: Weighting Schemes, Statistics and Models.
R package version 0.4-17, URL http://CRAN.R-project.org/package=spdep.
Breusch TS, Pagan AR (1980). “The Lagrange Multiplier Test and its Applications to Model
Specification in Econometrics.” Review of Economic Studies, 47, 239–253.
Cornwell C, Rupert P (1988). “Efficient Estimation with Panel Data: An Empirical Compar-
ison of Instrumental Variables Estimators.” Journal of Applied Econometrics, 3, 149–155. Cribari-Neto F (2004). “Asymptotic Inference Under Heteroskedasticity of Unknown Form.”
Computational Statistics & Data Analysis, 45, 215–233.
Croissant Y (2006). Ecdat: Data Sets for Econometrics. R package version 0.1-5, URL
http://CRAN.R-project.org/package=Ecdat.
De Hoyos RE, Sarafidis V (2006). “Testing for Cross-Sectional Dependence in Panel-Data
Models.” The Stata Journal, 6(4), 482–496.
Drukker DM (2003). “Testing for Serial Correlation in Linear Panel-Data Models.” The Stata
Journal, 3(2), 168–177.
Fox J (2002). An R and S-PLUS Companion to Applied Regression. Sage Publications,
Thousand Oaks, CA.
Fox J (2007). car: Companion to Applied Regression. R package version 1.2-5, URL http:
//CRAN.R-project.org/package=car.
Gourieroux C, Holly A, Monfort A (1982). “Likelihood Ratio Test, Wald Test, and Kuhn- Tucker Test in Linear Models with Inequality Constraints on the Regression Parameters.” Econometrica, 50, 63–80.
Greene WH (2003). Econometric Analysis. 5th edition. Prentice Hall.
Hausman JA (1978). “Specification Tests in Econometrics.” Econometrica, 46, 1251–1271.
Hausman JA, Taylor WE (1981). “Panel Data and Unobservable Individual Effects.” Econo- metrica, 49, 1377–1398.

42 Panel Data Econometrics in R: The plm Package
Holtz-Eakin D, Newey W, Rosen HS (1988). “Estimating Vector Autoregressions with Panel
Data.” Econometrica, 56, 1371–1395.
Honda Y (1985). “Testing the Error Components Model with Non-Normal Disturbances.”
Review of Economic Studies, 52, 681–690.
King ML, Wu PX (1997). “Locally Optimal One-Sided Tests for Multiparameter Hypotheses.”
Econometric Reviews, 33, 523–529.
Koenker R, Ng P (2007). sparsem: Sparse Linear Algebra. R package version 0.74, URL
http://CRAN.R-project.org/package=sparsem.
Laird NM, Ware JH (1982). “Random-Effects Models for Longitudinal Data.” Biometrics,
38, 963–974.
MacKinnon JG, White H (1985). “Some Heteroskedasticity-Consistent Covariance Matrix Estimators with Improved Finite Sample Properties.” Journal of Econometrics, 29, 305– 325.
Mundlak Y (1978). “On the Pooling of Time Series and Cross Section Data.” Econometrica, 46(1), 69–85.
Munnell AH (1990). “Why Has Productivity Growth Declined? Productivity and Public Investment.” New England Economic Review, pp. 3–22.
Nerlove M (1971). “Further Evidence on the Estimation of Dynamic Economic Relations from a Time-Series of Cross-Sections.” Econometrica, 39, 359–382.
Pesaran MH (2004). “General Diagnostic Tests for Cross Section Dependence in Panels.” Technical Report 1229, CESifo Working Paper Series.
Pinheiro J, Bates D, DebRoy S, Sarkar D (2007). nlme: Linear and Nonlinear Mixed Effects Models. R package version 3.1-86, URL http://CRAN.R-project.org/package=nlme.
Pinheiro JC, Bates D (2000). Mixed-Effects Models in S and S-PLUS. Springer-Verlag, New York.
R Development Core Team (2008). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0, URL http: //www.R-project.org/.
Stock J, Watson M (2006). “Heteroskedasticity-Robust Standard Errors for Fixed Effects Panel Data Regression.” Technical Report 0323, NBER Working Paper Series.
Swamy PAVB (1970). “Efficient Inference in a Random Coefficient Regression Model.” Econo- metrica, 38, 311–323.
Swamy PAVB, Arora SS (1972). “The Exact Finite Sample Properties of the Estimators of Coefficients in the Error Components Regression Models.” Econometrica, 40, 261–275.
Wallace TD, Hussain A (1969). “The Use of Error Components Models in Combining Cross Section with Time Series Data.” Econometrica, 37(1), 55–72.

Journal of Statistical Software 43 White H (1980). Asymptotic Theory for Econometricians. Academic Press, Orlando.
White H (1984). “A Heteroskedasticity-Consistent Covariance Matrix and a Direct Test for Heteroskedasticity.” Econometrica, 48, 817–838.
Windmeijer F (2005). “A Finite Sample Correction for the Variance of Linear Efficient Two- Steps GMM Estimators.” Journal of Econometrics, 126, 25–51.
Wooldridge JM (2002). Econometric Analysis of Cross-Section and Panel Data. MIT Press. Zeileis A (2004). “Econometric Computing with HC and HAC Covariance Matrix Estimators.”
Journal of Statistical Software, 11(10), 1–17. URL http://www.jstatsoft.org/v11/i10/. Zeileis A, Hothorn T (2002). “Diagnostic Checking in Regression Relationships.” R News,
2(3), 7–10. URL http://CRAN.R-project.org/doc/Rnews/.
Affiliation:
Yves Croissant
LET-ISH
Avenue Berthelot
F-69363 Lyon cedex 07, France Telephone: +33/4/78727249 Fax: +33/4/78727248
E-mail: yves.croissant@let.ish-lyon.cnrs.fr
Giovanni Millo
DiSES, Un. of Trieste and R&D Dept., Generali SpA Via Machiavelli 4
I-34131 Trieste, Italy
Telephone: +39/040/671184
Fax: +39/040/671160
E-mail: Giovanni_Millo@Generali.com
Journal of Statistical Software
published by the American Statistical Association
Volume 27, Issue 2 July 2008
http://www.jstatsoft.org/
http://www.amstat.org/
Submitted: 2007-06-06 Accepted: 2008-03-20