The formulation and estimation of random effects panel data models of trade

Munich Personal RePEc Archive

The formulation and estimation of random effects panel data models of trade

Matyas, Laszlo and Hornok, Cecilia and Pus, Daria

Central European University

17 February 2012

Online at https://mpra.ub.uni-muenchen.de/36789/

MPRA Paper No. 36789, posted 20 Feb 2012 13:47 UTC

Working Paper 2012/2 Central European University Department of Economics

The Formulation and Estimation of Random Effects Panel Data Models of Trade

Laszlo Matyas

¹

Central European University Cecilia Hornok

Central European University and

Daria Pus

Central European University

February 16, 2012

Abstract: The paper introduces for the most frequently used three-dimensional panel data sets several random effects model specifications. It derives appropriate estimation methods for the balanced and unbalanced cases. An application is also presented where the bilateral trade of 20 EU countries is analysed for the period 2001-2006. The differences between the fixed and random effects specifications are highlighted through this empirical exercise.

Key words: panel data, multidimensional panel data, random effects, error compo- nents model, trade model, gravity model.

JEL classification: C1, C2, C4, F17, F47.

Acknowledgement: Support by Australian Research Council grant DP110103824 is kindly acknowledged, as well as the hospitality of Monash University for Laszlo Matyas.

1 Corresponding author; Central European University, Department of Economics, Bu- dapest 1052, Nador u. 9, Hungary; matyas@ceu.hu

1. Introduction

The use of multidimensional panel data sets has received momentum the last few years. Especially, three dimensional data bases are becoming readily available and frequently used to analyze different types of economic flows, like capital flows (FDI) for example, or most predominantly trade relationships (for a recent reviews of the subject see Anderson [2010] or van Bergeijk and Brakman [2010]). Several model specifications have been proposed in the literature to deal with the heterogeneity of these types of data sets, but all of them considered these heterogeneity factors as fixed effects, i.e., fixed unknown parameters. As it is pretty well understood from the use of “usual” two dimensional panel data sets, the fixed effects formulations are more suited to deal with cases when the panel, at least in one dimension, is short.

On the other hand, for large data sets, the random effects specifications seems to be more suited, where the specific effects are considered as random variables, rather than parameters.

In this paper we present different types of random effects model specifications which mirror the fixed effects models used so far in the literature (some earlier versions were introduced in Davis[2002]), derive proper estimation methods for each of them and analyze their properties under some data problems. Finally, we present an interesting application.

2. Different Heterogeneity Formulations

The most widely used fixed effects model specifications have been proposed byBaltagi et al. [2003], Egger and Pfanffermayr [2003],Baldwin and Taglioni [2006], and Baier and Bergstrand[2007]. The straightforward direct generalization of the standard fixed effects panel data model (where the usual individuals are in fact the (ij) country pairs) takes into account bilateral interaction. The model specification is

y_ijt =β^′x_ijt+γ_ij+ε_ijt i = 1, . . . , N j = 1, . . . , N, t= 1, . . . , T

where the γij are the bilateral specific fixed effects. If the specification is used in a macro trade model, for example, with say 150 countries involved, this explicitly or implicitly, means the estimation of 150×150 = 22,500 parameters. This looks very much like a textbook over-specification case. Instead we propose, like in a standard panel data context, the use of the much more parsimonious random effects specification

yijt =β^′xijt+µij+εijt i= 1, . . . , N, j = 1, . . . , N, t = 1, . . . , T (1)

where E(µij) = 0, the random effects are pairwise uncorrelated, and E(µijµi^′j^′) =

σ²_µ i=i^′ and j =j^′ 0 otherwise

A natural extension of this model is to include time effects as well

yijt =β^′xijt+µij +λt +εijt i= 1, . . . , N j = 1, . . . , N, t= 1, . . . , T (2) where E(λt) = 0 and

E(λtλ^′_t) =

σ_λ² t=t^′ 0 otherwise

Another form of heterogeneity is to use individual-time-varying effects yijt =β^′xijt+αjt+εijt

