Improved Average Estimation in Seemingly Unrelated Regressions

Abstract

In this paper, we propose an efficient weighted average estimator in Seemingly Unrelated Regressions. This average estimator shrinks a generalized least squares (GLS) estimator towards a restricted GLS estimator, where the restrictions represent possible parameter homogeneity specifications. The shrinkage weight is inversely proportional to a weighted quadratic loss function. The approximate bias and second moment matrix of the average estimator using the large-sample approximations are provided. We give the conditions under which the average estimator dominates the GLS estimator on the basis of their mean squared errors. We illustrate our estimator by applying it to a cost system for United States (U.S.) Commercial banks, over the period from 2000 to 2018. Our results indicate that on average most of the banks have been operating under increasing returns to scale. We find that over the recent years, scale economies are a plausible reason for the growth in average size of banks and the tendency toward increasing scale is likely to continue

1. Introduction

Seemingly unrelated regressions (SUR) was introduced by (Zellner 1962) and is one of the econometric developments that has been widely used in applied work. The relative ease of estimation, applying a large class of modeling and testing problems, and the availability of data representing a sample of cross section units observed over several time periods are related to the popularity of this model. (Zellner 1962) proposed a generalized least squares (GLS) estimator for estimating the coefficients of a set of SUR and established that it yields, at least asymptotically, to more efficient estimators than those obtained by single-equation least squares. See the surveys by (Srivastava and Dwivedi 1979; Fiebig 2001) and the book by (Srivastava and Giles 1987) for a concise coverage of the literature in this area.

Shrinkage estimations in SUR was first introduced by (Zellner and Vandaele 1975), which extends the results of (James and Stein 1961) and (Sclove 1968) to multivariate regression equations and presents a technique of constructing an estimator whose risk is smaller than the risk of the GLS estimator. However, the resulting estimator depends on some unknown matrices and is not practical. (Srivastava 1973) investigates the properties of the estimator when consistent estimators are substituted for these unknown matrices. (Maddala 1991) reviewed the shrinkage estimators and showed that these estimators appear to perform better than both pooled and single-equation least squares estimators (see also Maddala and Hu 1996; Maddala et al. 1997; Choi and Li 2000). Maddala et al. (2001) show the superior properties of shrinkage estimators among single-equation estimators and various averaging estimators in a heterogeneous panel data model under error homoscedasticity framework. In univariate equation models, recently, (Hansen 2016) introduces shrinkage for general parametric models by shrinking maximum likelihood estimators (MLE) toward a restricted MLE. Hansen (2016) shows the dominance of the shrinkage estimator over the MLE in terms of having lower asymptotic risk when the shrinkage dimension exceeds two, using a local to zero asymptotic framework. Wang et al. (2019) propose a Mallow pooling averaging estimator for heterogeneous panel data models and conclude that the pooling averaging estimator is preferred when the panel is heterogenous and the signal-to-noise ratio is moderate or large. For more on averaging estimators see (Ullah and Wang 2013).

In the analysis of SUR, a question that often practitioners are faced is whether to assume parameter homogeneity or parameter heterogeneity. On one hand, the parameter heterogeneity assumption results in consistent estimators and violation of this assumption causes misleading estimates (see, for example, Robertson and Symons 1992; Pesaran and Smith 1995; Su and Chen 2013; Durlauf et al. 2001; Browning and Carro 2007). On the other hand, the parameter homogeneity assumption causes higher efficiency but could be at the cost of estimation bias and inconsistency of estimators, which is supported by an increasing number of studies due to a better forecast performance of the estimators under this assumption (see, for example, Maddala 1991; Maddala and Hu 1996; Baltagi and Griffin 1984; Baltagi et al. 2000; Hoogstrate et al. 2000). This question and the results of the mentioned research show the typical bias-variance trade-off that needs to be considered in choosing the restrictions. While efficiency is important, robustness is also critical, since researchers prefer as few ad hoc restrictions as possible. In the present scenario the efficient estimator depends on more stringent condition of homogeneity and therefore is less robust to the heterogeneity restriction. Therefore, this efficiency-robustness trade-off (bias-variance trade-off) calls for thorough examination. A natural approach is to consider a pre-test estimator, but it is proven unable to solve the efficiency-robustness issue (Leeb and Pötsche 2005).

A more useful approach considered here is an averaging estimator (Stein-type shrinkage), that is a weighted average of the robust and efficient (the unrestricted GLS and the restricted GLS) estimators. The weight is inversely related to a weighted quadratic loss function, which measures the weighted distance between the unrestricted and the restricted GLS estimators. The first and second moments of our proposed average estimator are derived using (Nagar 1959) large-sample approximations. Furthermore, we show the dominance properties in terms of mean squared error (MSE) of the estimators, which ensures that the proposed estimator is robust against arbitrary deviations from the restrictions. This is an advantage of our method relative to the “local asymptotic” argument that some previous studies rely on (see, for example, Hansen 2016). Our dominance property ensures that the proposed averaging estimator is robust against arbitrary deviations from the restrictions, while previous estimators in the literature consider mainly very small violations of the restrictions. Further discussion in this area has been generally theoretical, but we present here, an important empirical question and show the advantages of the averaging estimator over the previous estimators considered. Essentially, we apply our estimator to estimate cost efficiency of United States (U.S.) Commercial banks using a cost system method over the period from 2000 to 2018. Since bank size is an important factor of production environment, following the literature, we use it to partition bank technologies. However, these partitions are user-specified, and the estimates can be misleading because of false parameter heterogeneity assumptions. Therefore, we use the average estimator introduced in the paper to estimate the cost efficiency, as it optimally balances the trade-off between bias and variance efficiency of the restricted and the unrestricted GLS estimators. We find that on average majority of banks have been operating under increasing returns to scale over the sample period. We also find more signs of cost efficiency for Large Banks (banks with asset size more than $500 Million dollars) and Small Banks (banks with asset size less than $100 Million dollars) relative to Medium banks (banks with asset size between $500 and $100 Million dollars). This finding is important for gauging costs and benefits of any policy intervention to control the size of banks.

The paper is organized as follows. Section 2 describes the model and the assumptions. In Section 3, we introduce the estimators. We give the bias, MSE matrix and the risk of the average estimator using the large-sample approximations in Section 4. Section 5 reports some Monte Carlo simulations to evaluate the accuracy of the approximations. Results from our empirical example are presented in Section 6. Finally, Section 7 contains some concluding remarks, and proofs are given in Appendix A.

2. The Model and Notation

Consider the following m seemingly unrelated linear regressions