The corresponding random effects specification now is

yijt =β^′xijt+ujt+εijt (3) where E(ujt) = 0, the random effects are pairwise uncorrelated, and

E(u_iju_j^′_t^′) =

σ_u² j =j^′ and t=t^′ 0 otherwise

Or alternatively we can also have the following random effects specification

yijt =β^′xijt+vit+εijt (4) where E(vit) = 0, the random effects are pairwise uncorrelated, and

E(vitvi^′t^′) =

σ²_v i=i^′ and t =t^′ 0 otherwise

The random effects specification containing both the above forms of heterogeneity now is

yijt =β^′xijt +vit+ujt+εijt (5) The model specification which encompasses all above effects is

yijt =β^′xijt+γij+αit+αjt+εijt

The corresponding random effects specification now is

yijt =β^′xijt+µij +vit+ujt+εijt (6) where E(µij) = 0, E(ujt) = 0, E(vit) = 0, all random effects are pairwise uncorrelated, and

E(µijµi^′j^′) =

σ_µ² i=i^′ and j =j^′ 0 otherwise

E(ujtuj^′t^′) =

σ_u² j =j^′ and t =t^′ 0 otherwise

E(vitvi^′t^′) =

σ_v² i=i^′ and t =t^′ 0 otherwise

In order to estimate efficiently these random effects models their corresponding covariance matrices need to be derived

3. Covariance Matrices of the Different Random Effects Specifications The standard way to estimate these models is with the Feasible GLS (FGLS) estimator. First, we need to derive the covariance matrix of each of the models introduced in Section 2, then the unknown variance components of these matrices need to be estimated.

For model (1) let us denote

u^⋆_ijt =µij +ǫijt (7)

So for all t observations

u^⋆_ij =µij⊗lT +ǫij

E

u^⋆_iju^⋆′_ij

=E[(µij⊗lT) (µij⊗l^′_T)] +E ǫijǫ^′_ij

=σ_µ²JT +σ_ǫ²IT

wherelT is the (T ×1) vector of ones,JT is the (T ×T) matrix of ones and IT is the (T ×T) identity matrix. In all the paper matrix J will denote the matrix of ones, with the size in the index, and I the identity matrix, also with the size in the index.

Now for individuali

u^⋆_i =µi⊗lT +ǫi

E

u^⋆_iu^⋆′_i

=E[(µi⊗lT) (µ^′_i⊗l^′_T)] +E[ǫiǫ^′_i]

=σ_µ²I_N ⊗J_T +σ²_ǫI_{N T}

And combining all these results we get for the covariance matrix of model (1) u^⋆ =µ⊗lT +ǫ

E u^⋆u^⋆′

=E[(µ⊗l_T) (µ^′⊗l^′_T)] +E[ǫǫ^′]

=σ²_µI_N² ⊗JT +σ_ǫ²I_N²_T = Ω Deriving likewise the covariance matrix for model (2)

u^⋆_ij =µij ⊗lT +λ+ǫij

E

u^⋆_iju^⋆′_ij

=E

(µij⊗lT) (µij⊗lT)^′

+E[λλ^′] +E ǫijǫ^′_ij

=σ_µ²J_T +σ_λ²I_T +σ_ǫ²I_T and

u^⋆_i =µi⊗lT +lN ⊗λ+ǫi

E u^⋆_iu^⋆′_i

=E

(µ_i⊗l_T) (µ_i⊗l_T)^′ +E

(l_N ⊗λ) (l_N ⊗λ)^′

+E[ǫ_iǫ^′_i]

=σ_µ²IN ⊗JT +σ_λ²JN ⊗IT +σ²_ǫIN T

so we obtain

u^⋆ =µ⊗lT +l_N² ⊗λ+ǫ E

u^⋆u^⋆′

=E

(µ⊗lT) (µ⊗lT)^′ +E

(l_N² ⊗λ) (l_N² ⊗λ)^′

+E[ǫǫ^′]

=σ_µ²I_N² ⊗JT +σ_λ²J_N² ⊗IT +σ²_ǫI_N²_T = Ω

Let us turn now to models (3) and (4) which can be dealt with in the same way as they are completely symmetric

u^⋆_ijt =ujt+ǫijt (8)

u^⋆_ij =uj +ǫij

E u^⋆_iju^⋆′_ij

=E uju^′_j

+E ǫijǫ^′_ij

=σ²_uIT +σ_ǫ²IT

u^⋆_i =u+ǫi

E u^⋆_iu^⋆′_i

=E[uu^′] +E[ǫ_iǫ^′_i] =σ_u²I_{N T} +σ²_ǫI_{N T} u^⋆ =lN ⊗u+ǫ

E u^⋆u^⋆′

=E[(lN ⊗u) (l_N^′ ⊗u^′)] +E[ǫǫ^′] =σ_u²JN ⊗IN T +σ_ǫ²IN²T = Ω

Using the same approach, the covariance matrix for model (5) is u^⋆_ijt =ujt+vit +ǫijt

u^⋆_ij =uj+vi+ǫij

E u^⋆_iju^⋆′_ij

=E uju^′_j

+E[viv_i^′] +E ǫijǫ^′_ij

=σ²_uIT +σ²_vIT +σ²_ǫIT

u^⋆_i =lN ⊗vi+u+ǫi

E u^⋆_iu^⋆′_i

=E[(l_N ⊗v_i) (l^′_N ⊗v^′_i)] +E[uu^′] +E[ǫ_iǫ^′_i] =

=σ²_vJN ⊗IT +σ_u²IN T +σ²_ǫIN T

and so E u^⋆u^⋆′

=σ²_v(IN ⊗JN ⊗IT) +σ²_u(JN ⊗IN T) +σ_ǫ²IN² = Ω And finally the covariance matrix of the all encompassing model (6) is

u^⋆_ijt =µij+ujt+vit+ǫijt (9)

u^⋆_ij =µij⊗lT +uj +vi+ǫij

E u^⋆_iju^⋆′_ij

=E[(µ_ij⊗l_T) (µ_ij⊗l^′_T)] +E u_ju^′_j

+E[v_iv_i^′] +E ǫ_ijǫ^′_ij

=σ_µ²JT +σ_u²IT +σ_v²IT +σ_ǫ²IT

u^⋆_i =µi⊗lT +lN ⊗vi+u+ǫi

E u^⋆_iu^⋆′_i

=E[(µi⊗lT) (µ^′_i⊗l^′_T)] +E[(lN ⊗vi) (l^′_N ⊗v^′_i)] +E[uu^′] +E[ǫiǫ^′_i] =

=σ_µ²IN ⊗JT +σ²_uIN T +σ_v²JN ⊗IT +σ²_ǫIN T

and so E u^⋆u^⋆′

=σ_µ²(I_N² ⊗JT) +σ_u²(JN ⊗IN T) +σ_v²(IN ⊗JN ⊗IT) +σ_ǫ²I_N²_T = Ω

4. Estimation of the Variance Components and the Feasible GLS Estimator

Turning now to the estimation of the variance components of the different models, let us start with model (1)

Eh u^⋆2_ijti

=Eh

(µij +ǫijt)²i

=E µ²_ij

+E ǫ²_ijt

=σ²_µ+σ_ǫ² (10) and let us introduce the appropriate Within transformation

u^⋆_ijt,within=u^⋆_ijt−u¯^⋆ij =ǫijt−¯ǫij (11) where ¯ǫ_ij = 1/T P

tǫ_ijt and ¯u^∗_ij = 1/T P

tu^∗_ijt, so we get Eh

u^⋆_ijt−u¯^⋆ij

²i

=Eh

(ǫijt−¯ǫij)²i

=E



ǫ²_ijt−2ǫijt

1 T

T

X

t=1

ǫijt+ 1 T

T

X

t=1

ǫijt

!²



=E ǫ²_ijt

−2E

"

ǫijt

1 T

T

X

t=1

ǫijt

# +E



 1 T

T

X

t=1

ǫijt

!²