$${\mathit{y}}_i={\mathit{X}}_i{\mathit{\beta }}_i+{\mathit{u}}_i, i=1,2,\dots m,$$

(1)

where ${\mathit{y}}_i={(y_{i1},y_{i2},\dots ,y_{iT})}^{\prime }$ is a $T\times 1$ vector of observations on the dependent variable $y_{it},$ with T being the number of observations, ${\mathit{X}}_i$ is a $T\times k$ matrix of observations on the k vector of regressors including the intercept (that is, $x_{it,1}=1$)1, ${\mathit{\beta }}_i$ is a $k\times 1$ vector of unknown coefficients and ${\mathit{u}}_i={(u_{i1},u_{i2},\dots ,u_{iT})}^{\prime }$ is a $T\times 1$ vector of disturbances, for $i=1,2,\dots ,m.$ It is convenient to stack the m equations above in the following form:

$$\left[\begin{array}{c}{\mathit{y}}_1\\ {\mathit{y}}_2\\ ⋮\\ {\mathit{y}}_m\end{array}\right]=\left[\begin{array}{cccc}{\mathit{X}}_1& \mathbf{0}& \dots & \mathbf{0}\\ \mathbf{0}& {\mathit{X}}_2& \\ ⋮& \ddots & \\ \mathbf{0}& \dots & \mathbf{0}& {\mathit{X}}_m\end{array}\right]\left[\begin{array}{c}{\mathit{\beta }}_1\\ {\mathit{\beta }}_2\\ ⋮\\ {\mathit{\beta }}_m\end{array}\right]+\left[\begin{array}{c}{\mathit{u}}_1\\ {\mathit{u}}_2\\ ⋮\\ {\mathit{u}}_m\end{array}\right],$$

(2)

or compactly as

$$\underset{mT\times 1}{\mathit{y}}=\underset{mT\times mk}{\mathit{X}}\underset{mk\times 1}{\mathit{\beta }}+\underset{mT\times 1}{\mathit{u}}.$$

(3)

We assume,

Assumption 1.

The $mT\times 1$ vector of disturbances, $\mathit{u},$ has a zero conditional mean

$$\mathbb{E}(\mathit{u}|{\mathit{X}}_1,{\mathit{X}}_2,\dots ,{\mathit{X}}_m)=\mathbf{0}.$$

Assumption 2.

The disturbances are uncorrelated across observations but correlated across equations,

$$\mathbb{E}({\mathit{u}}_i{\mathit{u}}_j^{\prime }|{\mathit{X}}_1,{\mathit{X}}_2,\dots ,{\mathit{X}}_m)={\sigma }_{ij}{\mathit{I}}_T,$$

or

$$\begin{array}{l}\mathbb{E}(\mathit{u}{\mathit{u}}^{\prime }|{\mathit{X}}_1,{\mathit{X}}_2,\dots ,{\mathit{X}}_m)=\mathbf{\Omega }=\left[\begin{array}{cccc}{\sigma }_{11}{\mathit{I}}_T& {\sigma }_{12}{\mathit{I}}_T& \dots & {\sigma }_{1m}{\mathit{I}}_T\\ {\sigma }_{21}{\mathit{I}}_T& {\sigma }_{22}{\mathit{I}}_T& \dots & {\sigma }_{2m}{\mathit{I}}_T\\ ⋮\\ {\sigma }_{m1}{\mathit{I}}_T& {\sigma }_{m2}{\mathit{I}}_T& \dots & {\sigma }_{mm}{\mathit{I}}_T\end{array}\right]=\mathbf{\Sigma }\otimes {\mathit{I}}_T,\end{array}$$

where ${\mathit{I}}_T$ is the $T\times T$ identity matrix,

$$\underset{m\times m}{\mathbf{\Sigma }}=\left[\begin{array}{cccc}{\sigma }_{11}& {\sigma }_{12}& \dots & {\sigma }_{1m}\\ {\sigma }_{21}& {\sigma }_{22}& \dots & {\sigma }_{2m}\\ & & ⋮\\ {\sigma }_{m1}& {\sigma }_{m2}& \dots & {\sigma }_{mm}\end{array}\right],$$

and we assume Ω is positive definite.

Assumption 3.

The disturbances are normally distributed with mean zero and variance-covariance matrix $\mathbf{\Omega }.$

We define some notations below, which will be used in the following sections. So let

$$\underset{mT\times mT}{\mathit{Q}}={\mathbf{\Omega }}^{-1}({\mathit{I}}_{mT}-\mathbf{\Psi }),$$

(4)

where

$$\underset{mT\times mT}{\mathbf{\Psi }}=\mathit{X}{({\mathit{X}}^{\prime }{\mathbf{\Omega }}^{-1}\mathit{X})}^{-1}{\mathit{X}}^{\prime }{\mathbf{\Omega }}^{-1}.$$

(5)

If we partition $\mathit{Q}$ in the sub-matrices of $T\times T$ as below

$$\mathit{Q}=\left[\begin{array}{cccc}{\mathit{Q}}_{11}& {\mathit{Q}}_{12}& \dots & {\mathit{Q}}_{1m}\\ {\mathit{Q}}_{21}& {\mathit{Q}}_{22}& \dots & {\mathit{Q}}_{2m}\\ ⋮& ⋮& & ⋮\\ {\mathit{Q}}_{m1}& {\mathit{Q}}_{m2}& \dots & {\mathit{Q}}_{mm}\end{array}\right],$$

(6)

then we define

$$\mathbf{\Pi }=\left[\begin{array}{cccc}{\mathit{Q}}_{11}^{\prime }& {\mathit{Q}}_{12}^{\prime }& \dots & {\mathit{Q}}_{1m}^{\prime }\\ {\mathit{Q}}_{21}^{\prime }& {\mathit{Q}}_{22}^{\prime }& \dots & {\mathit{Q}}_{2m}^{\prime }\\ ⋮& ⋮& & ⋮\\ {\mathit{Q}}_{m1}^{\prime }& {\mathit{Q}}_{m2}^{\prime }& \dots & {\mathit{Q}}_{mm}^{\prime }\end{array}\right].$$

(7)

Also, we define

$$\underset{T\times T}{\mathbf{\Phi }}=\sum _{i=1}^m{\mathbf{\Psi }}_{ii},$$

(8)

where ${\mathbf{\Psi }}_{ii}$ is the ith diagonal $T\times T$ sub-matrix of $\mathbf{\Psi }$ which is partitioned as below