=σ_ǫ²− 2

Tσ_ǫ²+ 1

Tσ²_ǫ =σ_ǫ²− 1

Tσ_ǫ² =σ_ǫ²T −1 T

Let ˆu^∗ be the OLS residual of model (1) and ˆu^∗_within the Within transformation of this residual. Then we can estimate the variance components as

ˆ

σ²_ǫ = T

T −1uˆ^⋆_within^′ uˆ^⋆_within ˆ

σ²_µ= 1 N²T

N

X

i=1 N

X

j=1 T

X

t=1

ˆ

u^⋆_ijt² −σˆ_ǫ²

These estimators naturally should be adjusted to the actual degrees of freedom.

Continuing with model (2) Eh

u^⋆2_ijti

=Eh

(µ_ij +λ_t+ǫ_ijt)²i

=E µ²_ij

+E λ²_t

+E ǫ²_ijt

=σ_µ² +σ²_λ+σ_ǫ² E



 1 T

T

X

t=1

u^⋆ijt

!²

=E



 1 T

T

X

t=1

µij +λt+ǫijt

!²



=E µ²_ij

+ 1 T²E

" _T X

t=1

λ²_t

# + 1

T²E

" _T X

t=1

ǫ²_ijt

#

=σ_µ² + 1

Tσ_λ²+ 1 T σ_ǫ²

and Eh

u^⋆_ijt−u¯^⋆_ij−u¯^⋆_t+ ¯u^⋆2i

=Eh

(ǫ_ijt−¯ǫ_ij −¯ǫ_t+ ¯ǫ)²i

=E ǫ²_ijt

+E



 1 T

T

X

t=1

ǫijt

!²

+

+E









 1 N²

N

X

i=1 N

X

j=1

ǫijt





2



+E









 1 N²T

N

X

i=1 N

X

j=1 T

X

t=1

ǫijt





2



−

−2E

"

ǫijt

1 T

T

X

t=1

ǫijt

#

−2E



ǫijt

1 N²

N

X

i=1 N

X

j=1

ǫijt



+

+ 2E



ǫijt

1 N²T

N

X

i=1 N

X

j=1 T

X

t=1

ǫijt



+ 2E



 1 T

T

X

t=1

ǫijt· 1 N²

N

X

i=1 N

X

j=1

ǫijt



−

−2E



 1 T

T

X

t=1

ǫijt· 1 N²T

N

X

i=1 N

X

j=1 T

X

t=1

ǫijt



−2E



 1 N²

N

X

i=1 N

X

j=1

ǫijt· 1 N²T

N

X

i=1 N

X

j=1 T

X

t=1

ǫijt





=σ_ǫ²+ 1

Tσ²_ǫ + 1

N²σ_ǫ²+ 1

N²Tσ²_ǫ − 2

T σ_ǫ²− 2 N²σ_ǫ²+

+ 2

N²Tσ_ǫ²+ 2

N²Tσ²_ǫ − 2

N²Tσ_ǫ²− 2

N²Tσ²_ǫ =

=σ_ǫ²(N −1)(N + 1)(T −1) N²T

This leads to the estimation of the variance components ˆ

σ_ǫ² = N²T

(N −1)(N + 1)(T −1)uˆ^⋆′_withinuˆ^⋆within

ˆ

σ_µ² = 1 N²T(T −1)





N

X

i=1 N

X

j=1





T

X

t=1

ˆ u^⋆_ijt

!²

−

T

X

t=1

ˆ u^⋆_ijt²









ˆ

σ_λ² = 1 N²T

N

X

i=1 N

X

j=1 T

X

t=1

ˆ u^⋆_ijt²

−σˆ²_µ−σˆ²_ǫ

Turning now to models (3) and (4) Eh

u^⋆2_ijti

=Eh

(ujt+ǫijt)²i

=E u²_jt

+E ǫ²_ijt

=σ²_u+σ_ǫ² (12)

and the appropriate Within transformation now is

u^⋆_ijt,within =u^⋆_ijt−u¯^⋆jt =ǫijt−¯ǫjt (13) where ¯u^∗_jt = 1/NP

iu^∗_ijt and ¯ǫjt = 1/N P

iǫijt and Eh

u^⋆_ijt−u¯^⋆_jt²i

=Eh

(ǫ_ijt−¯ǫ_jt)²i

=E



ǫ²_ijt−2ǫijt

1 N

N

X

i=1

ǫijt+ 1 N

N

X

i=1

ǫijt

!²



=E ǫ²_ijt

−2E

"

ǫijt

1 N

N

X

i=1

ǫijt

# +E



 1 N

N

X

i=1

ǫijt

!²



=σ²_ǫ − 2

Nσ_ǫ²+ 1

Nσ_ǫ² =σ_ǫ²− 1

Nσ_ǫ² =σ²_ǫN −1 N And the estimators for the variance components are

ˆ

σ²_ǫ = N

N −1uˆ^⋆_within^′ uˆ^⋆_within ˆ

σ²_µ= 1 N²T

N

X

i=1 N

X

j=1 T

X

t=1

ˆ

u^⋆_ijt² −σˆ_ǫ² Now for model (5) the Within transformation is

u^∗_ijt,within = (u^∗_ijt−1/NX

i

u^∗_ijt−1/NX

j

u^∗_ijt+ 1/N²X

i

X

j

u^∗_ijt) (14) so we get

Eh

u^⋆_ijt−u¯^⋆jt−u¯^⋆it+ ¯u^⋆t

²i

=Eh

(ǫijt−¯ǫjt−¯ǫit + ¯ǫt)²i

=E ǫ²_ijt

+E



 1 N²

N

X

i=1

ǫijt

!²

+E





 1 N²





N

X

j=1

ǫijt





2



+E





 1 N⁴





N

X

i=1 N

X

j=1

ǫijt





2



−

−2E

"

ǫijt

1 N

N

X

i=1

ǫijt

#

−2E



ǫijt

1 N

N

X

j=1

ǫijt



+ 2E



ǫijt

1 N²

N

X

i=1 N

X

j=1

ǫijt



+

+ 2E



 1 N²

N

X

i=1

ǫijt N

X

j=1

ǫijt



−2E



 1 N³

N

X

i=1

ǫijt N

X

i=1 N

X

j=1

ǫijt



−2E



 1 N³

N

X

j=1

ǫijt N

X

i=1 N

X

j=1

ǫijt



=

=σ_ǫ²+ 1

Nσ_ǫ²+ 1

Nσ_ǫ²+ 1

N²σ_ǫ²− 2

Nσ²_ǫ − 2

Nσ_ǫ²+ 2

N²σ²_ǫ + 2

N²σ_ǫ²− 2

N²σ_ǫ²− 2 N²σ_ǫ² =

=σ_ǫ²

1− 2 N + 1

N²

=σ_ǫ²

N²−2N + 1 N²

=σ_ǫ²(N −1)² N²

(15)

And, also,

Eh u^⋆2_ijti

=Eh

(u_jt+v_it +ǫ_ijt)²i

=σ²_u+σ_v²+σ_ǫ² E



 1 N

N

X

i=1

u^⋆_ijt

!²

=E



 1 N

N

X

i=1

(ujt+vit +ǫijt)

!²



=E u²_jt

+ 1 N²E

" _N X

i=1

v_it²

# + 1

N²E

" _N X

i=1

ǫ²_ijt

#

=σ_u²+ 1

Nσ_v²+ 1 Nσ²_ǫ

(16)

The estimators of the variance components therefore are ˆ

σ_ǫ² = N²

(N −1)²uˆ^⋆_within^′ uˆ^⋆_within ˆ

σ_u² = 1 N²T(N −1)





N

X

j=1 T

X

t=1





N

X

i=1

ˆ uijt⋆

!²

−

N

X

i=1

ˆ uijt⋆2









ˆ

σ_v² = 1 N²T

N

X

i=1 N

X

j=1 T

X

t=1

ˆ

u^⋆2_ijt−σˆ²_ǫ −σˆ²_u