$$\mathbf{\Psi }=\left[\begin{array}{cccc}{\mathbf{\Psi }}_{11}& {\mathbf{\Psi }}_{12}& \dots & {\mathbf{\Psi }}_{1m}\\ {\mathbf{\Psi }}_{21}& {\mathbf{\Psi }}_{22}& \dots & {\mathbf{\Psi }}_{2m}\\ ⋮& ⋮& & ⋮\\ {\mathbf{\Psi }}_{m1}& {\mathbf{\Psi }}_{m2}& \dots & {\mathbf{\Psi }}_{mm}\end{array}\right].$$

(9)

3. Estimators

Our goal is to estimate the vector of slope parameters, $\mathit{\beta },$ in Equation (3). We consider three estimators of the slope parameters. The first estimator is the (Zellner 1962) GLS estimator (the unrestricted GLS estimator), which is the standard estimator in SUR. The second estimator is a restricted GLS estimator that ignores the slope parameters heterogeneity and estimates a pooled model. The third estimator, called the average estimator, is a weighted average of the restricted and the unrestricted GLS estimators where the weight is proportional to a weighted quadratic loss function.

3.1. Unrestricted Estimator

The typical estimator of the slope parameters in SUR is a feasible GLS estimator defined as

$$\begin{array}{ll}\hfill \widehat{\mathit{\beta }}& ={({\mathit{X}}^{\prime }{\widehat{\mathbf{\Omega }}}^{-1}\mathit{X})}^{-1}{\mathit{X}}^{\prime }{\widehat{\mathbf{\Omega }}}^{-1}\mathit{y}=\mathit{\beta }+{({\mathit{X}}^{\prime }{\widehat{\mathbf{\Omega }}}^{-1}\mathit{X})}^{-1}{\mathit{X}}^{\prime }{\widehat{\mathbf{\Omega }}}^{-1}\mathit{u},\hfill \end{array}$$

(10)

where $\widehat{\mathbf{\Omega }}$ is an estimator of $\mathbf{\Omega },$ which can be calculated as

$$\widehat{\mathbf{\Omega }}=\widehat{\mathbf{\Sigma }}\otimes {\mathit{I}}_T,$$

(11)

where $\widehat{\mathbf{\Sigma }}$ is an estimator of $\mathbf{\Sigma },$ such that its $(i,j)$th element, $s_{ij},$ estimates ${\sigma }_{ij}$, using a single-equation estimator of ${\mathit{\beta }}_i$, defined as ${\breve{\mathit{\beta }}}_i={({\mathit{X}}_i^{\prime }{\mathit{X}}_i)}^{-1}{\mathit{X}}_i^{\prime }{\mathit{y}}_i,\mathrm{for}i=1,2,\dots ,m$. Hence, $s_{ij}$ is equal to

$$s_{ij}={({\mathit{y}}_i-{\mathit{X}}_i{\breve{\mathit{\beta }}}_i)}^{\prime }({\mathit{y}}_j-{\mathit{X}}_j{\breve{\mathit{\beta }}}_j)/T={\mathit{u}}_i^{\prime }{\mathit{M}}_i{\mathit{M}}_j{\mathit{u}}_j/T,$$

(12)

where ${\mathit{M}}_i={\mathit{I}}_T-{\mathit{X}}_i{({\mathit{X}}_i^{\prime }{\mathit{X}}_i)}^{-1}{\mathit{X}}_i^{\prime }$ is an idempotent projection matrix.

3.2. Restricted Estimator

The restricted estimator is defined under the parameter homogeneity assumption across equations, which can be written as

$${\mathit{\beta }}_1={\mathit{\beta }}_2=\cdots ={\mathit{\beta }}_m=\overline{\mathit{\beta }},$$

where $\overline{\mathit{\beta }}$ is a weighted average of the slope parameters, ${\mathit{\beta }}_i$’s, defined as

$$\overline{\mathit{\beta }}={({\mathit{J}}^{\prime }{\mathit{X}}^{\prime }{\mathbf{\Omega }}^{-1}\mathit{X}\mathit{J})}^{-1}{\mathit{J}}^{\prime }{\mathit{X}}^{\prime }{\mathbf{\Omega }}^{-1}\mathit{X}\mathit{\beta },$$

(13)

in which ${\mathit{J}}^{\prime }=({\mathit{I}}_k,{\mathit{I}}_k,\dots ,{\mathit{I}}_k)$ is a $k\times mk$ matrix, where ${\mathit{I}}_k$ denotes the $k\times k$ identity matrix.

Equivalently, the parameter homogeneity assumption can be formulated as a restriction matrix as

$$\left[\begin{array}{c}{\mathit{\beta }}_1-\overline{\mathit{\beta }}\\ {\mathit{\beta }}_2-\overline{\mathit{\beta }}\\ ⋮\\ {\mathit{\beta }}_m-\overline{\mathit{\beta }}\end{array}\right]=\mathit{\beta }-\mathit{J}\overline{\mathit{\beta }}=({\mathit{I}}_{mk}-\mathit{J}{({\mathit{J}}^{\prime }{\mathit{X}}^{\prime }{\mathbf{\Omega }}^{-1}\mathit{X}\mathit{J})}^{-1}{\mathit{J}}^{\prime }{\mathit{X}}^{\prime }{\mathbf{\Omega }}^{-1}\mathit{X})\mathit{\beta }=\mathit{R}\mathit{\beta }=\mathbf{0},$$

(14)

where $\mathit{R}={\mathit{I}}_{mk}-\mathit{J}{({\mathit{J}}^{\prime }{\mathit{X}}^{\prime }{\mathbf{\Omega }}^{-1}\mathit{X}\mathit{J})}^{-1}{\mathit{J}}^{\prime }{\mathit{X}}^{\prime }{\mathbf{\Omega }}^{-1}\mathit{X}$ is an idempotent matrix.

Hence, we can derive the restricted estimator from the following minimization

$$\underset{\mathrm{s}.\mathrm{t}.\mathit{\beta }}{\mathrm{Minimize}}{(\mathit{y}-\mathit{X}\mathit{\beta })}^{\prime }{\mathbf{\Omega }}^{-1}(\mathit{y}-\mathit{X}\mathit{\beta }),\mathrm{subject}\mathrm{to}\mathit{R}\mathit{\beta }=\mathbf{0}.$$

The solution to the above minimization can be formulated as a feasible restricted GLS estimator in below