Finally, to derive the estimators of the variance components for model (6), we need first the appropriate Within transformation

u^∗_ijt,within= (u^∗_ijt−1/T X

t

u^∗_ijt−1/NX

i

u^∗_ijt−1/NX

j

u^∗_ijt+ 1/N²X

i

X

j

u^∗_ijt

+ 1/(N T)X

i

X

t

u^∗_ijt+ 1/(N T)X

j

X

t

u^∗_ijt−1/(N²T)X

i

X

j

X

t

u^∗_ijt) Carrying out the derivation as earlier, we get to the following estimators

ˆ

σ_ǫ² = N²T

N(N −1)(T −1) + 1uˆ^⋆_within^′ uˆ^⋆_within ˆ

σ_v² = 1 N²T(N −1)







N

X

i=1 T

X

t=1











N

X

j=1

ˆ u^⋆





2

−

N

X

j=1

ˆ u^⋆2













ˆ

σ_u² = 1 N²T(N −1)





N

X

j=1 T

X

t=1





N

X

i=1

ˆ u^⋆

!²

−

N

X

i=1

ˆ u^⋆2









ˆ

σ²_µ = 1 N²T

N

X

i=1 N

X

j=1 T

X

t=1

ˆ

u^⋆2_ijt−σˆ²_ǫ −σˆ_v²−σˆ_u²

Now we have all the tools to properly use the FGLS estimators.

5. Unbalanced Data

Like in the case of the usual panel data models, just more frequently, one may be faced with a situation when the data at hand is unbalance. In our framework of analysis this means that for all models (1)-(6) in general t = 1, . . . , Tij, P

i

P

jTij = T and Tij often is not equal to Ti^′j^′. For this unbalanced data case, as we did when the data was balanced, we need to derive the covariance matrices of the models and the appropriate estimators for the variance components.

For model (1), using decomposition (7) we get u^⋆_ij =µij ⊗lTij +ǫij

E

u^⋆_iju^⋆′_ij

=Eh

µ_ij ⊗l_Tij

µ_ij⊗l_Tij^′i +E

ǫ_ijǫ^′_ij

=

=σ_µ²JTij +σ_ǫ²ITij

and u^⋆_i = ˜µi+ǫi

E

u^⋆_iu^⋆′_i

=E

˜ µiµ˜i

′

+E[ǫiǫ^′_i]

=σ_µ²A+σ²_ǫIPN j=1Tij

where µ˜i =











































 µi1

... µi1

µi2

... µi2

... µiN

... µiN













































, A=











I_Ti1 0 . . . 0 0 I_Ti2 . . . 0 ... ... . .. ... 0 0 . . . ITiN











of size

N

X

j=1

Tij ×

N

X

j=1

Tij

and finally for the complete model

u^⋆ = ˜µ+ǫ E

u^⋆u^⋆′

=E[˜µ˜µ^′] +E[ǫǫ^′]

=σ²_µB+σ_ǫ²IT

where µ˜=

































































 µ11

... µ11

µ12

... µ¹²

... µij

... µij

... µ_{N N}

... µN N



































































, B =











JT₁₁ 0 . . . 0 0 JT₁₂ . . . 0 ... ... . .. ... 0 0 . . . JTN N











of size (T ×T)

Continuing with model (2)

u^⋆_ij =µij ⊗lTij +λ+ǫij

E

u^⋆_iju^⋆′_ij

=Eh

µij ⊗lTij

µij⊗lTij

^′i

+E[λλ^′] +E ǫijǫ^′_ij

=σ_µ²JTij +σ_λ²ITij +σ_ǫ²ITij

u^⋆_i = ˜µ_i+ ˜λ_i+ǫ_i where

λ˜^′_i = (λ1, λ2, . . . , λ_Ti1, . . . , λ1, λ2, . . . , λ_TiN) E

u^⋆_iu^⋆′_i

=E

˜ µiµ˜i

′ +Eh

λ˜iλ˜i

′i

+E[ǫiǫ^′_i]

=σ_µ²A+σ²_λDi+σ_ǫ²IPN j=1Tij

u^⋆ = ˜µ+ ˜λ+ǫ E

u^⋆u^⋆′

=E[˜µ˜µ^′] +Eh

˜λλ˜^′i

+E[ǫǫ^′]

=σ_µ²B+σ²_λE+σ_ǫ²I_T

with

E(E¹¹, E¹², . . . , E¹N, . . . , EN1, EN2, . . . , EN N)

Eij =











MT₁₁×Tij

MT₁₂×Tij

... MTN N×Tij











and Di =











IT_i1 MT_i1×T_i2 . . . MT_i1×TiN

MT_i2×T_i1 IT_i1 . . . MT_i2×TiN

... ... . .. ... MTiN×T_i1 MTiN×T_i2 . . . ITiN











where

MTij×Tlj =











1 0 . . . 0 0 . . . 0 0 1 . . . 0 0 . . . 0

... ... . .. ... ... ...

0 0 . . . 1 0 . . . 0











if Tlj > Tij

and

MTij×Tlj =

























1 0 . . . 0 0 1 . . . 0 ... ... . .. ...

0 0 . . . 1 0 0 . . . 0 ... ... . .. ...

0 0 . . . 0

























if Tlj < Tij

Doing the same exercise for model (3) using decomposition (8) we end up with u^⋆_ij =uj +ǫij

E u^⋆_iju^⋆′_ij

=E uju^′_j

+E ǫijǫ^′_ij

=σ_u²ITij +σ_ǫ²ITij

u^⋆_i =u+ǫi

E u^⋆_iu^⋆′_i

=E[uu^′] +E[ǫiǫ^′_i] =σ²_uIPN

j=1Tij +σ_ǫ²IPN j=1Tij

u^⋆ = ˜u+ǫ and so for the complete model we get

E u^⋆u^⋆′

=E[˜uu˜^′] +E[ǫǫ^′] =σ²_uC+σ_ǫ²IT

where

˜

u^′ = (u11, . . . , u1T₁₁, . . . , uN1, . . . , uN T_1N, . . . , u11, . . . , u1T_N1, . . . , uN1, . . . , uN TN N)

C = (C¹, C², C³)

C1 =



























































IT₁₁ 0 . . . 0

0 IT₁₂ . . . 0

... ... . .. ...

0 0 . . . IT_1N

MT₂₁×T₁₁ 0 . . . 0 0 MT₂₂×T₁₂ . . . 0 ... ... . .. ... 0 0 . . . M_T2N×T1N

... ... . .. ... MTN1×T11 0 . . . 0 0 MTN2×T₁₂ . . . 0 ... ... . .. ... 0 0 . . . MTN N×T_1N



























































C² =



























































MT₁₁×T₂₁ 0 . . . 0 . . . 0 MT₁₂×T₂₂ . . . 0 . . . ... ... . .. ... . . . 0 0 . . . M_T1N×T2N . . .

IT₂₁ 0 . . . 0 . . .

0 IT₂₂ . . . 0 . . .

... ... . .. ... . . .

0 0 . . . IT_2N . . .

... ... . .. ... . . . MT_N1×T₂₁ 0 . . . 0 . . . 0 MT_N2×T₂₂ . . . 0 . . . ... ... . .. ... . . . 0 0 . . . M_{TN N}^×_T1N . . .



























































C3 =



























































MT₁₁×TN1 0 . . . 0 0 MT₁₂×T_N2 . . . 0 ... ... . .. ... 0 0 . . . MT_1N×TN N

MT₂₁×T_N1 0 . . . 0 0 MT₂₂×T_N2 . . . 0 ... ... . .. ... 0 0 . . . MT_2N×TN N

... ... . .. ...

ITN1 0 . . . 0

0 ITN2 . . . 0

... ... . .. ...

0 0 . . . ITN N



























































Let us now turn to model (4). Following the same steps as above, we get for the covariance matrix (σ_v²D+σ²_ǫIT) where

D=





D¹ 0 . . . 0 0 D2 . . . 0 0 0 . . . DN