$$\begin{array}{ll}\hfill \tilde{\mathit{\beta }}& =\widehat{\mathit{\beta }}-{({\mathit{X}}^{\prime }{\widehat{\mathbf{\Omega }}}^{-1}\mathit{X})}^{-1}{\widehat{\mathit{R}}}^{\prime }{\left[\widehat{\mathit{R}}{({\mathit{X}}^{\prime }{\widehat{\mathbf{\Omega }}}^{-1}\mathit{X})}^{-1}{\widehat{\mathit{R}}}^{\prime }\right]}^{-1}\widehat{\mathit{R}}\widehat{\mathit{\beta }}=({\mathit{I}}_{mk}-\widehat{\mathit{R}})\widehat{\mathit{\beta }}\hfill \\ \hfill & =\mathit{J}{({\mathit{J}}^{\prime }{\mathit{X}}^{\prime }{\widehat{\mathbf{\Omega }}}^{-1}\mathit{X}\mathit{J})}^{-1}{\mathit{J}}^{\prime }{\mathit{X}}^{\prime }{\widehat{\mathbf{\Omega }}}^{-1}\mathit{X}\widehat{\mathit{\beta }},\hfill \end{array}$$

(15)

where2 $\widehat{\mathit{R}}={\mathit{I}}_{mk}-\mathit{J}{({\mathit{J}}^{\prime }{\mathit{X}}^{\prime }{\widehat{\mathbf{\Omega }}}^{-1}\mathit{X}\mathit{J})}^{-1}{\mathit{J}}^{\prime }{\mathit{X}}^{\prime }{\widehat{\mathbf{\Omega }}}^{-1}\mathit{X}$ is an estimate of $\mathit{R}.$

3.3. Average Estimator

We define the average estimator as below

$${\widehat{\mathit{\beta }}}_A=(1-\frac{\tau }{D})\widehat{\mathit{\beta }}+\frac{\tau }{D}\tilde{\mathit{\beta }},$$

(16)

where D is a weighted quadratic loss function defined as

$$D={(\widehat{\mathit{\beta }}-\tilde{\mathit{\beta }})}^{\prime }\mathit{W}(\widehat{\mathit{\beta }}-\tilde{\mathit{\beta }}),$$

(17)

with $\mathit{W}$ an arbitrary symmetric positive definite weight matrix with elements of order $O(T)$, and $\tau$ is a positive characterizing parameter. We will defer describing the optimal choice for this parameter in the next section.3

The idea behind the average estimator defined above is that when the difference between the restricted and the unrestricted GLS estimators is small (D is small), the average estimator gives a higher weight to the restricted GLS estimator, as it is the most efficient estimator. However, when the difference between the restricted and the unrestricted GLS estimators is substantial, the bias of the restricted GLS estimator, resulting from ignoring the parameter heterogeneity, could be more than its variance efficiency gain, so the average estimator assigns a higher weight to the unrestricted GLS estimator.

4. Large-Sample Approximate Bias and MSE

We employ the large-sample approximations method developed by (Nagar 1959), to analyze the bias, mean squared error matrix (MSEM) and risk of the average estimator.

Theorem 1.

Under Assumptions 1–3, the bias of the average estimator up to order $O(T^{-1})$ is

$$\begin{array}{ll}\hfill Bias({\widehat{\mathit{\beta }}}_A)& =\mathbb{E}({\widehat{\mathit{\beta }}}_A-\mathit{\beta })=-\frac{\tau }{\varphi }\mathit{R}\mathit{\beta },\hfill \end{array}$$

(18)

and the MSEM of the average estimator up to order $O(T^{-2})$ is

$$\begin{array}{ll}\hfill \mathit{MSEM}({\widehat{\mathit{\beta }}}_A)& =\mathbb{E}\left[({\widehat{\mathit{\beta }}}_A-\mathit{\beta }){({\widehat{\mathit{\beta }}}_A-\mathit{\beta })}^{\prime }\right]=\mathit{MSEM}(\widehat{\mathit{\beta }})+\frac{{\tau }^2}{{\varphi }^2}\mathit{R}\mathit{\beta }{\mathit{\beta }}^{\prime }{\mathit{R}}^{\prime }-\frac{2\tau }{\varphi }\mathit{R}{({\mathit{X}}^{\prime }{\mathbf{\Omega }}^{-1}\mathit{X})}^{-1}{\mathit{R}}^{\prime }\hfill \\ & +2\frac{\tau }{{\varphi }^2}[\mathit{R}\mathit{\beta }{\mathit{\beta }}^{\prime }{\mathit{R}}^{\prime }\mathit{W}\mathit{R}{({\mathit{X}}^{\prime }{\mathbf{\Omega }}^{-1}\mathit{X})}^{-1}{\mathit{R}}^{\prime }+\mathit{R}{({\mathit{X}}^{\prime }{\mathbf{\Omega }}^{-1}\mathit{X})}^{-1}{\mathit{R}}^{\prime }\mathit{W}\mathit{R}\mathit{\beta }{\mathit{\beta }}^{\prime }{\mathit{R}}^{\prime }],\hfill \end{array}$$

(19)

and for the symmetric positive definite weight matrix $\mathit{W}$ of order $O(T)$, the risk of the average estimator up to order $O(T^{-1})$ is