Models (5) and (6) can be dealt with together using decomposition (9) u^⋆_ij =µij ⊗lT +uj+vi+ǫij

E u^⋆_iju^⋆′_ij

=Eh

µij ⊗lTij

µij ⊗l_T^′_iji +E

uju^′_j

+E[viv^′_i] +E ǫijǫ^′_ij

=σ_µ²JTij +σ_u²ITij +σ_v²ITij +σ_ǫ²ITij

u^⋆_i = ˜µ_i+ ˜v_i+u+ǫ_i E u^⋆_iu^⋆′_i

=E

˜ µiµ˜i

′ +E

˜ viv˜i

′

+E[uu^′] +E[ǫiǫ^′_i]

=σ_µ²A+σ²_uIPN

j=1Tij +σ²_vDi+σ_ǫ²IPN j=1Tij

u^⋆ = ˜µ+ ˜v+ ˜u+ǫ where v˜i

′ = (vi1, vi2, . . . , viT_i1, vi1, vi2, . . . , viT_i2, . . . , vi1, vi2, . . . , viTiN)

˜

v^′ = ( ˜v1,v˜2, . . . ,v˜_N,) E u^⋆u^⋆′

=E[˜µ˜µ^′] +E[˜v˜v^′] +E[˜uu˜^′] +E[ǫǫ^′] =

=σ_µ²B+σ_u²C+σ_v²D+σ²_ǫIT

For model (5) the appropriate covariance matrix is the same with B= 0.

Now that we derived the covariance matrices for unbalanced data it is time to turn to the estimation of the variance components. Using (10) and (11)

E



 1 N²

N

X

i=1 N

X

j=1

u^⋆_ijt−u¯^⋆ij

²



= 1 N²

N

X

i=1 N

X

j=1

Eh

(ǫijt−¯ǫij)²i

= 1 N²

N

X

i=1 N

X

j=1

E





ǫ²_ijt−2ǫijt

1 Tij

Tij

X

t=1

ǫijt+



 1 Tij

Tij

X

t=1

ǫijt





2





= 1 N²

N

X

i=1 N

X

j=1





E ǫ²_ijt

−2E



ǫijt

1 Tij

Tij

X

t=1

ǫijt



+E









 1 Tij

Tij

X

t=1

ǫijt





2











= 1 N²

N

X

i=1 N

X

j=1

σ²_ǫ − 2 Tij

σ_ǫ²+ 1 Tij

σ_ǫ²

=σ_ǫ² 1 N²

N

X

i=1 N

X

j=1

Tij −1 Tij

so for the variance components we get the following estimators ˆ

σ_ǫ² = N² PN

i=1

PN j=1

Tij−1 Tij

ˆ

u^⋆_within^′ uˆ^⋆_within

ˆ σ_µ² = 1

T

N

X

i=1 N

X

j=1 Tij

X

t=1

uˆ^⋆2_ijt−σˆ_ǫ²

For model (3) (and similarly for model (4)), using (12) and (13) and using the same derivations as there we get

ˆ

σ²_ǫ = N

N −1uˆ^⋆_within^′ uˆ^⋆_within ˆ

σ_u² = 1 T

N

X

i=1 N

X

j=1 Tij

X

t=1

ˆ

u^⋆2_ijt−σˆ_ǫ²

Turning now to model (5), as (14) and (15) are the same in the unbalanced case we get

ˆ

σ_ǫ² = N²

(N −1)²uˆ^⋆_within^′ uˆ^⋆_within ˆ

σ_u² = 1 N −1





N

X

i=1

1 PN

j=1Tij N

X

j=1 Tij

X

t=1

1 N

N

X

i=1

uˆ^⋆_ijt

!²

− 1 T

N

X

i=1 N

X

j=1 Tij

X

t=1

ˆ u^⋆_ijt²





σˆ_v² = 1 T

N

X

i=1 N

X

j=1 Tij

X

t=1

ˆ

u^⋆_ijt² −σˆ_ǫ²−σˆ_u²