$$\begin{array}{ll}\hfill \mathit{Risk}({\widehat{\mathit{\beta }}}_A)& =\mathbb{E}\left[{({\widehat{\mathit{\beta }}}_A-\mathit{\beta })}^{\prime }\mathit{W}({\widehat{\mathit{\beta }}}_A-\mathit{\beta })\right]=\mathit{Risk}(\widehat{\mathit{\beta }})+\frac{\tau }{\varphi }2[\tau -2[\mathit{tr}(\mathit{P})-2\frac{{\varphi }_p}{\varphi }]\hfill \\ \hfill & \le \mathit{Risk}(\widehat{\mathit{\beta }})+\frac{\tau }{\varphi }2[\tau -2[\mathit{tr}(\mathit{P})-2{\lambda }_{max}(\mathit{P})],\hfill \end{array}$$

(20)

where $\varphi ={\mathit{\beta }}^{\prime }{\mathit{R}}^{\prime }\mathit{W}\mathit{R}\mathit{\beta }=O(T),$${\varphi }_p={\mathit{\beta }}^{\prime }{\mathit{R}}^{\prime }{\mathit{W}}^{1/2}\mathit{P}{\mathit{W}}^{1/2}\mathit{R}\mathit{\beta }=O(T),$ ${\lambda }_{max}(.)$ denotes the maximum eigenvalue, and $\mathit{P}={\mathit{W}}^{1/2}\mathit{R}{({\mathit{X}}^{\prime }{\mathbf{\Omega }}^{-1}\mathit{X})}^{-1}{\mathit{R}}^{\prime }{\mathit{W}}^{1/2}.$

Proof. 

Appendix A.  ☐

We note that,

$$\mathrm{MSEM}(\widehat{\mathit{\beta }})=(1+\frac{m}{T}){({\mathit{X}}^{\prime }{\mathbf{\Omega }}^{-1}\mathit{X})}^{-1}-\frac{1}{T}{({\mathit{X}}^{\prime }{\mathbf{\Omega }}^{-1}\mathit{X})}^{-1}\mathit{H}{({\mathit{X}}^{\prime }{\mathbf{\Omega }}^{-1}\mathit{X})}^{-1}+o(T^{-2}),$$

(21)

see (Srivastava 1970) for a proof, hence

$$\begin{array}{ll}\hfill \mathrm{Risk}(\widehat{\mathit{\beta }})& =(1+\frac{m}{T})\mathrm{tr}\left[\mathit{W}{({\mathit{X}}^{\prime }{\mathbf{\Omega }}^{-1}\mathit{X})}^{-1}\right]-\frac{1}{T}\mathrm{tr}\left[\mathit{W}{({\mathit{X}}^{\prime }{\mathbf{\Omega }}^{-1}\mathit{X})}^{-1}\mathit{H}{({\mathit{X}}^{\prime }{\mathbf{\Omega }}^{-1}\mathit{X})}^{-1}\right]+o(T^{-1}),\hfill \end{array}$$

(22)

where $\mathit{H}={\mathit{X}}^{\prime }({\mathbf{\Sigma }}^{-1}\otimes \mathbf{\Phi })\mathit{X}-{\mathit{X}}^{\prime }\mathbf{\Pi }\mathit{X}.$

From Theorem 1, it follows that the average estimator dominates the unrestricted GLS estimator in terms of having a smaller risk, when the second term on the right-hand side of Equation (20) is negative, which will hold when

$$0<\tau <2[\mathrm{tr}(\mathit{P})-2\frac{{\varphi }_p}{\varphi }],$$

(23)

given $\mathrm{tr}(\mathit{P})>2{\varphi }_p/\varphi .$ As the upper bound of the condition above depends on the slope parameters, one could replace it with an infimum value. Let d be defined as

$$d=\frac{\mathrm{tr}(\mathit{P})}{{\lambda }_{max}(\mathit{P})},$$

(24)

which lies in the range $d\in [0,(m-1)k],$ as $\mathit{P}$ is a non-zero positive semi-definite matrix. Therefore, when $d>2$, an infimum value for the upper bound is $2[\mathrm{tr}(\mathit{P})-2{\lambda }_{max}(\mathit{P})].$ 4 Therefore, given $d>2,$ an equivalent condition for the condition in (23) can be written as

$$0<\tau \le 2[\mathrm{tr}(\mathit{P})-2{\lambda }_{max}(\mathit{P})].$$

(25)

In other words, when $d>2,$ and $\tau$ satisfies the condition in Equation (25), the risk of the average estimator is less than the risk of the unrestricted GLS estimator up to the order of interest. In addition, as the choice of the characteristic parameter is user-specified, its optimal value, ${\tau }_{opt}$, that minimizes the upper bound of the risk of the average estimator (the last term in Equation (20)), up to order $O(T^{-1})$, provided $d>2$, is

$${\tau }_{opt}=\mathrm{tr}(\mathit{P})-2{\lambda }_{max}(\mathit{P}).$$

(26)

Since the optimal $\tau$ depends on the unknown value of $\mathbf{\Omega }$, one could substitute it with its estimated value, and use an estimate of ${\tau }_{opt}$, as below

$${\widehat{\tau }}_{opt}=\mathrm{tr}(\widehat{\mathit{P}})-2{\lambda }_{max}(\widehat{\mathit{P}}),$$

where $\widehat{\mathit{P}}={\mathit{W}}^{1/2}\widehat{\mathit{R}}{({\mathit{X}}^{\prime }{\widehat{\mathbf{\Omega }}}^{-1}\mathit{X})}^{-1}{\widehat{\mathit{R}}}^{\prime }{\mathit{W}}^{1/2},$ is an estimate of $\mathit{P}$.

Corollary 1.

Under Assumptions 1–3, when $d>2,$ then up to order $O(T^{-1})$ we have

$$\mathit{Risk}({\widehat{\widehat{\mathit{\beta }}}}_A)\le \mathit{Risk}(\widehat{\mathit{\beta }})-\frac{{\tau }_{opt}^2}{\varphi }<\mathit{Risk}(\widehat{\mathit{\beta }}),$$

(27)

where ${\widehat{\widehat{\mathit{\beta }}}}_A$ is the average estimator with ${\widehat{\tau }}_{opt}.$

Proof. 

Appendix A.  ☐

Two arbitrary choices of $\mathit{W}$ are $T{\mathit{I}}_{mk},$ and ${\mathit{X}}^{\prime }{\mathbf{\Omega }}^{-1}\mathit{X},$ where the former one in the risk gives the mean squared error and the latter one, provides the (in-sample) mean squared forecast error (MSFE).

Corollary 2.

Under Assumptions 1–3, when $0<\tau <2[(m-1)k-2]$, then up to order $O(T^{-1})$, we have

$$\mathit{MSFE}({\widehat{\mathit{\beta }}}_A)=\mathit{MSFE}(\widehat{\mathit{\beta }})+\frac{\tau }{\varphi }2[\tau -2[(m-1)k-2]2]<\mathit{MSFE}(\widehat{\mathit{\beta }}).$$

The optimal value of τ that minimizes the MSFE of the average estimator, provided $(m-1)k>2$, is

$${\tau }_{opt,F}=(m-1)k-2,$$

(28)

and the associated optimal MSFE of the average estimator up to order $O(T^{-1})$ is

$${\mathit{MSFE}}_{opt}({\widehat{\mathit{\beta }}}_A)=\mathit{MSFE}(\widehat{\mathit{\beta }})-\frac{{[(m-1)k-2]}^2}{\varphi },$$

(29)

where

$$\mathit{MSFE}(\widehat{\mathit{\beta }})=(1+\frac{m}{T})mk-\frac{1}{T}\mathit{tr}\left[\mathit{H}{({\mathit{X}}^{\prime }{\mathbf{\Omega }}^{-1}\mathit{X})}^{-1}\right].$$

Proof. 

Appendix A. ☐

5. Monte Carlo Simulation

The results below are the simulation results of the model of Section 2, where $x_{it,1}=1$ and the remaining regressors are independently generated from the standard normal distribution. The sample size varies from $T=100,m=3,6,\mathrm{and}k=3,5$, leading to four combinations of m, and k. $u_{1t}$ is generated as $IIDN(0,1),$ while $u_{it}=c{u}_{1t}+v_{it},\mathrm{for}i=2,\dots ,m,$ where $v_{it}\sim IIDN(0,1)$ and $c=0.5.$ We consider two DGPs for generating ${\mathit{\beta }}_i$’s, the first one is under a complete heterogeneity in coefficients where we assume that

$$\mathrm{DGP}1:{\mathit{\beta }}_i=\overline{\mathit{\beta }}+(i\times \delta )/m,i=1,2,\dots ,m,$$

with $\overline{\mathit{\beta }}={(1,1,\dots ,1)}^{\prime },$ and the second DGP is under a weak heterogeneity where we assume that

$$\mathrm{DGP}2:{\mathit{\beta }}_{i1},{\mathit{\beta }}_{i2}=\left\{\begin{array}{cc}1+(i\times \delta )/m,\hfill & \mathrm{if}i=1,\dots ,[m/2]\hfill \\ 1.2,\hfill & \mathrm{if}i=[m/2]+1,\dots ,m\hfill \end{array}\right.,{\mathit{\beta }}_{il}=2,l\in \{3,\dots ,k\},$$

where $[m/2]$ denotes the largest integer value that is smaller than $m/2,$ and $\delta$ takes values on a 10-point grid on $[0,1]$.

The results of 1000 monte carlo simulations are given in Figure 1, Figure 2, Figure 3 and Figure 4, where the vertical axes measure the relative mean squared error (RMSE) of the unrestricted GLS estimator, the restricted GLS estimator and the average estimator to the unrestricted GLS estimator. Hence, the RMSE of the unrestricted GLS estimator is equal to one. The horizontal axes measure the degree of parameter heterogeneity, $\delta$, which is set between zero and one with $0.1$ grid value.

The Monte Carlo results support our theoretical findings of the previous section. The figures show that the RMSE of the average estimator for the whole parameter heterogeneity is below that of the unrestricted estimator. This shows the superiority of the average estimator relative to the unrestricted GLS estimator.

The RMSE of the average estimator in DGP1 of a complete heterogeneous SUR, is smaller than that of the restricted GLS estimator except for very small values of parameter heterogeneity ($\delta$). This is expected because as $\delta$ takes higher values, the bias of the restricted GLS estimator increases, which then results in higher MSE. In DGP2 where the SUR is characterized by some degrees of homogeneity, the RMSE of the restricted GLS estimator remains smaller than that of the unrestricted GLS estimator for larger values of $\delta$ relative to DGP1. In this case, the unrestricted GLS estimator can be inferior to the restricted GLS estimator even with the presence of weak degrees of heterogeneity. This is because although the unrestricted GLS estimator is unbiased, it is inefficient, especially under small sample sizes, and a high number of regressors. In contrast, the restricted GLS estimator properly makes the use of cross equation variations and thus provides a more accurate results.

In general, we find that the average estimator performs robustly well in SUR with various degrees of heterogeneity. When there is a strong heterogeneity, the average estimator prevails. When there is a relatively weak heterogeneity, the average estimator tends to gain more from the efficiency of the restricted GLS estimator by assigning a high weight to this estimator and thus still remains one of the best choices.

6. Application: Returns to Scale in US Banking Industry

In this section we apply the average estimator studied in the previous sections to regressions of the cost system for U.S. commercial banks. We are interested in estimating the returns to scale (RTS) for these banks over the past recent years.

Over the past few years, the number of U.S. commercial banks fell by almost 70%, where in 1984 the total number of U.S. commercial banks was 14,391 and dropped to 4773 in 2018. Over the same period of time, the average asset value of U.S. banks (adjusted for inflation), which is also a measure of bank size, increased by about ten times, from 140 million dollars in 1984 to 1400 million dollars in 2018 (See Figure 5). To support this bank size expansion, bank executives and analysts claim that due to the changes in regulation (such as the permission of interstate branching and combination of banks) and because of technological and financial innovations (such as communication technologies, the securitization and sale of bank loans) over the past few years, the cost of production for larger banks has reduced and encouraged banks to grow larger and/or merge.

On the other hand, critics contend that this decrease in the number of operating banks, and having banks with large assets not only impact the market competition, but also result in agency problems and disproportionate benefits of government policies in favor of large banks. In particular, the financial crisis of 2007 focused attention on large financial institutions considered as “too-big-to-fail”. These together have brought attention of policy makers for regulatory limits on bank size. However, any policy intervention needs to consider the potential efficiency benefits of operating at a large scale. Therefore, estimation of scale economies and RTS is essential for analyzing the costs and benefits of any policy intervention to control the size of banks.

The estimation of scale economies and RTS for U.S. banking industry has stimulated a substantial body of studies. Older empirical studies that used data from the 1980s and 1990s did not find scale economies in banking industry except for very small banks. But recent research that used data from the 2000s, and more modern methods for estimating the banking models, has found considerably more evidence of scale economies in banking. These studies include (Hughes et al. 1996, 2000, 2001; Berger and Mester 1997; Hughes and Mester 1998, 2013; Wheelock and Wilson 2001; Feng and Serletis 2009). Most of the studies in the literature partition banks based on their asset size in groups and estimate each group independently from the other groups. However, not only there is no reason to believe these categories are set based on banks underlying technology, at the same time, it is hard to believe that these groups are not affected by some unknown factors that could have resulted in correlations between groups. Therefore, there are two issues that need to be carefully considered by researchers. First, estimating cost efficiency using all observations as one group has the advantage of smaller variance but at the same time, it means ignoring the potential heterogeneity bias due to difference in technology. Second, partitioning banks in different groups and using single-equation estimators are inefficient compared to pooled estimators which ignore heterogeneity. Hence, there is a trade-off between bias and variance efficiency between these two estimators. As the average estimator introduced in the previous sections results in the optimal balance between bias and variance efficiency, we recommend using this estimator in the estimation of the returns to scale for banking industry to obtain robust and efficient estimators.

6.1. The Model

We follow the so called “intermediation approach” framework of (Sealey and Lindley 1977), which is broadly employed in the literature. According to this approach, a bank’s balance sheet is assumed to capture the essential structure of a bank’s core business. Inputs are considered to be liabilities (core deposits and purchased funds), physical capital and labor. The inputs result in the bank’s productions which are assets (other than the physical, includes loans and trading securities).

With regard to variable specification, we define five inputs and five outputs that are the ones used in the literature. We define the following output quantities: consumer loans ($y_1$), real estate loans ($y_2$), loans to business and other institutions ($y_3$), federal funds sold and securities purchased under agreements to resell ($y_4$), and other assets ($y_5$). The input variables are: labor quantities ($x_1$), premises and fixed assets ($x_2$), purchased funds ($x_3$), interest-bearing transaction accounts ($x_4$), and non-transaction accounts ($x_5$). For each input $x_j$, its price $w_j$ is obtained by dividing its total expenses by the corresponding input quantities.

For modeling the cost of banking industry, we consider a translog cost function and normalize it, so that the homogeneity (in input prices) property is automatically satisfied. We allow for individual (fixed) effects by adding intercepts in each regression, to control for specific group characteristics, heterogeneity in skills and so on. Hence, the cost equation for each group $i=1,2,\dots ,m,$ is considered as

$$\begin{array}{l}\ln (C_{it}/w_{5,it})={\beta }_{0,i}+{\displaystyle \sum _{n=1}^4}{\beta }_{n,i}\ln (w_{n,it}/w_{5,it})+{\displaystyle \sum _{l=1}^5}{\gamma }_{l,i}\ln (y_{l,it})+\frac{1}{2}{\displaystyle \sum _{n=1}^5}{\displaystyle \sum _{l=1}^5}{\gamma }_{nl,i}\ln (y_{n,it})\ln (y_{l,it})\hfill \\ +\frac{1}{2}{\displaystyle \sum _{n=1}^4}{\displaystyle \sum _{l=1}^4}{\eta }_{nl,i}\ln (w_{n,it}/w_{5,it})\ln (w_{l,it}/w_{5,it})+{\displaystyle \sum _{l=1}^4}{\displaystyle \sum _{n=1}^5}{\delta }_{ln}\ln (w_{l,it}/w_{5,it})\ln (y_{n,it})+u_{it},t=1,2,\dots ,T_i,\hfill \end{array}$$

(30)

where $T_i$ is the number of observations in group i (the number of banks operating within group i), and $C_{it}$ is the total cost of bank t, in group i, defined as

$$C_{it}=w_{1,it}{x}_{1,it}+w_{2,it}{x}_{2,it}+w_{3,it}{x}_{3,it}+w_{4,it}{x}_{4,it}+w_{5,it}{x}_{5,it},t=1,2,\dots ,T_i,\mathrm{and}i=1,\dots ,m.$$

(31)

The cost function is symmetric which requires the imposition of the following restrictions on the parameters as below

$$\begin{array}{ll}\hfill {\eta }_{jq,i}& ={\eta }_{qj,i}\hfill \\ \hfill {\gamma }_{nl,i}& ={\gamma }_{ln,i}.\hfill \end{array}$$

(32)

RTS is defined as the inverse of the sum of cost elasticities. If we define the output elasticity of the model for output j of bank t in group i, as $Ec{y}_{j,it}=\frac{\partial \ln (C_{it})}{\partial \ln (y_{j,it})}$ and the sum of cost elasticities as $Ec{y}_{it}={\sum }_{j=1}^5\frac{\partial \ln (C_{it})}{\partial \ln (y_{j,it})},$ then RTS of bank t in group i is defined as $RT{S}_{it}=1/Ec{y}_{it}.$ Also, the RTS for group i is a vector of $T_i\times 1$ defined as

$$RT{S}_i=\left[\begin{array}{c}RT{S}_{i1}\\ RT{S}_{i2}\\ ⋮\\ RT{S}_{i{T}_i}\end{array}\right],$$

(33)

which is used for calculating mean, quartiles, and deciles of RTS for group i with $T_i$ banks, see Section 6.3.

A bank with $RTS>1$, has increasing returns to scale, that is for one percent increase in all outputs, cost is increased by less than one percent, and the bank is operating below its efficient scale size ($RTS=1$) when $RTS<1$.

6.2. The Data

The data we use is obtained from the Reports of Income and Condition (Call Reports)5, over the period from 2000 to 2018. We omit observations where negative values for assets, equity, outputs, and prices are reported. The summary of the data for years 2000 and 2018 is in

Table 1. Summary Statistics.

Variable	Min		Max
	2000	2018	2000	2018
C	136.49	151.78	18,369,359.04	26,757,643.39
$w_1$	7015.14	15,321.43	163,358.47	261,240.11
$w_2$	6.75	9.61	12375.00	43,469.79
$w_3$	3.20	0.02	206.53	48.82
$w_4$	0.26	0.02	265.36	153.96
$w_5$	0.46	0.06	49.21	16.82
$y_1$	1.25	1.00	42,638,250.00	173,922,000.00
$y_2$	0.25	170.50	168,465,250.00	686,161,250.00
$y_3$	45.75	125.50	178,056,500.00	457,517,750.00
$y_4$	89.50	0.25	144,188,250.00	703,099,250.00
$y_5$	39.00	42.50	86,346,000.00	704,384,250.00
Variable	Mean		STD
	2000	2018	2000	2018
C	24,454.49	43,803.47	331,091.51	643,872.50
$w_1$	26,175.32	49,386.16	7084.42	15,332.39
$w_2$	230.31	265.78	351.89	920.70
$w_3$	33.22	3.51	6.38	2.75
$w_4$	14.79	2.43	9.01	3.87
$w_5$	22.87	3.46	4.79	1.90
$y_1$	57,014.23	268,717.11	754,295.15	4,722,816.78
$y_2$	175,843.57	881,915.93	2,961,928.59	15,632,429.24
$y_3$	206,836.22	929,043.01	2,423,963.16	10,526,805.47
$y_4$	192,820.92	862,701.61	2,985,314.66	16,392,980.08
$y_5$	85,019.93	842,793.39	1,554,035.20	16,868,329.90

Note: All variables except $w_1$–$w_5,$ are measured in thousands of dollars.

6.

Following (Feng and Serletis 2009) and others in the literature, we classify the banks into three groups which is mainly based on the standard asset size categories that are used by the Federal Financial Institutions Examination Council (FFIEC). Banks with over $\$500$ million in total assets are classified as Large banks, banks with assets between $\$100$ million and $\$500$ million are classified as Medium banks, and banks with under $\$100$ million in assets are classified as Small banks. In order to have a consistent partitions over time, the asset size caps in each year are justified upward by the growth in the CPI.

Table 2. Banks Asset Size Classes.

Bank Groups	Asset Size (in Millions of Dollars of the Year)	Number of Banks	Share of Banks
	2000
Large Banks	Assets ≥ 500	739	8.9
Medium Banks	100 ≤ Assets < 500	2946	35.5
Small Banks	Assets < 100	4620	55.6
	2018
Large Banks	Assets ≥ 729	988	19.8
Medium Banks	146≤ Assets <729	2298	46.1
Small Banks	Assets < 146	1704	34.1

presents the number and share of banks in each group with their corresponding asset ranges for years 2000 and 2018.

6.3. Estimation

We estimate model of Equation (30) using the average estimation method developed in the previous sections for each year separately. Basically, our SUR at each year consists of three cost equations representing Large, Medium, and Small bank groups, and the observations for each regression are the operating banks data under each bank group. Since the sample size for each group is different, we face a SUR with unequal number of observations and to estimate the variance-covariance matrix ($\mathbf{\Omega }$), we consider the following procedures recommended in the literature (See Schmidt 1977; Baltagi et al. 1989):

• Ignore the extra observations in estimating $\mathbf{\Omega }$;
• Use the extra observations to estimate variances. This procedure has the disadvantage of producing estimates of $\mathbf{\Omega }$ that are not positive definite;
• Use the extra observations to estimate variances, and modifying the estimates of covariances using the method of (Srivastava and Zaatar 1973);
• Use all observations in estimation, following the method of (Hocking and Smith 1968).

It is known in the literature that the results of the above procedures are much the same. Likewise, we find that the procedures above, generate similar results, so we only report the results of method 3.

After estimating Equation (30), we obtain the sum of cost elasticities for each bank by

$$Ec{y}_{it}=\sum _{j=1}^5\frac{\partial \ln (C_{it})}{\partial \ln (y_{j,it})}=\sum _{j=1}^5{\gamma }_{j,i}+\sum _{j=1}^5\sum _{n=1}^5{\gamma }_{nj,i}\ln (y_{n,it})+\sum _{j=1}^5\sum _{n=1}^4{\delta }_{nj}\ln (w_{n,it}/w_{5,it}),$$

(34)

where $Ec{y}_{it}$ is the sum of cost elasticity of bank t in group i, and the parameters are replaced with their estimated values. Then, the RTS is calculated following Equation (33). We also obtain the RTS using the unrestricted GLS estimator, and the restricted GLS estimator.

The results of years 2000, 2009 and 20187 are reported in

Table 3. Summary of estimates of returns to scale.

Bank Size	Estimator	Estimates of RTS
		${\mathit{D}}_{10}$	${\mathit{Q}}_{25}$	${\mathit{Q}}_{50}$	${\mathit{Q}}_{75}$	${\mathit{D}}_{90}$	Mean
	2000
	Restricted	0.979	0.988	0.997	1.008	1.019	0.998
Large Banks	Unrestricted	0.971	0.985	0.998	1.009	1.021	0.996
	Average Estimator	0.980	0.989	0.998	1.007	1.018	0.998
Medium Banks	Unrestricted	0.969	0.992	1.017	1.043	1.071	1.019
	Average Estimator	0.995	1.002	1.011	1.020	1.028	1.011
Small Banks	Unrestricted	0.997	1.014	1.034	1.055	1.075	1.035
	Average Estimator	1.006	1.014	1.024	1.034	1.044	1.025
	2009
	Restricted	1.007	1.016	1.026	1.038	1.051	1.028
Large Banks	Unrestricted	1.001	1.025	1.053	1.080	1.108	1.055
	Average Estimator	1.004	1.026	1.051	1.075	1.101	1.052
Medium Banks	Unrestricted	0.956	0.990	1.032	1.079	1.125	1.038
	Average Estimator	0.961	0.993	1.031	1.076	1.118	1.037
Small Banks	Unrestricted	0.979	1.013	1.049	1.085	1.124	1.050
	Average Estimator	0.983	1.015	1.048	1.082	1.118	1.049
	2018
	Restricted	1.009	1.021	1.034	1.049	1.066	1.037
Large Banks	Unrestricted	1.009	1.052	1.100	1.140	1.184	1.098
	Average Estimator	1.010	1.052	1.094	1.132	1.172	1.093
Medium Banks	Unrestricted	0.915	0.950	0.988	1.030	1.070	0.989
	Average Estimator	0.919	0.953	0.989	1.029	1.068	0.990
Small Banks	Unrestricted	0.960	0.997	1.042	1.104	1.166	1.054
	Average Estimator	0.962	0.998	1.040	1.098	1.158	1.051

Note: Decile, quartile and mean estimates of RTS for the restricted, unrestricted, and average estimators. $D_{10}$, $D_{90}$, $Q_{25}$, $Q_{50}$, and $Q_{75}$ are the lower decile, the upper decile, the lower quartile, median and upper quartile, respectively.

, which presents the extreme deciles, quartiles and means of our estimated RTS using the three estimation methods. The patterns of mean, extreme deciles and quartiles for all years are plotted in Figure 6. The results base on the restricted estimator over the most recent years suggest increasing RTS at each decile. However, these results are not economically reasonable. On the other hand, we find evidence of decreasing RTS for almost 50% of Small and Medium banks using the unrestricted estimator over the most recent years. These contradicting results show the importance of the average estimator which is used to respond to this model uncertainty.

Comparing the RTS of banks using the average estimator over the sample period shows that, on average majority of banks have increasing returns to scale. In most recent years, the results exhibit more signs of cost efficiency for Large and Small banks, such that all of Large banks have increasing RTS and only less than 25% of Small banks have exhausted their cost efficiency. However, we find that more than 50% of Medium banks are operating under decreasing returns to scale near the end of the sample. The results are consistent with some recent studies (e.g., References Feng and Serletis 2009; Hughes and Mester 2013; Wheelock and Wilson 2001; Henderson et al. 2015; Mailkov et al. 2015) although we are not aware of any study from 2011 to 2018.

7. Conclusions

In this paper, we introduce an averaging estimator for a SUR model. The introduced estimator is a weighed average of an unrestricted GLS estimator which is the (Zellner 1962) estimator and a restricted GLS estimator. The weight is inversely related to a quadratic loss function which measures the weighted distance between the unrestricted and the restricted GLS estimators. The bias, MSE matrix, and risk of the average estimator using the large-sample approximations of (Nagar 1959) are derived. The superiority conditions of the average estimator in terms of the weighted mean squared error is given for any user-specific symmetric positive definite weight matrix, and is not limited to the case where the weight is the inverse of the variance-covariance matrix of the unrestricted GLS estimator.

We also provide some Monte Carlo results which support our theoretical claims. Finally, as our estimator is motivated by economic theory, we use U.S. Commercial banking data, and estimate a cost system for the banking industry to show how our estimator can be used in the applied work. We also estimate the cost system using single-equation least squares and a pooled estimator, and compare them with our proposed average estimator. We found more reliable estimation results with the cost system using our average estimator than the other estimators. We found that on average majority of banks have been operating under increasing returns to scale over the sample period from 2000 to 2018.