# Specification tests for normal/gamma and stable/gamma stochastic frontier
models based on empirical transforms

Christos K. Papadimitriou111Department of Economics, National and Kapodistrian
University of Athens, Greece  , Simos G. Meintanis††footnotemark:  222Pure and
Applied Analytics, North–WestUniversity, South Africa 333Correspondence:Simos
G. Meintanis (on sabbatical leave from the University of Athens),
<EMAIL_ADDRESS>, Mike G. Tsionas444Department of Economics, Lancaster
University, UK

###### Abstract

Goodness–of–fit tests for the distribution of the composed error term in a
Stochastic Frontier Model (SFM) are suggested. We focus on the case of a
normal/gamma SFM and the heavy–tailed stable/gamma SFM. In the first case we
use the moment generating function as tool while in the latter case the
characteristic function of the error term and formulate our test statistics as
weighted integrals of properly standardized data. The new tests are
consistent, and for the normal/gamma SFM are shown to have an intrinsic
relation to moment–based tests. The finite–sample behavior of a resampling
version of the new test for the normal/gamma SFM is investigated by Monte
Carlo simulations, while the corresponding test for the stable/gamma SFM is
applied on several real–data sets.
Keywords: Goodness–of–fit tests; Moment generating function; Characteristic
function; Normal/gamma distribution; Stable distribution

##  1 Introduction

Let $Y$ denote the maximum log–output obtainable from a (row) vector of
log–outputs $X$, and consider the stochastic frontier production model (SFM),

$Y=X\beta+\varepsilon,$ (1)

where $\beta$ is an unknown (column) vector of parameters, and $\varepsilon$
denotes the error term. As it is standard practice in the context of SFMs, the
error term is decomposed as $\varepsilon=v-u$, where the random component $v$
is intended to capture purely random effects, while $u\geq 0$ captures the
technical inefficiency of the underlying firm. In this connection, the
exponential and half/normal distribution are the most commonly used laws for
the one-sided error component $u$. However more complicated distributions,
such as the Gamma distribution have also been studied.

Earlier work on goodness–of–fit tests for certain aspects of the SFM include
Schmidt and Lin (1984), Coelli (1995), Lee (1983), Kopp and Mullahy (1990),
and Bera and Mallick (2002). The work herein however is closer to Wang et al.
(2011) because as in that paper we also suggest omnibus tests, i.e. tests that
have non–trivial power against arbitrary deviations from the null hypothesis,
and use resampling in order to capture the extra variability introduced at the
estimation of parameters step preceding the test. We also share common tools
with Chen and Wang (2012) who use the characteristic function for testing
distributional specifications in the context of SFMs. Herein we follow the
approach of Meintanis and Papadimitriou (2022) in suggesting specification
tests for SFMs based on distributional transforms such as the moment
generating function (MGF) and the characteristic function (CF). In this
connection we introduce a two–level generalization of the methods of Meintanis
and Papadimitriou (2022), first from the normal/exponential SFM to the
normal/gamma SFM using the MGF, and then to the even more general stable/gamma
SFM. Because the latter model is heavy–tailed we employ the CF rather than the
MGF, both in the estimation step as well as in formulating our tests. We
finally mention a significant most recent contribution by Cheng et al. (2022)
that propose versions of the classical distribution–function based
Kolmogorov–Smirnov and Cramér–von Mises goodness–of–fit tests for the law of
the inefficiency component $u$. Their tests are more general than ours in that
they do not assume a specific parametric form for the production function and
that they are not tailored, but apply to any distribution under test and allow
for heteroskedasticity. On the other hand the tests of Cheng et al. (2022) can
not be applied to our setting because they assume that all other parameters
besides scale are fixed, while our tests are fully composite in the sense that
all distributional parameters are estimated from the data at hand. In our
methods we employ as tools the MGF and the CF which are much more easier to
apply in the situations considered herein, namely for a normal/gamma SFM as
well as for the stable/gamma SFM because the corresponding distribution
functions are quite complicated. A further advantage of the methods suggested
herein is that they allow one to distinguish between the former model
involving the normal distribution and a purely heavy–tailed stable/gamma SFM.

The rest of the paper runs as follows. In Section 2 we introduce the test
statistic for the normal/gamma SFM, show its consistency and discuss the role
of the weight function. In Section 3 we consider estimation of model
parameters and obtain an interesting limit statistic. Simulation results for a
resampling version of the test for the normal/gamma null hypothesis are
presented in Section 4. In Section 5 we consider the more general problem of
testing the stable/gamma SFM, while Section 6 contains application of the
proposed test to several data sets from the literature. Finally, conclusions
and outlook are presented in Section 7. Some technical material is deferred to
the Appendices A.1–A.5.

##  2 Test for the normal/gamma SFM

###  2.1 Test statistic

Consider the SFM in equation (1) where $\varepsilon=v-u$, where $u$ and $v$
are distributed independently of each other. In this context and on the basis
of data $(X_{j},Y_{j}),\ j=1,...,n$, we wish to test the null hypothesis

$\begin{split}{\cal{H}}_{0}:\;&\textrm{Model (\ref{SFM}) holds true with
}v\sim{\cal{N}}(0,\sigma_{v}^{2}),\ u\sim{\cal{G}}(p,c),\textrm{ for some
}\sigma_{v}^{2},p,c>0,\end{split}$ (2)

where ${\cal{N}}(\mu,\sigma^{2})$, denotes a normal distribution with mean
$\mu$ and variance $\sigma^{2}$, and ${\cal{G}}(p,c)$ denotes a gamma
distribution with shape parameter $p$ and scale parameter $c$.

Recall in this connection that the moment generating function (MGF) of a
generic random variable, say $Z$, is defined by
$M_{Z}(t)=\mathbb{E}({e^{tZ}}),\>t\in\mathbb{R}$, while if $Z\geq 0$, then its
Laplace transform is given by $L_{Z}(t)=\mathbb{E}(e^{-tZ}),\ t>0$. For later
use we also mention that if $Z\sim{\cal{N}}(\mu,\sigma^{2})$ with density
$(\sigma\sqrt{2\pi})^{-1}e^{-(1/2)((z-\mu)/\sigma)^{2}}$ then its MGF is given
by $M_{Z}(t)=e^{t\mu+\frac{1}{2}\sigma^{2}t^{2}}$, while if
$Z\sim{\cal{G}}(p,c)$ with density $z^{p-1}/(c^{p}{\rm{\Gamma}}(p))e^{-z/c},\
z>0$, then $L_{Z}(t)=(1+ct)^{-p}$, where
${\rm{\Gamma}}(p)=\int_{0}^{\infty}z^{p-1}e^{-z}{\rm{d}}z$ is the gamma
function.

Based on the independence of $u$ and $v$, it readily follows that under the
null hypothesis ${\cal{H}}_{0}$, the MGF of the composed error term
$\varepsilon$ may be computed as

$M_{\varepsilon}(t)=M_{v}(t)L_{u}(t)=\frac{e^{\frac{1}{2}\sigma_{v}^{2}t^{2}}}{(1+ct)^{p}},\
t>0,$ (3)

and thus satisfies the differential equation

$(1+ct)M^{\prime}_{\varepsilon}(t)+[cp-\sigma^{2}_{v}t(1+ct)]M_{\varepsilon}(t)=0;$
(4)

see also Henze et al. (2012).

Consequently the MGF of the standardized composed error term
$\widetilde{\varepsilon}=\varepsilon/c$ is given by

$M_{\widetilde{\varepsilon}}(t)=\frac{e^{\frac{1}{2}\lambda
t^{2}}}{(1+t)^{p}},$ (5)

and thus satisfies the differential equation

$(1+t)M^{\prime}_{\widetilde{\varepsilon}}(t)+[p-\lambda
t(1+t)]M_{\widetilde{\varepsilon}}(t)=0,\ t>0,$ (6)

where $\lambda=\frac{\sigma^{2}_{v}}{c^{2}}$.

Now let $\widehat{\varepsilon}_{j}=Y_{j}-\widehat{\beta}X_{j}$, be the
residuals of the SFM and write $r_{j}=\widehat{\varepsilon}_{j}/\widehat{c},\
j=1,\ldots,n$, for the respective standardized residuals, involving suitable
parameter estimates
($\widehat{\beta},\widehat{\sigma}_{v}^{2},\widehat{c},\widehat{p}$), obtained
under the null hypothesis ${\cal{H}}_{0}$. Then the left– hand side of
equation (6) may be estimated by

$D_{n}(t):=(1+t)M^{\prime}_{n}(t)+[\widehat{p}-\widehat{\lambda}t(1+t)]M_{n}(t),\
t>0,$ (7)

where

$M_{n}(t)=\frac{1}{n}\sum_{j=1}^{n}e^{tr_{j}},\
M^{\prime}_{n}(t)=\frac{1}{n}\sum_{j=1}^{n}r_{j}e^{tr_{j}},$ (8)

stand for the empirical MGF based on $r_{j}$ and its derivative, respectively.
In view of (6), we should have for large $n$, $D_{n}(t)\approx 0$ under the
null hypothesis, identically in $t$. Thus it is reasonable to reject the null
hypothesis ${\cal{H}}_{0}$ of a normal/gamma SFM for large values of the test
statistic

$T_{n,w}=n\int_{0}^{\infty}D^{2}_{n}(t)\;w(t)\;{\rm{d}}t,$ (9)

where $w(t)>0$ is an integrable weight function.

Due to the uniqueness of the solution of the differential equation figuring in
(4) and the uniqueness of MGFs, it follows that the test statistic $T_{n,w}$
defined by (9) satisfies

$T_{n,w}\rightarrow\infty,\ \mbox{a.s. as}\ n\to\infty,$

under alternatives, and thus the test which rejects the null hypothesis
${\cal{H}}_{0}$ for large values of $T_{n,w}$ has asymptotic power one against
arbitrary deviations from ${\cal{H}}_{0}$. The argument leading to this
property may be found in Appendix A.1.

The choice $w(t)=e^{-\gamma t^{2}},\gamma>0$, is particularly appealing from
the computational point of view as it leads after straightforward algebra to
the explicit formula

$\begin{split}T_{n,\gamma}&=\frac{1}{n}\sum_{j,k=1}^{n}\Bigg{\\{}\frac{r_{j}r_{k}+r_{j}(\widehat{p}-\widehat{\lambda})-\widehat{p}\widehat{\lambda}}{\gamma}+\frac{r_{j}r_{k}+\widehat{\lambda}\big{(}\widehat{\lambda}-4r_{j}-2\widehat{p}\big{)}}{4\gamma^{2}}(r_{j}+r_{k})\\\
&\hskip
50.00008pt+\frac{\widehat{\lambda}(\widehat{\lambda}-r_{j})}{4\gamma^{3}}\big{[}(r_{j}+r_{k})^{2}+4\gamma\big{]}+\frac{\widehat{\lambda}^{2}}{16\gamma^{4}}(r_{j}+r_{k})\big{[}(r_{j}+r_{k})^{2}+10\gamma\big{]}\\\
&\hskip
50.00008pt+\frac{\sqrt{\pi}}{\sqrt{\gamma}}e^{\frac{(r_{j}+r_{k})^{2}}{4\gamma}}\;{\rm{\Phi}}\Big{(}\frac{r_{j}+r_{k}}{\sqrt{2\gamma}}\Big{)}\Bigg{[}(r_{j}r_{k}+2\widehat{p}r_{j}+\widehat{p}^{2})+\frac{r_{j}r_{k}+r_{j}(\widehat{p}-\widehat{\lambda})-\widehat{p}\widehat{\lambda}}{\gamma}(r_{j}+r_{k})\\\
&\hskip
50.00008pt+\frac{r_{j}r_{k}+\widehat{\lambda}\big{(}\widehat{\lambda}-4r_{j}-2\widehat{p}\big{)}}{4\gamma^{2}}\big{[}(r_{j}+r_{k})^{2}+2\gamma\big{]}+\frac{\widehat{\lambda}(\widehat{\lambda}-r_{j})}{4\gamma^{3}}(r_{j}+r_{k})\big{[}(r_{j}+r_{k})^{2}+6\gamma\big{]}\\\
&\hskip
50.00008pt+\frac{\widehat{\lambda}^{2}}{16\gamma^{4}}\big{[}12\gamma^{2}+12\gamma(r_{j}+r_{k})^{2}+(r_{j}+r_{k})^{4}\big{]}\Bigg{]}\Bigg{\\}},\end{split}$
(10)

where ${\rm{\Phi}}(\cdot)$ denotes the standard normal distribution function.
Notice that we write $T_{n,\gamma}$ for the test statistic in (9)
corresponding to the weight function $w(t)=e^{-\gamma t^{2}}$, and thus make
the dependence of the test statistic on $\gamma>0$ explicit. Further details
on the above computational formula for $T_{n,\gamma}$ are provided in Appendix
A.2.

###  2.2 The role of the weight function

Since the parameters of the model are unknown we need to efficiently estimate
them from the data $(X_{j},Y_{j}),\ j=1,\ldots,n$. Maximum likelihood
estimators appear to be the most common, however drawbacks emerge when a
closed form for the density of the random variable is not available as in the
present case; see Tsionas (2012) for alternative estimation methods for SFMs
whose likelihood function is not available in closed form. The corrected least
squares (CLS) estimation is the method employed herein. This method will be
further specified in the next subsection.

At this point, we will illustrate the role that the weight function
$e^{-\gamma t^{2}}$ plays in the test statistic $T_{n,\gamma}$. To this end,
we compute the first five moments of the composed normal/gamma error
$\varepsilon$ in (1) as

$\begin{split}&\mu_{1}=-cp,\mu_{2}=\sigma^{2}_{v}+c^{2}p(p+1),\mu_{3}=-c^{3}p(p+1)(p+2)-3cp\sigma^{2}_{v},\\\
&\mu_{4}=3\sigma^{4}_{v}+c^{4}p(p+1)(p+2)(p+3)+6\sigma^{2}_{v}c^{2}p(p+1),\\\
&\mu_{5}=-15\lambda^{2}p-10\lambda
p(p+1)(p+2)-p(p+1)(p+2)(p+3)(p+4),\end{split}$ (11)

where $\mu_{k}=\mathbb{E}(\varepsilon^{k})$.

From (11) we readily obtain that under the null hypothesis ${\cal{H}}_{0}$ of
a normal/gamma SFM, the moments of the standardized error
$\widetilde{\mu}_{k}=\mathbb{E}(\widetilde{\varepsilon}^{k})$ satisfy the
moment equations

${\cal{M}}_{k}=0,\ k=1,2,3,4,5,$ (12)

where

$\begin{split}&{\cal{M}}_{1}=\widetilde{\mu}_{1}+p,\;{\cal{M}}_{2}=\widetilde{\mu}_{2}+(p+1)\widetilde{\mu}_{1}-\lambda,\;{\cal{M}}_{3}=\widetilde{\mu}_{3}+(p+2)\widetilde{\mu}_{2}-2\lambda(\widetilde{\mu}_{1}+1),\\\
&{\cal{M}}_{4}=\widetilde{\mu}_{4}+(p+3)\widetilde{\mu}_{3}-3\lambda(\widetilde{\mu}_{2}+2\widetilde{\mu}_{1}),\;{\cal{M}}_{5}=\widetilde{\mu}_{5}+(p+4)\widetilde{\mu}_{4}-4\lambda(\widetilde{\mu}_{3}+3\widetilde{\mu}_{2}).\end{split}$
(13)

Now consider the expansion
$e^{z}=1+z+z^{2}/2+z^{3}/6+z^{4}/24+{\rm{o}}(z^{4}),\ z\to 0$, of the
exponential function figuring in (8), and then starting from (7) and after
some long but straightforward algebra we obtain

$D_{n}(t)={\cal{M}}_{1,n}+t{\cal{M}}_{2,n}+\frac{t^{2}}{2}{\cal{M}}_{3,n}+\frac{t^{3}}{6}{\cal{M}}_{4,n}+\frac{t^{4}}{24}{\cal{M}}_{5,n}+{\rm{o}}(t^{4}),$
(14)

$t\to 0$, where ${\cal{M}}_{k,n}$ are the empirical moment equations resulting
from (13) if we replace the population moments $\widetilde{\mu}_{k}$ by the
sample moments $\widetilde{m}_{k}=n^{-1}\sum_{j=1}^{n}r^{k}_{j}$ of the
standardized residuals, $k=1,2,3,4,5$, and the parameters $(\lambda,p,c)$ by
estimators $(\widehat{\lambda},\widehat{p},\widehat{c})$.

A little reflection on (14) illustrates the role of the weight function and of
the tuning parameter $\gamma$. Specifically the weight function provides the
particular functional form of the weights by means of which each empirical
moment equation figuring in (14) enters the test statistic $T_{n,\gamma}$.
Consequently the expansion in (14) underlines the fact that a higher value of
$\gamma$ results in a reduced impact of higher moments and vice versa, smaller
$\gamma$ allows these higher moments to have a more significant impact on
$T_{n,\gamma}$. In this connection caution should be exercised when choosing
$\gamma$: A value that is too large allows only lower moments to enter the
test statistic which in turn will adversely affect the power of
$T_{n,\gamma}$, while on the other hand, a $\gamma$–value which is too small
results in a test statistic prone to numerical error due to divergence of the
exponential terms involved (note that for $\gamma=0$, the value of
$T_{n,\gamma}$ diverges).

##  3 CLS estimation and a limit statistic

As already mentioned we use the CLS estimator; see Kopp and Mullahy (1993).
Specifically estimates of the regression coefficients are obtained by using
the first moment equation in (13), while the estimates of the distributional
parameters $(\sigma^{2}_{v},p,c)$ are obtained by using, in addition to the
first moment equation ${\cal{M}}_{1}$, the next three equations, i.e.
${\cal{M}}_{2},{\cal{M}}_{3}$ and ${\cal{M}}_{4}$, figuring in (13) and
thereby solving the system of four equations so produced (see Appendix A.3).

Referring again to (14), notice that the use of CLS estimators of the
parameters, leads to ${\cal{M}}_{k,n}=0,\ k=1,2,3,4$. Therefore inserting (14)
in the test statistic equation (9), and integrating $D^{2}_{n}(t)$
term–by–term we deduce

$\displaystyle T_{n,\gamma}$ $\displaystyle=$ $\displaystyle
n\left(\frac{{\cal{M}}_{5,n}}{24}\right)^{2}{\cal{I}}_{8,\gamma}+{\rm{o}}(\gamma^{-9/2}),\
\gamma\to\infty,$ (15)

where

${\cal{I}}_{k,\gamma}:=\int_{0}^{\infty}t^{k}e^{-\gamma
t^{2}}{\rm{d}}t=\frac{{\rm{\Gamma}}(\frac{k+1}{2})}{2\gamma^{\frac{k+1}{2}}},$

and consequently

$\displaystyle\lim_{\gamma\to\infty}\frac{2\gamma^{9/2}}{{\rm{\Gamma}}(9/2)}T_{n,\gamma}$
$\displaystyle=$ $\displaystyle n\left(\frac{{\cal{M}}_{5,n}}{24}\right)^{2}.$

The last equation shows that conditionally on CLS estimation, the test
statistic $T_{n,\gamma}$ is dominated by the fifth empirical moment equation,
in the limiting case $\gamma\to\infty$. Also note that, on account of (13) and
as $n\to\infty$, this empirical moment equation vanishes under the null
hypothesis ${\cal{H}}_{0}$ of a SFM with a normal/gamma composed error
distribution.

##  4 Simulations

In this section we present the results of Monte Carlo study for the new test
statistic given by equation (10)555For simulations we used Matlab software
R2020b version.. We consider the normal/gamma SFM, whereby
$\varepsilon\sim{\cal{NG}}(\sigma^{2}_{v},p,c)$ denotes a normal/gamma
distribution with the indicated parameters.

We consider the simple location SFM in equation (1),
$Y_{j}=\beta+\varepsilon_{j}$. The parameter vector
$\Theta:=(\beta,\sigma^{2}_{v},p,c)$ is estimated by the CLS method mentioned
in Section 3. The number of Monte Carlo replications is $M=1,000$, with sample
size $n=50,100,200$, and nominal level of significance $5\%$. For reasons
explained in Section 2.2, we avoided to choose values of the tuning parameter
$\gamma$ close to zero, choosing to implement our simulations for $\gamma=4$,
$\gamma=6$ and $\gamma=8$.

Since the parameters of the model are considered unknown we employ a
parametric bootstrap version of the tests in order to obtain critical values
of the test statistic, say $T$ for simplicity. Specifically, we implement the
so-called “warp speed” bootstrap procedure in order to approximate the
bootstrap critical values in the Monte Carlo study (see Giacomini et al.
(2013) and Chang and Hall (2015)). This procedure is carried out by performing
the following steps:

  * (S1)

Draw a Monte Carlo sample $\\{Y^{(m)}_{j},X^{(m)}_{j}\\},\ j=1,...,n$, compute
the estimator–vector $\widehat{\Theta}^{(m)}$, where
$\widehat{\Theta}^{(m)}=(\widehat{\beta}^{(m)},\widehat{\sigma}^{2(m)}_{v},\widehat{p}^{(m)},\widehat{c}^{(m)})$

  * (S2)

On the basis of $\widehat{\Theta}^{(m)}$ calculate the residuals
$\widehat{\varepsilon}^{(m)}_{j}$ and the corresponding test statistic
$T_{m}=T(\widehat{\varepsilon}^{(m)}_{1},\ldots,\widehat{\varepsilon}^{(m)}_{n})$

  * (S3)

Generate i.i.d. bootstrap errors $\varepsilon_{j}^{(m)},\ j=1,\ldots,n$, where
$\varepsilon_{j}^{(m)}=v_{j}^{(m)}-u_{j}^{(m)}$, with
$v_{j}^{(m)}\sim{\cal{N}}(0,\widehat{\sigma}_{v}^{2})$ and
$u_{j}^{(m)}\sim{\rm{Gamma}}(\widehat{p},\widehat{c})$, and independent.

  * (S4)

Define the bootstrap observations
$Y_{j}^{(m)}=\widehat{\beta}^{(m)}X_{j}+\varepsilon_{j}^{(m)}$, $j=1,...,n$.

  * (S5)

Based on $\\{Y_{j}^{(m)},X_{j}\\}$ compute the bootstrap estimator
$\widehat{\Theta}_{b}^{(m)}=(\widehat{\beta}_{b}^{(m)},\widehat{\sigma}_{b,v}^{2(m)},\widehat{p}_{b}^{(m)},\widehat{c}_{b}^{(m)})$,
and the corresponding bootstrap residuals, say,
$\widehat{\epsilon}_{j}^{(m)}$, $j=1,...,n$.

  * (S6)

Compute the test statistic
$\widehat{T}_{m}:=T(\widehat{\epsilon}^{(m)}_{1},\ldots,\widehat{\epsilon}^{(m)}_{n})$

  * (S7)

Repeat steps (S1)–(S6), for $m=1,...,M$, and collect test–statistic values
$T_{m}$ and bootstrap statistic values $\widehat{T}_{m},\ m=1,...,M$.

  * (S8)

Set the critical point, say $C$, equal to $\widehat{T}_{(M-0.05M)}$, where
$\widehat{T}_{(m)}$, $m=1,\ldots,M$, denote the order statistics corresponding
to $\widehat{T}_{m}$.

  * (S9)

Obtain the empirical rejection rate $M^{-1}\sum_{m=1}^{M}{\bf{1}}_{T_{m}>C}$.

In the Table 1 and Table 2 the size results are presented for the test
statistic $T_{n,\gamma}$, corresponding to the ${\cal{NG}}(1,p,1)$ null
hypothesis for $p=1,2,3$. The power of the test is also obtained against an
alternative mixture hypothesis $0.7{\cal{NG}}(1,1,1)+0.3{\cal{NG}}(1,p,1)$,
whereby an observation is drawn from a ${\cal{NG}}(1,1,1)$ distribution (resp.
${\cal{NG}}(1,p,1)$ distribution) with probability 0.7 (resp. 0.3), for
$p=0.5,2,3$.

Table 1: Size of the test for the normal/gamma null hypothesis ${\cal{NG}}(1,p,1)$ with sample size $n$ | $n$ |  $\gamma$=4 |  $\gamma$=6 |  $\gamma$=8
---|---|---|---|---
$p$=1.0 | 50 | 5.5 | 6.3 | 6.1
| 100 | 5.1 | 4.6 | 4.4
| 200 | 5.2 | 4.7 | 5.7
$p$=2.0 | 50 | 3.9 | 3.8 | 3.4
| 100 | 5.7 | 6.2 | 6.2
| 200 | 6.0 | 6.7 | 7.1
$p$=3.0 | 50 | 5.5 | 6.3 | 6.1
| 100 | 5.1 | 4.6 | 4.4
| 200 | 5.2 | 4.7 | 5.7
Table 2: Power of the test against mixture alternatives of the form $0.7{\cal{NG}}(1,1,1)+0.3{\cal{NG}}(1,p,1)$ with sample size $n$ | $n$ |  $\gamma$=4 |  $\gamma$=6 |  $\gamma$=8
---|---|---|---|---
$p=0.5$ | 50 | 12.2 | 14.3 | 14.9
| 100 | 17.9 | 19.8 | 20.0
| 200 | 38.1 | 46.8 | 49.2
$p=2.0$ | 50 | 33.0 | 35.1 | 35.3
| 100 | 50.5 | 52.5 | 52.9
| 200 | 62.3 | 65.4 | 68.1
$p=3.0$ | 50 | 55.3 | 54.9 | 55.3
| 100 | 61.2 | 62.1 | 62.7
| 200 | 74.4 | 75.7 | 76.5

From Table 1 we observe that the empirical size of the test $T_{n,\gamma}$
varies with the value of the tuning parameter $\gamma$, but it is overall,
with a few exceptions, close to the $5\%$ nominal size of the test. Moreover,
looking at Table 2 we observe that the power of the test progressively
increases along with the sample size $n$. Also the value of the tuning
parameter $\gamma$ is seen to exert a certain effect on the power of the test,
which is however mostly mild.

##  5 Testing for the stable/gamma model

###  5.1 Test statistic

We now extend our test from the normal/gamma to the more general case of
testing the stable/gamma SFM whereby the two–sided error $v$ follows a
symmetric stable distribution, rather than a normal distribution. Specifically
we wish to test the null hypothesis

$\begin{split}{\Upsilon_{0}}:\;&\textrm{Model (\ref{SFM}) holds true with
}v\sim{\cal{S}}(\alpha,\kappa),\ u\sim{\cal{G}}(p,c),\textrm{ for some
}0<\alpha\leq 2,\textrm{and}\ \kappa,p,c>0,\end{split}$ (16)

where ${\cal{S}}(\alpha,\kappa)$, denotes a symmetric stable distribution with
tail–index $\alpha$ and scale parameter $\kappa>0$.

In this connection recall that the CF of a generic random variable $Z$ is
defined by $\varphi_{Z}(t)=\mathbb{E}({e^{{\rm{i}}tZ}}),\>t\in\mathbb{R}$,
${\rm{i}}=\sqrt{-1}$, and that the CF of a symmetric stable distribution is
given by $e^{-\kappa^{\alpha}|t|^{\alpha}}$666When we write stable
distribution, we actually mean the symmetric stable distribution. Note that
for $\alpha=2$, this distribution reduces to the normal distribution with
variance equal to $2\kappa^{2}$, which is the only stable distribution with
finite variance. Otherwise if $\alpha<2$, the stable distribution has finite
moments only of order strictly less than $\alpha$. Thus the mean of the
distribution is equal to zero provided that $\alpha>1$. Accordingly, and since
stable laws have potential for application mostly, if not exclusively, in the
finite–mean situation, hereafter we only consider the case of $\alpha>1$. For
further properties of this distribution we refer to Paolella (2007).

In order to introduce our test we note that under the null hypothesis of a
stable/gamma SFM, the CF of the standardized composed error term
$\widetilde{\varepsilon}=\varepsilon/c$ is given by

$\varphi_{\widetilde{\varepsilon}}(t)=\frac{e^{-\lambda|t|^{\alpha}}}{(1+{\rm{i}}t)^{p}},$

and consequently it satisfies the differential equation

$\Delta(t):=(1+{\rm{i}}t)\varphi^{\prime}_{\widetilde{\varepsilon}}(t)+\left[{\rm{i}}p+\alpha\lambda|t|^{\alpha-1}(1+{\rm{i}}t){\rm{sgn}}(t)\right]\varphi_{\widetilde{\varepsilon}}(t)=0,$

where $\lambda=\left(\kappa/c\right)^{\alpha}$, and ${\rm{sgn}}(t)=-1,0,1$ if
$t<0$, $t=0$, and $t>0$, respectively; see also Meintanis (2005). Thus arguing
as before we suggest the test statistic

${\cal{T}}_{n,w}=n\int_{-\infty}^{\infty}\left|\Delta_{n}(t)\right|^{2}\;w(t)\;{\rm{d}}t,$
(17)

where

$\Delta_{n}(t)=(1+{\rm{i}}t)\varphi^{\prime}_{n}(t)+\left[{\rm{i}}\widehat{p}+\widehat{\alpha}\widehat{\lambda}|t|^{\widehat{\alpha}-1}(1+{\rm{i}}t){\rm{sgn}}(t)\right]\varphi_{n}(t),$

with

$\varphi_{n}(t)=\frac{1}{n}\sum_{j=1}^{n}e^{{\rm{i}}tr_{j}}:=C_{n}(t)+{\rm{i}}S_{n}(t),$
$\varphi^{\prime}_{n}(t)=\frac{1}{n}\sum_{j=1}^{n}{\rm{i}}\>r_{j}e^{{\rm{i}}tr_{j}}:=C^{\prime}_{n}(t)+{\rm{i}}S^{\prime}_{n}(t),$

being the empirical CF and its derivative, respectively, where

$C_{n}(t)=\frac{1}{n}\sum_{j=1}^{n}\cos tr_{j},\
S_{n}(t)=\frac{1}{n}\sum_{j=1}^{n}\sin tr_{j},$

denote the real and imaginary part, respectively, of $\varphi_{n}(t)$, and
$C^{\prime}_{n}(t)$ and $S^{\prime}_{n}(t)$ the corresponding parts of
$\varphi_{n}^{\prime}(t)$. Also we write $|z|=\sqrt{z^{2}_{1}+z^{2}_{2}}$, for
the modulus of a complex number $z=z_{1}+{\rm{i}}z_{2}$. The consistency of
the test that rejects the null hypothesis of a stable/gamma SFM for large
values of the test statistic ${\cal{T}}_{n,w}$ figuring in (17) may be shown
by analogous arguments as in the case of the test in equation (9); see also
Meintanis (2005).

Following straightforward computations we obtain

$\begin{split}|\Delta_{n}(t)|^{2}&=(1+t^{2})|\varphi_{n}^{\prime}(t)|^{2}+\left[(\widehat{\alpha}\widehat{\lambda}|t|^{\widehat{\alpha}-1})^{2}+(\widehat{\alpha}\widehat{\lambda}|t|^{\widehat{\alpha}}+\widehat{p})^{2}\right]|\varphi_{n}(t)|^{2}+2\widehat{p}\left(S^{\prime}_{n}(t)C_{n}(t)-C^{\prime}_{n}(t)S_{n}(t)\right)\\\
&\hskip
10.00002pt+2\left[\widehat{\alpha}\widehat{\lambda}|t|^{\widehat{\alpha}-1}{\rm{sgn}}(t)+t(\widehat{p}+\widehat{\alpha}\widehat{\lambda}|t|^{\widehat{\alpha}})\right](C_{n}^{\prime}(t)C_{n}(t)+S_{n}^{\prime}(t)S_{n}(t))\end{split}$
(18)

Now if we consider $w(t)=e^{-\gamma t^{2}}$ as weight function, then the test
statistic may be explicitly written by making use of the integrals

$I_{\nu,\gamma}(z)=\int_{0}^{\infty}t^{\nu}\cos(tz)e^{-\gamma
t^{2}}{\rm{d}}t,\ \nu>-2,$

and

$J_{\nu,\gamma}(z)=\int_{0}^{\infty}t^{\nu}\sin(tz)e^{-\gamma
t^{2}}{\rm{d}}t,\ \nu>-1.$

Specifically following a long but otherwise straightforward computation we
have from (17) and (18),

$\begin{split}{\cal{T}}_{n,\gamma}&=\frac{1}{n}\sum_{j,k=1}^{n}\Big{\\{}2r_{j}r_{k}I_{2,\gamma}(r_{j}-r_{k})+2\big{(}r_{j}r_{k}+2\widehat{p}r_{j}+\widehat{p}^{2}\big{)}I_{0,\gamma}(r_{j}-r_{k})+4\widehat{p}\widehat{\alpha}\widehat{\lambda}I_{\widehat{\alpha},\gamma}(r_{j}-r_{k})+4\widehat{p}r_{k}J_{1,\gamma}(r_{j}-r_{k})\\\
&\hskip
70.0001pt+2\widehat{\alpha}^{2}\widehat{\lambda}^{2}\Big{(}I_{2\widehat{\alpha},\gamma}(r_{j}-r_{k})+I_{2(\widehat{\alpha}-1),\gamma}(r_{j}-r_{k})\Big{)}+4\widehat{\alpha}\widehat{\lambda}r_{j}\Big{(}J_{\widehat{\alpha}+1,\gamma}(r_{j}-r_{k})+J_{\widehat{\alpha}-1,\gamma}(r_{j}-r_{k})\Big{\\}}\end{split}$
(19)

Further details on the above computational formula for ${\cal{T}}_{n,\gamma}$
and for the integrals $I_{\nu,\gamma}(\cdot),J_{\nu,\gamma}(\cdot)$, are
provided in Appendices A.4 and A.5, respectively.

###  5.2 Parameter estimation

Following Tsionas (2012) we suggest CF–based estimators for the parameter
vector $\Theta:=(\kappa,\alpha,p,c)$. Note in this connection that the CF of
the stable/gamma error term $\varepsilon$ is given by
$\varphi_{\varepsilon}(t)=e^{-\kappa^{\alpha}|t|^{\alpha}}/(1+{\rm{i}}ct)^{p}$.
Thus using the OLS estimates for $\widehat{\beta}$, then under the null
hypothesis of a stable/gamma SFM, our method consists of estimating $\Theta$
by fitting an empirical version of this CF to its population counterpart, in
the logarithmic scale777the principal value of the complex logarithm is used.
Specifically we define the estimator as

$\widehat{\Theta}:=\underset{\Theta}{\mathrm{argmin}}\Big{\\{}\sum_{t\in\mathbb{T}}\left|-\kappa^{\alpha}|t|^{\alpha}-p\log(1+{\rm{i}}ct)-\log\widetilde{\varphi}_{n}(t)\right|^{2}\Big{\\}},$

where
$\widehat{\Theta}:=(\widehat{\kappa},\widehat{\alpha},\widehat{p},\widehat{c})$,
$\widetilde{\varphi}_{n}(t)=n^{-1}\sum_{j=1}^{n}e^{{\rm{i}}t\frac{\widehat{\varepsilon}_{j}}{{\rm{mad}}(\widehat{\varepsilon})}}$,
with ${\rm{mad}}(\widehat{\varepsilon})$ being the mean absolute deviation of
the OLS residuals $(\widehat{\varepsilon}_{j},\ j=1,...,n)$, and $\mathbb{T}$
is a finite grid of points in the neighborhood of zero. In particular we set
${\mathbb{T}}=\\{\pm 0.01,\pm 0.005\\}$.

##  6 Data and empirical results

###  6.1 Data

In order to examine the behaviour of the proposed test statistic defined in
equation (19) we use four datasets from previous studies – which are also
employed by Horrace and Wang (2022) – whereby we consider the stable/gamma
distribution as the null hypothesis of the composed error term. The datasets
are in Christensen and Greene (1976) on US electric utility companies, in
Tsionas (2006) and Feng and Serletis (2009) on US bank costs (considering a
cost function), and in Polachek and Yoon (1996) on the earnings of workers
(considering a production function). Note that the above mentioned studies
assume normal/half-normal or normal/exponential distribution for the composed
error component $\varepsilon$. These data are publicly available in the
Journal of Applied Econometrics Data
Archive.888http://qed.econ.queensu.ca/jae/

In order to apply the test we consider the parametric bootstrap version of the
test with $B=100$ bootstrap samples from the hypothesized stable/gamma law
with estimated parameters obtained as described in Subsection 5.2. (We refer
to Cornea–Madeira and Davidson (2015) for the validity of the parametric
bootstrap in the case of heavy–tailed distributions such as the stable
distributions). The nominal level of the test is set to $5\%$ and use
$\gamma=1$ as value of the tuning parameter.999We also tried $\gamma=0.5$ and
$\gamma=2$, leading to the same conclusions

For consistency and comparison reasons we follow Tsionas (2002) and take the
negatives of $Y$ and $X$ in the SFM of equation (1).

###  6.2 Empirical results

First we turn our attention to a widely used dataset in stochastic frontier
context, that of the cross-sectional case of 123 US electric utility companies
in 1970 used by Christensen and Greene (1976) and Tsionas (2002). We consider
the following functional form for the cost function

$\displaystyle-\log{\frac{C_{j}}{PF_{j}}}$ $\displaystyle=$
$\displaystyle-\beta_{0j}-\beta_{1j}\log{Q_{j}}-\beta_{2j}\log^{2}{Q_{j}}-\beta_{3j}\log{\frac{PL_{j}}{PF_{j}}}-\beta_{4j}\log{\frac{PK_{j}}{PF_{j}}}+v_{j}+u_{j}$

where $Q_{j}$ and $C_{j}$ are the output and the cost of the $j^{\rm{th}}$
firm respectively, and $PL,PK,PF$ denote the prices of the three production
factors, namely labour, capital, and fuel.

The parameter estimates, as well as the p-value of the test statistic are
presented in Table 3. The slope and intercept estimates are very close to
those of Table I in Tsionas (2002), while the p-value of the test statistic is
clearly non–significant for a stable/gamma null hypothesis $\Upsilon_{0}$. At
the same time though, the estimate of the tail–index $\alpha$ borders the case
of normality $\alpha=2$, and thus our findings may be interpreted as being in
agreement with those presented by Horrace and Wang (2022) implying a
medium–tailed normal/gamma SFM.

Table 3: Parameter estimates and p-value of the test for the dataset in Christensen and Green (1976) for the stable/gamma null hypothesis. Sample size $n=123$. Critical values obtained by $B=100$ bootstrap samples. $\widehat{\beta}_{0}$ | -7.205 | $\widehat{\kappa}$ | 1.1249
---|---|---|---
$\widehat{\beta}_{1}$ | 0.386 | $\widehat{\alpha}$ | 1.9999
$\widehat{\beta}_{2}$ | 0.032 | $\widehat{p}$ | 0.2223
$\widehat{\beta}_{3}$ | 0.246 | $\widehat{c}$ | 0.1942
$\widehat{\beta}_{4}$ | 0.079 |  p–value=0.49

We next revisit the panel data set in Polachek and Yoon (1996), the Panel
Study of Income Dynamics (PSID), with the following specification of the
production function

$\displaystyle\log W_{jt}$ $\displaystyle=$
$\displaystyle\beta_{0j}+\beta_{1j}EDU_{jt}+\beta_{2j}EX_{jt}+\beta_{3j}EX_{jt}^{2}+\beta_{4j}TEN_{jt}+\beta_{5j}TEN_{jt}^{2}+v_{jt}-u_{jt}$

where $W_{jt},EDU_{jt},EX_{jt},TEN_{jt}$ represent the wage, education and
experience in years and tenure in months, respectively, for $838$ workers in
the years 1969–1984.

The corresponding results for the stable/gamma null hypothesis are presented
in Table 4. These results are clearly significant for a number of years
(namely for the years 1970, 1980, 1982), while in other cases lead to non-
rejection, but the p-value is either borderline or the estimate of the
tail–index is too close to normality (namely for the years 1969, 1972, 1973,
1974, 1976, 1977, 1983, and 1984). However there exist also (sort of speak)
“clear–cut” cases for the presence of heavy tails, such as those corresponding
to the years 1971, 1975, 1978, 1979, and 1981.101010Note in this connection
that the tails of the stable distribution are quite sensitive to the value of
the tail index and that a value, say, $\alpha=1.98$ implies the presence of a
genuinely heavy–tailed distribution.

In the third real–data example, we employ the dataset of large111111The
commercial banks included in the Report of Condition and Income (Call Report)
whose total assets exceeded one billion dollars for at least three (not
necessarily consecutive) years. US banks’ costs used in Tsionas (2006). The
cost function is given as follows

$-\log C_{jt}=X_{jt}\beta_{j}+v_{jt}+u_{jt},$

where $C_{jt}$ is the cost of the $j^{\rm{th}}$ bank in year $t$, and $X_{jt}$
contains a constant, the logarithms of the outputs (installment loans,
installment loans (to individuals for personal/household expenses), real
estate loans, business loans, federal funds sold and securities purchased
under agreements to resell, other assets (assets that cannot be properly
included in any other asset items in the balance sheet)) and the logarithms of
the relevant prices of the inputs (labour, capital, purchased funds, interest-
bearing deposits in total transaction accounts, interest-bearing deposits in
total non-transaction accounts); relevant prices calculated by dividing total
expenses of each input by the corresponding quantity. The results for the
$128$ banks for $12$ years (1989–2000) are presented in Table 5. The bootstrap
p–values are clearly non–significant for $\Upsilon_{0}$ with each year but
there exist a number of cases where the estimate of the tail–index borders
normality whereas, again, for several other years we can not reject a
stable/gamma SFM involving a sub–Gaussian stable law.

In our last application we apply the proposed test statistic on the subsample
XII in Feng and Serletis (2009) containing data of $1262$ banks for 8 years
(1998–2005). We consider the following functional form for the cost function

$-\ln{(C_{jt})}=-\beta_{0j}-\beta_{1j}\ln{PL_{jt}}-\beta_{2j}\ln{PF_{jt}}-\beta_{3j}\ln{PK_{jt}}-\beta_{4j}\ln{CL_{jt}}-\beta_{5j}\ln{NCL_{jt}}-\beta_{6j}\ln{SEC_{jt}}+v_{j}+u_{j},$

where $C_{jt}$ is the total cost of the $j^{\rm th}$ firm in year $t$, $PL$ is
the wage rate for labor, $PF$ is the interest rate for borrowed funds and $PK$
the price of physical capital, $CL$ denotes the consumer loans, $NCL$ the non
consumer loans, and $SC$ gives securities including all non-financial assets.
The corresponding results are reported in Table 6, leading to safely assume
(or not reject) a stable/gamma null hypothesis for the years 2000 and 2004.
For the remaining years the results are borderline for $\Upsilon_{0}$ either
in the p-value or in the estimate of the tail–index $\alpha$.

The conclusions drawn by the three last applications are in partial agreement
with those of Horrace and Wang (2022) obtained by an entirely non–parametric
method, with some discrepancies potentially be attributed to the lack of fit
of the particular heavy–tailed law under test. In this connection the stable
law advocated herein is one option while the Student–t distribution may also
be considered as a potential heavy tailed model; see for instance Wheat et al.
(2019). The two families however are generally separate families of
distributions having only the Cauchy law in common, a law that is extremely
heavy–tailed which makes it an unlikely candidate for the distribution of the
pure random error component $v$.

Table 4: Parameter estimates and p-value of the test for the dataset in Polacheck and Yoon (1996) for the stable/gamma null hypothesis. Sample size $n=838$. Critical values obtained by $B=100$ bootstrap samples.  1969 | $\hat{\beta}_{0}$ | -1.238 | $\hat{\kappa}$ | 1.1260 | 1977 | $\hat{\beta}_{0}$ | -0.743 | $\hat{\kappa}$ | 1.1415
---|---|---|---|---|---|---|---|---|---
$\hat{\beta}_{1}$ | 0.066 | $\hat{\alpha}$ | 1.9999 | $\hat{\beta}_{1}$ | 0.049 | $\hat{\alpha}$ | 1.9995
$\hat{\beta}_{2}$ | 0.043 | $\hat{p}$ | 0.2222 | $\hat{\beta}_{2}$ | 0.025 | $\hat{p}$ | 0.2132
$\hat{\beta}_{3}$ | -0.0008 | $\hat{c}$ | 0.1941 | $\hat{\beta}_{3}$ | -0.0004 | $\hat{c}$ | 0.1939
$\hat{\beta}_{4}$ | 0.005 |  p–value=0.17 | $\hat{\beta}_{4}$ | 0.001 |  p–value=0.17
$\hat{\beta}_{5}$ | -3.3E-5 |  |  | $\hat{\beta}_{5}$ | -2.2E-6 |  |
1970 | $\hat{\beta}_{0}$ | -1.401 | $\hat{\kappa}$ | 0.9924 | 1978 | $\hat{\beta}_{0}$ | -0.683 | $\hat{\kappa}$ | 1.0834
$\hat{\beta}_{1}$ | 0.061 | $\hat{\alpha}$ | 1.8893 | $\hat{\beta}_{1}$ | 0.048 | $\hat{\alpha}$ | 1.9893
$\hat{\beta}_{2}$ | 0.060 | $\hat{p}$ | 0.1905 | $\hat{\beta}_{2}$ | 0.025 | $\hat{p}$ | 0.2292
$\hat{\beta}_{3}$ | -0.0012 | $\hat{c}$ | 0.1923 | $\hat{\beta}_{3}$ | -0.0004 | $\hat{c}$ | 0.1819
$\hat{\beta}_{4}$ | 0.007 |  p–value=0.01 | $\hat{\beta}_{4}$ | 0.001 |  p–value=0.15
$\hat{\beta}_{5}$ | -4.7E-5 |  |  | $\hat{\beta}_{5}$ | -3.3E-6 |  |
1971 | $\hat{\beta}_{0}$ | -1.050 | $\hat{\kappa}$ | 0.8845 | 1979 | $\hat{\beta}_{0}$ | -0.722 | $\hat{\kappa}$ | 0.9807
$\hat{\beta}_{1}$ | 0.067 | $\hat{\alpha}$ | 1.9840 | $\hat{\beta}_{1}$ | 0.055 | $\hat{\alpha}$ | 1.9674
$\hat{\beta}_{2}$ | 0.038 | $\hat{p}$ | 0.1700 | $\hat{\beta}_{2}$ | 0.018 | $\hat{p}$ | 0.1843
$\hat{\beta}_{3}$ | -0.0007 | $\hat{c}$ | 0.1893 | $\hat{\beta}_{3}$ | -0.0003 | $\hat{c}$ | 0.1885
$\hat{\beta}_{4}$ | -0.007 |  p–value=0.41 | $\hat{\beta}_{4}$ | 0.001 |  p–value=0.12
$\hat{\beta}_{5}$ | 5.2E-5 |  |  | $\hat{\beta}_{5}$ | -3.6E-6 |  |
1972 | $\hat{\beta}_{0}$ | -1.135 | $\hat{\kappa}$ | 1.0503 | 1980 | $\hat{\beta}_{0}$ | -0.630 | $\hat{\kappa}$ | 1.0848
$\hat{\beta}_{1}$ | 0.052 | $\hat{\alpha}$ | 1.9646 | $\hat{\beta}_{1}$ | 0.046 | $\hat{\alpha}$ | 1.8738
$\hat{\beta}_{2}$ | 0.046 | $\hat{p}$ | 0.1871 | $\hat{\beta}_{2}$ | 0.020 | $\hat{p}$ | 0.1871
$\hat{\beta}_{3}$ | -0.0008 | $\hat{c}$ | 0.1786 | $\hat{\beta}_{3}$ | -0.0003 | $\hat{c}$ | 0.2137
$\hat{\beta}_{4}$ | 0.006 |  p–value=0.07 | $\hat{\beta}_{4}$ | 0.001 |  p–value=0.00
$\hat{\beta}_{5}$ | -3.4E-5 |  |  | $\hat{\beta}_{5}$ | -2.1E-6 |  |
1973 | $\hat{\beta}_{0}$ | -0.996 | $\hat{\kappa}$ | 0.9685 | 1981 | $\hat{\beta}_{0}$ | -0.860 | $\hat{\kappa}$ | 0.9282
$\hat{\beta}_{1}$ | 0.058 | $\hat{\alpha}$ | 1.9938 | $\hat{\beta}_{1}$ | 0.055 | $\hat{\alpha}$ | 1.9882
$\hat{\beta}_{2}$ | 0.027 | $\hat{p}$ | 0.1921 | $\hat{\beta}_{2}$ | 0.029 | $\hat{p}$ | 0.1921
$\hat{\beta}_{3}$ | -0.0004 | $\hat{c}$ | 0.1971 | $\hat{\beta}_{3}$ | -0.0005 | $\hat{c}$ | 0.1950
$\hat{\beta}_{4}$ | 0.006 |  p–value=0.35 | $\hat{\beta}_{4}$ | 0.001 |  p–value=0.46
$\hat{\beta}_{5}$ | -3.3E-5 |  |  | $\hat{\beta}_{5}$ | -1.2E-6 |  |
1974 | $\hat{\beta}_{0}$ | -1.028 | $\hat{\kappa}$ | 0.9639 | 1982 | $\hat{\beta}_{0}$ | -0.703 | $\hat{\kappa}$ | 0.8860
$\hat{\beta}_{1}$ | 0.056 | $\hat{\alpha}$ | 1.9982 | $\hat{\beta}_{1}$ | 0.049 | $\hat{\alpha}$ | 1.8185
$\hat{\beta}_{2}$ | 0.037 | $\hat{p}$ | 0.1826 | $\hat{\beta}_{2}$ | 0.018 | $\hat{p}$ | 0.1847
$\hat{\beta}_{3}$ | -0.0007 | $\hat{c}$ | 0.1742 | $\hat{\beta}_{3}$ | -0.0003 | $\hat{c}$ | 0.2095
$\hat{\beta}_{4}$ | 0.005 |  p–value=0.61 | $\hat{\beta}_{4}$ | 0.001 |  p–value=0.00
$\hat{\beta}_{5}$ | -3.0E-5 |  |  | $\hat{\beta}_{5}$ | -3.2E-6 |  |
1975 | $\hat{\beta}_{0}$ | -0.715 | $\hat{\kappa}$ | 1.0294 | 1983 | $\hat{\beta}_{0}$ | -0.882 | $\hat{\kappa}$ | 1.0596
$\hat{\beta}_{1}$ | 0.048 | $\hat{\alpha}$ | 1.9845 | $\hat{\beta}_{1}$ | 0.049 | $\hat{\alpha}$ | 1.9748
$\hat{\beta}_{2}$ | 0.023 | $\hat{p}$ | 0.1932 | $\hat{\beta}_{2}$ | 0.031 | $\hat{p}$ | 0.2154
$\hat{\beta}_{3}$ | -0.0004 | $\hat{c}$ | 0.2154 | $\hat{\beta}_{3}$ | -0.0005 | $\hat{c}$ | 0.1842
$\hat{\beta}_{4}$ | 0.021 |  p–value=0.15 | $\hat{\beta}_{4}$ | 0.002 |  p–value=0.06
$\hat{\beta}_{5}$ | -5.4E-4 |  |  | $\hat{\beta}_{5}$ | -3.5E-6 |  |
1976 | $\hat{\beta}_{0}$ | -0.0869 | $\hat{\kappa}$ | 1.1047 | 1984 | $\hat{\beta}_{0}$ | -0.638 | $\hat{\kappa}$ | 1.0571
$\hat{\beta}_{1}$ | 0.051 | $\hat{\alpha}$ | 1.9976 | $\hat{\beta}_{1}$ | 0.049 | $\hat{\alpha}$ | 1.9972
$\hat{\beta}_{2}$ | 0.029 | $\hat{p}$ | 0.1918 | $\hat{\beta}_{2}$ | 0.014 | $\hat{p}$ | 0.2190
$\hat{\beta}_{3}$ | -0.0005 | $\hat{c}$ | 0.2309 | $\hat{\beta}_{3}$ | -0.0002 | $\hat{c}$ | 0.1843
$\hat{\beta}_{4}$ | 0.003 |  p–value=0.22 | $\hat{\beta}_{4}$ | 0.001 |  p–value=0.30
$\hat{\beta}_{5}$ | -1.6E-5 |  |  | $\hat{\beta}_{5}$ | -2.2E-6 |  |
Table 5: Parameter estimates and p-value of the test for the dataset in Tsionas (2006) for the stable/gamma null hypothesis. Sample size $n=128$. Critical values obtained by $B=100$ bootstrap samples. 1989 | $\hat{\beta}_{0}$ | 3.078 | $\hat{\kappa}$ | 1.1251 | 1995 | $\hat{\beta}_{0}$ | 2.671 | $\hat{\kappa}$ | 1.0291
---|---|---|---|---|---|---|---|---|---
$\hat{\beta}_{1}$ | 0.092 | $\hat{\alpha}$ | 1.9999 | $\hat{\beta}_{1}$ | 0.237 | $\hat{\alpha}$ | 1.9846
$\hat{\beta}_{2}$ | -0.006 | $\hat{p}$ | 0.2223 | $\hat{\beta}_{2}$ | -0.020 | $\hat{p}$ | 0.1932
$\hat{\beta}_{3}$ | 0.181 | $\hat{c}$ | 0.1942 | $\hat{\beta}_{3}$ | 0.172 | $\hat{c}$ | 0.2154
$\hat{\beta}_{4}$ | 0.218 |  p–value=0.83 | $\hat{\beta}_{4}$ | 0.041 |  p–value=0.79
$\hat{\beta}_{5}$ | 0.008 |  |  | $\hat{\beta}_{5}$ | 0.015 |  |
$\hat{\beta}_{6}$ | 0.260 |  |  | $\hat{\beta}_{6}$ | 0.242 |  |
$\hat{\beta}_{7}$ | 0.121 |  |  | $\hat{\beta}_{7}$ | 0.057 |  |
$\hat{\beta}_{8}$ | 0.059 |  |  | $\hat{\beta}_{8}$ | 0.079 |  |
$\hat{\beta}_{9}$ | 0.497 |  |  | $\hat{\beta}_{9}$ | 0.527 |  |
1990 | $\hat{\beta}_{0}$ | 3.438 | $\hat{\kappa}$ | 0.9909 | 1996 | $\hat{\beta}_{0}$ | 0.938 | $\hat{\kappa}$ | 1.1040
$\hat{\beta}_{1}$ | -0.042 | $\hat{\alpha}$ | 1.8924 | $\hat{\beta}_{1}$ | 0.529 | $\hat{\alpha}$ | 1.9977
$\hat{\beta}_{2}$ | 0.038 | $\hat{p}$ | 0.1907 | $\hat{\beta}_{2}$ | 0.047 | $\hat{p}$ | 0.1919
$\hat{\beta}_{3}$ | -0.017 | $\hat{c}$ | 0.1925 | $\hat{\beta}_{3}$ | 0.220 | $\hat{c}$ | 0.2311
$\hat{\beta}_{4}$ | 0.309 |  p–value=0.43 | $\hat{\beta}_{4}$ | 0.072 |  p–value=0.78
$\hat{\beta}_{5}$ | 0.024 |  |  | $\hat{\beta}_{5}$ | 0.047 |  |
$\hat{\beta}_{6}$ | 0.300 |  |  | $\hat{\beta}_{6}$ | 0.288 |  |
$\hat{\beta}_{7}$ | 0.129 |  |  | $\hat{\beta}_{7}$ | 0.035 |  |
$\hat{\beta}_{8}$ | 0.056 |  |  | $\hat{\beta}_{8}$ | 0.069 |  |
$\hat{\beta}_{9}$ | 0.463 |  |  | $\hat{\beta}_{9}$ | 0.435 |  |
1991 | $\hat{\beta}_{0}$ | 1.955 | $\hat{\kappa}$ | 0.8843 | 1997 | $\hat{\beta}_{0}$ | 0.828 | $\hat{\kappa}$ | 1.1408
$\hat{\beta}_{1}$ | 0.260 | $\hat{\alpha}$ | 1.9841 | $\hat{\beta}_{1}$ | 0.489 | $\hat{\alpha}$ | 1.9995
$\hat{\beta}_{2}$ | 0.046 | $\hat{p}$ | 0.1701 | $\hat{\beta}_{2}$ | -0.063 | $\hat{p}$ | 0.2133
$\hat{\beta}_{3}$ | 0.115 | $\hat{c}$ | 0.1893 | $\hat{\beta}_{3}$ | 0.149 | $\hat{c}$ | 0.1940
$\hat{\beta}_{4}$ | 0.366 |  p–value=0.94 | $\hat{\beta}_{4}$ | 0.034 |  p–value=0.74
$\hat{\beta}_{5}$ | 0.061 |  |  | $\hat{\beta}_{5}$ | 0.009 |  |
$\hat{\beta}_{6}$ | 0.277 |  |  | $\hat{\beta}_{6}$ | 0.316 |  |
$\hat{\beta}_{7}$ | 0.080 |  |  | $\hat{\beta}_{7}$ | 0.064 |  |
$\hat{\beta}_{8}$ | 0.042 |  |  | $\hat{\beta}_{8}$ | 0.054 |  |
$\hat{\beta}_{9}$ | 0.497 |  |  | $\hat{\beta}_{9}$ | 0.473 |  |
1992 | $\hat{\beta}_{0}$ | 2.782 | $\hat{\kappa}$ | 1.0492 | 1998 | $\hat{\beta}_{0}$ | 2.627 | $\hat{\kappa}$ | 1.0830
$\hat{\beta}_{1}$ | 0.239 | $\hat{\alpha}$ | 1.9653 | $\hat{\beta}_{1}$ | 0.334 | $\hat{\alpha}$ | 1.9894
$\hat{\beta}_{2}$ | 0.032 | $\hat{p}$ | 0.1874 | $\hat{\beta}_{2}$ | 0.039 | $\hat{p}$ | 0.2294
$\hat{\beta}_{3}$ | 0.130 | $\hat{c}$ | 0.1789 | $\hat{\beta}_{3}$ | 0.284 | $\hat{c}$ | 0.1820
$\hat{\beta}_{4}$ | 0.278 |  p–value=0.70 | $\hat{\beta}_{4}$ | 0.131 |  p–value=0.72
$\hat{\beta}_{5}$ | 0.011 |  |  | $\hat{\beta}_{5}$ | 0.006 |  |
$\hat{\beta}_{6}$ | 0.255 |  |  | $\hat{\beta}_{6}$ | 0.358 |  |
$\hat{\beta}_{7}$ | 0.061 |  |  | $\hat{\beta}_{7}$ | 0.080 |  |
$\hat{\beta}_{8}$ | 0.035 |  |  | $\hat{\beta}_{8}$ | 0.015 |  |
$\hat{\beta}_{9}$ | 0.549 |  |  | $\hat{\beta}_{9}$ | 0.358 |  |
1993 | $\hat{\beta}_{0}$ | 3.268 | $\hat{\kappa}$ | 0.9680 | 1999 | $\hat{\beta}_{0}$ | 2.356 | $\hat{\kappa}$ | 0.9800
$\hat{\beta}_{1}$ | 0.150 | $\hat{\alpha}$ | 1.9883 | $\hat{\beta}_{1}$ | 0.288 | $\hat{\alpha}$ | 1.9678
$\hat{\beta}_{2}$ | 0.032 | $\hat{p}$ | 0.1921 | $\hat{\beta}_{2}$ | 0.046 | $\hat{p}$ | 0.1844
$\hat{\beta}_{3}$ | 0.298 | $\hat{c}$ | 0.1972 | $\hat{\beta}_{3}$ | 0.287 | $\hat{c}$ | 0.1886
$\hat{\beta}_{4}$ | 0.248 |  p–value=0.89 | $\hat{\beta}_{4}$ | 0.092 |  p–value=0.74
$\hat{\beta}_{5}$ | 0.018 |  |  | $\hat{\beta}_{5}$ | 0.001 |  |
$\hat{\beta}_{6}$ | 0.253 |  |  | $\hat{\beta}_{6}$ | 0.336 |  |
$\hat{\beta}_{7}$ | 0.081 |  |  | $\hat{\beta}_{7}$ | 0.087 |  |
$\hat{\beta}_{8}$ | 0.059 |  |  | $\hat{\beta}_{8}$ | 0.029 |  |
$\hat{\beta}_{9}$ | 0.515 |  |  | $\hat{\beta}_{9}$ | 0.417 |  |
1994 | $\hat{\beta}_{0}$ | 2.841 | $\hat{\kappa}$ | 0.9637 | 2000 | $\hat{\beta}_{0}$ | 1.088 | $\hat{\kappa}$ | 1.0847
$\hat{\beta}_{1}$ | 0.194 | $\hat{\alpha}$ | 1.9982 | $\hat{\beta}_{1}$ | 0.435 | $\hat{\alpha}$ | 1.8740
$\hat{\beta}_{2}$ | 0.074 | $\hat{p}$ | 0.1826 | $\hat{\beta}_{2}$ | -0.048 | $\hat{p}$ | 0.1871
$\hat{\beta}_{3}$ | 0.124 | $\hat{c}$ | 0.1743 | $\hat{\beta}_{3}$ | 0.320 | $\hat{c}$ | 0.2137
$\hat{\beta}_{4}$ | 0.281 |  p–value=0.95 | $\hat{\beta}_{4}$ | 0.064 |  p–value=0.30
$\hat{\beta}_{5}$ | 0.034 |  |  | $\hat{\beta}_{5}$ | -0.009 |  |
$\hat{\beta}_{6}$ | 0.254 |  |  | $\hat{\beta}_{6}$ | 0.374 |  |
$\hat{\beta}_{7}$ | 0.059 |  |  | $\hat{\beta}_{7}$ | 0.128 |  |
$\hat{\beta}_{8}$ | 0.015 |  |  | $\hat{\beta}_{8}$ | 0.030 |  |
$\hat{\beta}_{9}$ | 0.562 |  |  | $\hat{\beta}_{9}$ | 0.364 |  |
Table 6: Parameter estimates and p-value of the test for the dataset in Feng and Serletis (2009) for the stable/gamma null hypothesis. Sample size $n=1262$. Critical values obtained by $B=100$ bootstrap samples. 1998 | $\hat{\beta}_{0}$ | 0.149 | $\hat{\kappa}$ | 1.1253 | 2002 | $\hat{\beta}_{0}$ | -0.237 | $\hat{\kappa}$ | 0.9684
---|---|---|---|---|---|---|---|---|---
$\hat{\beta}_{1}$ | 0.087 | $\hat{\alpha}$ | 1.9999 | $\hat{\beta}_{1}$ | 0.182 | $\hat{\alpha}$ | 1.9982
$\hat{\beta}_{2}$ | 0.287 | $\hat{p}$ | 0.2222 | $\hat{\beta}_{2}$ | 0.295 | $\hat{p}$ | 0.1826
$\hat{\beta}_{3}$ | 0.024 | $\hat{c}$ | 0.1941 | $\hat{\beta}_{3}$ | 0.021 | $\hat{c}$ | 0.1972
$\hat{\beta}_{4}$ | 0.067 |  p–value=0.10 | $\hat{\beta}_{4}$ | 0.078 |  p–value=0.24
$\hat{\beta}_{5}$ | 0.523 |  |  | $\hat{\beta}_{5}$ | 0.512 |  |
$\hat{\beta}_{6}$ | 0.262 |  |  | $\hat{\beta}_{6}$ | 0.270 |  |
1999 | $\hat{\beta}_{0}$ | -0.011 | $\hat{\kappa}$ | 0.9914 | 2003 | $\hat{\beta}_{0}$ | -0.014 | $\hat{\kappa}$ | 0.9639
$\hat{\beta}_{1}$ | 0.254 | $\hat{\alpha}$ | 1.8915 | $\hat{\beta}_{1}$ | 0.185 | $\hat{\alpha}$ | 1.9982
$\hat{\beta}_{2}$ | 0.302 | $\hat{p}$ | 0.1907 | $\hat{\beta}_{2}$ | 0.261 | $\hat{p}$ | 0.1826
$\hat{\beta}_{3}$ | 0.030 | $\hat{c}$ | 0.1925 | $\hat{\beta}_{3}$ | 0.042 | $\hat{c}$ | 0.1742
$\hat{\beta}_{4}$ | 0.092 |  p–value=0.00 | $\hat{\beta}_{4}$ | 0.098 |  p–value=0.40
$\hat{\beta}_{5}$ | 0.490 |  |  | $\hat{\beta}_{5}$ | 0.485 |  |
$\hat{\beta}_{6}$ | 0.248 |  |  | $\hat{\beta}_{6}$ | 0.252 |  |
2000 | $\hat{\beta}_{0}$ | 0.065 | $\hat{\kappa}$ | 0.8844 | 2004 | $\hat{\beta}_{0}$ | -0.058 | $\hat{\kappa}$ | 1.0295
$\hat{\beta}_{1}$ | 0.112 | $\hat{\alpha}$ | 1.9841 | $\hat{\beta}_{1}$ | 0.135 | $\hat{\alpha}$ | 1.9845
$\hat{\beta}_{2}$ | 0.289 | $\hat{p}$ | 0.1700 | $\hat{\beta}_{2}$ | 0.313 | $\hat{p}$ | 0.1932
$\hat{\beta}_{3}$ | 0.029 | $\hat{c}$ | 0.1893 | $\hat{\beta}_{3}$ | 0.028 | $\hat{c}$ | 0.2154
$\hat{\beta}_{4}$ | 0.072 |  p–value=0.25 | $\hat{\beta}_{4}$ | 0.078 |  p–value=0.15
$\hat{\beta}_{5}$ | 0.520 |  |  | $\hat{\beta}_{5}$ | 0.524 |  |
$\hat{\beta}_{6}$ | 0.262 |  |  | $\hat{\beta}_{6}$ | 0.264 |  |
2001 | $\hat{\beta}_{0}$ | -0.311 | $\hat{\kappa}$ | 1.0496 | 2005 | $\hat{\beta}_{0}$ | 0.028 | $\hat{\kappa}$ | 1.1042
$\hat{\beta}_{1}$ | 0.236 | $\hat{\alpha}$ | 1.9651 | $\hat{\beta}_{1}$ | 0.271 | $\hat{\alpha}$ | 1.9977
$\hat{\beta}_{2}$ | 0.281 | $\hat{p}$ | 0.1873 | $\hat{\beta}_{2}$ | 0.293 | $\hat{p}$ | 0.1919
$\hat{\beta}_{3}$ | 0.045 | $\hat{c}$ | 0.1789 | $\hat{\beta}_{3}$ | 0.041 | $\hat{c}$ | 0.2310
$\hat{\beta}_{4}$ | 0.099 |  p–value=0.05 | $\hat{\beta}_{4}$ | 0.104 |  p–value=0.11
$\hat{\beta}_{5}$ | 0.498 |  |  | $\hat{\beta}_{5}$ | 0.468 |  |
$\hat{\beta}_{6}$ | 0.254 |  |  | $\hat{\beta}_{6}$ | 0.240 |  |

##  7 Conclusions

We propose computationally convenient goodness–of–fit tests in SFMs based on
empirical versions of distributional transforms (the moment generating and the
characteristic function) of the law of the composed error term
$\varepsilon=v-u$. Conditionally on correct specification of the production
function (or cost function), the interpretation of test results depends on the
postulated pure random error term: If we are confident enough about the law of
$v$ then the test result has a direct and definite bearing on the law of
inefficiency component $u$, while otherwise rejection of the null hypothesis
may also be due to poor fit of the postulated law of $v$. In either case our
tests are omnibus, meaning that they are sensitive to arbitrary deviations
from the null hypothesis. In finite–samples, this property is investigated by
means of bootstrap versions of the test statistics. Moreover a series of
applications confirms earlier findings regarding the presence, for selective
cases, of heavy–tails in the empirical modelling of production.

## References

  * Bera and Mallick (2002) Bera AK, Mallick NC (2002) Information matrix tests for the composed error frontier model. In: Balakrishnan N(ed) Advances on methodological and applied aspects of probability and statistics. Gordon and Breach Science Publishers, London
  * Chang and Hall (2015) Chang J, Hall P (2015) Double–bootstrap methods that use a single double–bootstrap simulation. Biometrika 102:203–214
  * Chen and Wang (2012) Chen YT, Wang HJ (2012) Centered-Residuals-Based Moment Estimator and Test for Stochastic Frontier Models. Econom Rev 31:625–653
  * Cheng et al. (2022) Cheng MY, Wang S, Xia, L, Zhang, X (2022) Testing specification of distribution in stochastic frontier analysis. J Econometr https://doi.org/10.1016/j.econom.2022.03.002
  * Christensen and Greene (1976) Christensen LR, Greene WH (1976) Economies of scale in US electric power generation. J of Political Econ 84:655–676
  * Coelli (1995) Coelli T (1995) Estimators and hypothesis tests for a stochastic frontier function: A Monte Carlo analysis. J Product Anal 6:247–268
  * CD (2015) Cornea–Madeira A, Davidson R (2015) A parameteric bootstrap for heavy–tailed distributions. Econometr Theor 31: 449–470
  * EHM (2012) Henze N, Meintanis SG, Ebner B (2012) Goodness–of–fit tests for the Gamma distribution based on the empirical Laplace transform. Commun Statist–Theor Meth 41: 1543–1556
  * Feng and Serletis (2009) Feng G, Serletis A (2009) Efficiency and productivity of the US banking industry, 1998-–2005: Evidence from the Fourier cost function satisfying global regularity conditions. J of Applied Econometr 24:105–138
  * Feu (1989) Feuerverger A (1989) On the empirical saddlepoint approximation. Biometrika 76:457–464
  * Giacomini et al. (2013) Giacomini R, Politis DN, White H (2013) A warp-speed method for conducting Monte Carlo experiments involving bootstrap estimators. Econometr Theor 29:567–589
  * Gradshteyn and Ryzhik (1994) Gradshteyn IS, Ryzhik IM (1994) Tables of Integrals, Series, and Products 6th edn. Academic Press, New York
  * Horrace and Wang (2022) Horrace WC, Wang Y (2022) Nonparametric Tests of Tail Behavior in Stochastic Frontier Models. J of Applied Econom 37:537–562
  * Jiang (2010) Jiang J (2010) Large Sample Techniques for Statistics. Springer, New York
  * Kopp and Mullahy (1990) Kopp RJ, Mullahy J (1990) Moment-based estimation and testing of stochastic frontier models. J Econom 46:165–183
  * Kopp and Mullahy (1993) Kopp RJ, Mullahy J (1993) Least squares estimation of econometric frontier models: Consistent estimation and inference. Scand J Econom 95:125–132
  * Lee (1983) Lee LF (1983) A test for distributional assumptions for the stochastic frontier functions. J Econom 22:245–267
  * Mei (2005) Meintanis SG (2005) Consistent tests for symmetric stability with finite mean based on the empirical characteristic function. J Statist Plann Inferen 128: 373-380.
  * Meintanis and Papadimitriou (2022) Meintanis SG, Papadimitriou CK. (2022) Goodness–of–fit tests for stochastic frontier models based on the characteristic function. J. Product. Anal. 57: 285-296.
  * Paolella (2007) Paolella MS (2007) Intermediate Probability. A Computational Approach. John Wiley & Sons, New York
  * Polachek and Yoon (1996) Polacheck SW, Yoon BJ (1996) Panel estimates of a two‐-tiered earnings frontier. J of Applied Econometr 11:169–178
  * Schmidt and Lin (1984) Schmidt P, Lin T (1984) Simple tests of alternative specifications in stochastic frontier models. J Econom 24:349–361
  * Severini (2005) Severini TA (2005) Elements of Distribution Theory. (Vol. 17) Cambridge University Press, Cambridge
  * Tsionas (2002) Tsionas EG (2002) Stochastic frontier models with random coefficients. J of Applied Econom 17:127–147
  * Tsionas (2006) Tsionas EG (2006) Inference in dynamic stochastic frontier models. J of Applied Econom 21:669–676
  * Tsionas (2012) Tsionas EG (2012) Maximum likelihood estimation of stochastic frontier models by the Fourier transform. J Econometr 170:234–248
  * Wang et al. (2011) Wang WS, Amsler C, Schmidt P (2011) Goodness of fit tests in stochastic frontier models. J Product Anal 35:95–118
  * Wheat et al. (2019) Wheat P, Stead AD, Greene WH (2019) Robust stochastic frontier analysis: a Student’s $t$–half normal model with application to highway maintenance costs in England. J Product Anal 51:21–38

##  Appendix A.1 Consistency of the MGF–based test

Assume that the law of the composed error $\varepsilon$ has a finite and
differentiable MGF. Assume also that
$(\widehat{\sigma}^{2}_{v},\widehat{c},\widehat{p})\to(\sigma_{v}^{2},c,p)\in(0,\infty)^{3}$,
almost surely, i.e. that the parameter estimates attain positive and finite
almost sure limits even under alternatives. Then we have by the pointwise
consistency of the empirical MGF (see Feuerverger, (1989)) that,
$D^{2}_{n}(t)\to D^{2}(t)$, where
$D(t)=((1+t)/c)M^{\prime}_{\varepsilon}(t/c)+[p-\lambda
t(1+t)]M_{\varepsilon}(t/c)$. If we further assume that the weight function
$w>0$ is such that

$\int_{0}^{\infty}t^{4}e^{tx}w(t){\rm{d}}t<\infty,\ \forall x\in\mathbb{R},$

it then follows by Fatou’s lemma (see Jiang (2010), §A.2.3) and (9) that

$\liminf_{n\to\infty}\frac{T_{n,w}}{n}\geq\int_{0}^{\infty}D^{2}(t)w(t){\rm{d}}t:=\Delta_{w},$

whereby a simple change of variables

$\Delta_{w}=c^{-1}\int_{0}^{\infty}\widetilde{D}^{2}(t)w(t){\rm{d}}t,$

with

$\widetilde{D}(t)=(1+ct)M^{\prime}_{\varepsilon}(t)+[cp-\sigma^{2}_{v}t(1+ct)]M_{\varepsilon}(t).$

Clearly $\Delta_{w}>0$ unless the integrand above vanishes identically in $t$.
However since $w>0$ and due to the uniqueness of solution of the differential
equation in (4) (with initial condition $M_{\varepsilon}(0)=1$), this
integrand vanishes only if the MGF of $\varepsilon$ is given by (3). In turn,
by the uniqueness of MGFs the last implication is equivalent to the null
hypothesis ${\cal{H}}_{0}$ figuring in (2), i.e. to a normal/gamma SFM with
the indicated parameters. Otherwise under alternatives, $T_{n,w}\to\infty$,
a.s. as $n\to\infty$, and thus the test which rejects ${\cal{H}}_{0}$ for
large values of $T_{n,w}$ is consistent.

######  Remark 1.

Consistency of the MGF–based test critically depends on the one–to–one
correspondence between distributions and MGFs. However strictly speaking this
uniqueness assumes that the underlying MGF is defined in a neighborhood of
zero (see Severini, (2005), Theorem 4.9). Consequently one needs to integrate
in the test statistic (9) over an interval containing the origin. Here for
convenience we integrate over $(0,\infty)$. A corresponding statistic of the
type $n\int_{t_{0}}^{\infty}D^{2}_{n}(t)e^{-\gamma t^{2}}{\rm{d}}t$, with
$t_{0}<0$, was also tried in simulations, but the results produced are very
similar and therefore not reported.

##  Appendix A.2 Computation of the MGF–based statistic

Starting from equation (7) we compute

$\displaystyle D^{2}_{n}(t)$ $\displaystyle=$ $\displaystyle(1+t)^{2}M^{\prime
2}_{n}(t)+[\widehat{p}-\widehat{\lambda}t(1+t)]^{2}M^{2}_{n}(t)+2(1+t)[\widehat{p}-\widehat{\lambda}t(1+t)]M_{n}(t)M^{\prime}_{n}(t)$
$\displaystyle=$
$\displaystyle(1+t)^{2}\frac{1}{n^{2}}\sum_{j,k=1}^{n}r_{j}r_{k}e^{t(r_{j}+r_{k})}+[\widehat{p}-\widehat{\lambda}t(1+t)]^{2}\frac{1}{n^{2}}\sum_{j,k=1}^{n}e^{t(r_{j}+r_{k})}+2(1+t)[\widehat{p}-\widehat{\lambda}t(1+t)]\frac{1}{n^{2}}\sum_{j,k=1}^{n}r_{j}e^{t(r_{j}+r_{k})}$

Now plugging the above expression for $D_{n}^{2}(t)$ into the test statistic
(9) and integrating term-by-term we arrive, after some grouping, at the
reported expression of $T_{n,\gamma}$, by making use of the integrals

$\int_{0}^{\infty}t^{k}e^{tx}e^{-\gamma t^{2}}dt,\ k=0,1,2,3,4.$

These integrals may be obtained in closed form by use of the standard normal
distribution function $\rm\Phi(\cdot)$. By way of example for $k=0,1,2$ we
have

$\int_{0}^{\infty}e^{tx}e^{-\gamma
t^{2}}dt=\sqrt{\frac{\pi}{\gamma}}{\rm{\Phi}}\left(\frac{x}{\sqrt{2\gamma}}\right)e^{\frac{x^{2}}{4\gamma}},$
$\int_{0}^{\infty}t\>e^{tx}e^{-\gamma
t^{2}}dt=\frac{1}{2\gamma}+\sqrt{\frac{\pi}{4\gamma^{3}}}{\rm{\Phi}}\left(\frac{x}{\sqrt{2\gamma}}\right)x\>e^{\frac{x^{2}}{4\gamma}},$

and

$\int_{0}^{\infty}t^{2}\>e^{tx}e^{-\gamma
t^{2}}dt=\frac{x}{4\gamma^{2}}+\sqrt{\frac{\pi}{16\gamma^{5}}}{\rm{\Phi}}\left(\frac{x}{\sqrt{2\gamma}}\right)(x^{2}+2\gamma)\>e^{\frac{x^{2}}{4\gamma}},$

respectively. For $k=3$ and $k=4$ the corresponding formulae are slightly more
complicated but still manageable.

##  Appendix A.3 CLS estimation equations

In order to obtain the CLS estimators described in Kopp and Mullahy (1993), we
consider the SFM with a corrected error term

$\displaystyle Y_{j}=X_{j}+v-(u-\mathbb{E}(u))=X_{j}+v-\delta\text{ , where
}\mathbb{E}(u)=cp$

and we solve the following system of equations:

$\displaystyle\frac{1}{n}\sum^{n}_{j=1}X^{\top}_{j}(Y_{j}-X_{j}\beta+cp)=0$
$\displaystyle\frac{1}{n}\sum^{n}_{j=1}(Y_{j}-X_{j}\beta+cp)^{2}-(\sigma_{v}^{2}+c^{2}p)=0$
$\displaystyle\frac{1}{n}\sum^{n}_{j=1}(Y_{j}-X_{j}\beta+cp)^{3}+2c^{3}p=0$
$\displaystyle\frac{1}{n}\sum^{n}_{j=1}(Y_{j}-X_{j}\beta+cp)^{4}-(3\sigma_{v}^{4}+6\sigma_{v}^{2}c^{2}p+3p(2+p)c^{4})=0$

##  Appendix A.4 Computation of the CF–based statistic

We follow analogous steps as those in Appendix A.2. Specifically we plug the
expression for $\Delta_{n}^{2}(t)$ given by (18) into the test statistic (17)
along with the formulae

$|\varphi_{n}(t)|^{2}=\frac{1}{n^{2}}\sum_{j,k=1}^{n}\cos t(r_{j}-r_{k})$
$|\varphi^{\prime}_{n}(t)|^{2}=\frac{1}{n^{2}}\sum_{j,k=1}^{n}r_{j}r_{k}\cos
t(r_{j}-r_{k})$

$S^{\prime}_{n}(t)C_{n}(t)-C^{\prime}_{n}(t)S_{n}(t)=\frac{1}{n^{2}}\sum_{j,k=1}^{n}r_{j}\cos
t(r_{j}-r_{k})$
$C^{\prime}_{n}(t)C_{n}(t)+S^{\prime}_{n}(t)S_{n}(t)=\frac{1}{n^{2}}\sum_{j,k=1}^{n}r_{j}\sin
t(r_{j}-r_{k})$

and then integrate term-by-term.

##  Appendix A.5 Integrals involved in the CF–test

The integral $I_{\nu,\gamma}$ may be computed as

$I_{\nu,\gamma}(z)=\frac{1}{2\gamma^{\frac{\nu+1}{2}}}{\rm{\Gamma}}\left(\frac{\nu+1}{2}\right)\
_{1}F_{1}\left(\frac{\nu+1}{2};\frac{1}{2};-\frac{z^{2}}{4\gamma}\right),$

while the integral $J_{\nu,\gamma}$ as

$J_{\nu,\gamma}(z)=\frac{z}{2\gamma^{1+\frac{\nu}{2}}}{\rm{\Gamma}}\left(1+\frac{\nu}{2}\right)\
_{1}F_{1}\left(1+\frac{\nu}{2},\frac{3}{2},-\frac{z^{2}}{4\gamma}\right),$

where 1$F_{1}$ stands for the Kummer confluent hypergeometric function; see
Gradshteyn and Ryzhik (1994), 3.952, eqn. 7 and 8 (p. 529), and §9.2 (p.
1084).


# Modeling Evacuee Behavior for Robot-Guided Emergency Evacuation

Mollik Nayyar1 and Alan R. Wagner2 This material is based upon work supported
by the National Science Foundation under Grant Number CNS-1830390 and
IIS-2045146. Any opinions, findings, and conclusions or recommendations
expressed in this material are those of the authors and do not necessarily
reflect the views of the National Science Foundation.1Mollik Nayyar is a
graduate student in aerospace engineering at the Robot Ethics and Aerial
Vehicles lab, The Pennsylvania State University<EMAIL_ADDRESS>R. Wagner
is an assistant professor in the Department of Aerospace Engineering at The
Pennsylvania State University<EMAIL_ADDRESS>
###### Abstract

This paper considers the problem of developing suitable behavior models of
human evacuees during a robot-guided emergency evacuation. We describe our
recent research developing behavior models of evacuees and potential future
uses of these models. This paper considers how behavior models can contribute
to the development and design of emergency evacuation simulations in order to
improve social navigation during an evacuation.

##  I INTRODUCTION

Our research is attempting to develop robots that are capable of helping
people evacuate during an emergency [1]. In order to accomplish this task the
robot must not only navigate social environments such as hallways, entry ways
and exits, it may also need to avoid collisions or, perhaps, even feint
collisions in order get people to move out of its way. Managing the collision
avoidance and social navigation problems in such environments is extremely
challenging for a variety of reasons [2]. First, emergencies are rare events
and difficult to predict. Hence, running real world experiments in which
robots attempt to guide evacuees during actual emergencies is not practical.
Moreover, emergencies can be unique. Earthquakes can damage buildings making
navigation maps unusable. Fires may require evacuation to distant or
improvised exits and avoidance of certain areas of a building. Finally, people
do not always react rationally during an emergency [3, 4]. Depending on the
type of emergency, humans may not evacuate at all, freeze and remain
motionless, or blindly comply with obviously incorrect information. Given the
wide range of potential environments, emergencies, robot designs, and
evacuees, evaluating the effectiveness of a robot’s guidance requires the use
of simulation experiments. Although simulation modeling of evacuations [5, 6,
7, 8] and human behavior modeling during emergencies is an active area of
research [3, 9, 10, 11], relatively few groups are exploring the development
of robot emergency guides [12, 13, 14].

Yet using simulation to evaluate robot designs, behaviors, or evacuation
strategies is also problematic for a number of reasons. For instance, the
behavior of other evacuees has a strong influence on the experiment’s human
subject (Fig. 1)[15]. In fact, during an evacuation, the behavior of other
evacuees is often a determining factor impacting when and how quickly a person
evacuates [4]. Our previous simulation experiments demonstrate that people
often follow the crowd, in spite of the robot’s suggestions [15]. Moreover, an
important aspect of designing a robot that guides people during an emergency
is understanding how people will react and respond to the robot [16]. Yet,
simulations may not generate the same visceral response as actual emergencies
[17]. Specifically, will evacuees follow an emergency guide robot? If so, how
closely will they follow the robot? Will they follow it through closed doors?
Will they hold the door open for the robot? Will they follow the robot if it
signals a change in direction? Clearly we have an ethical obligation to
thoroughly evaluate emergency evacuation robots prior to their deployment
[18].

Our research is currently working to address these problems. We are currently
in the process of developing behavior models of evacuees that we will then use
to create more realistic simulation environments (discussed more below). We
are also developing realistic simulations in virtual reality to improve the
realism of the experience for the purpose of better understanding social
navigation during an emergency. The remainder of this paper presents our
preliminary data and ongoing research on this topic.

##  II MODELING EVACUEE BEHAVIOR

This section details our method for creating and using behavior models of
evacuees.

###  II-A Creating Behavior Model

Given that evacuee behavior is strongly influence by the other people
experiencing the emergency, in order to begin to improve the accuracy of an
emergency simulation one must also simulate the behavior of the other evacuees
[8]. To do this we are currently running human subject experiments that
capture and record an evacuee’s behavior to a robot’s guidance directions when
an alarm has sounded in a physical environment. Specifically, we have naive
human subjects arrive at the experimental site (Fig. 2). These subjects are
briefly introduced to a guidance robot which leads them to a cubicle in order
to read and comment on an article. While reading the article a smoke alarm is
sounded. The subject has not been informed that a smoke alarm will go off
during the experiment. The subject then chooses whether to follow the robot’s
guidance directions to an unfamiliar exit or may exit the way they came.

Figure 1: A simulation experiment is depicted in which the subject chooses to
follow the robot during an emergency or groups of people running towards a
different exit. Our work shows that most people follow the crowd.  Figure 2:
The subject enters the experiment from Exit 1. Upon entering they are
introduced to the robot near Exit 1. The robot leads the person to point P1.
The subject sits at a cubicle in the Task Room. If the robot system under test
is a shepherding system, when the alarm sounds the robot leads the person from
point P1 through point P2 to Exit 2 stopping at point P3. If the robot system
under test is the handoff system, three individual robots are placed at points
P1, P2, and P3 again directing the subject to Exit 2.  Figure 3: A subject
following the robot to an exit.

We are using this general experimental paradigm to test two different types of
emergency evacuation robot systems. The first is a shepherding robot which is
a single robot that guides individuals or groups to exits. The second system
is a handoff system in which robots are deployed to the important evacuation
decision points (e.g. hall intersections). Each robot remains at that position
and points to either the next robot along an evacuation path or to the desired
exit. These two types of robot systems are being tested with either single
evacuees or groups of evacuees. We are recording video data capturing the
person’s movement through the environment. Soon we will explore some
variations of a simple evacuation. Specifically, we intend to capture the
subject’s behavior when the robot’s guidance directions change mid-evacuation,
when the robot must guide through a closed door, when confederates ignore the
robot’s directions, and when the robot instructs the subject to shelter in
place.

Nearly seventy human subjects have participated in our experiment to date
(Fig. 3). Importantly, our intent is to capture the descriptive statistics
representing the human subject’s behavior. This is in contrast to traditional
human-robot interaction (HRI) studies which tend to focus on hypothesis
testing of a limited number of hypotheses. These descriptive statistics
include measures of the person’s speed, local distance from the robot,
evacuation time, and the response time of the participant between the sounding
of smoke alarm and evacuation action. These measures will allow us to create
parameterized behavior models representing each human subject’s behavior
during the evacuation [19]. Although the form of these models has yet to be
determined, we envision some form of parameterized finite state machine that
can be used in conjunction with a specific set of experimental parameters to
generate motion commands representing prototypical evacuees behavior (average
motions), exemplar behavior (behavior tied to specific individual subjects),
or variations which randomize elements of the subject’s behavior.

###  II-B Using Behavior Models

As mentioned, we intend to use these behavior models as a tool to facilitate
the evaluation of emergency evacuation robots. In the simulation environment
pictured in Figure 1, the evacuating crowd runs at a fixed speed to an out of
sight exit. While this type of crowd behavior is possible, it is unlikely
given that most emergencies generate an initial sense of confusion accompanied
with information seeking behaviors [3].

The most straightforward use of these models is to vary the type of
environment (number of exits, number of corridors, etc.) and to use the
behavior models to estimate the evacuation time or other evacuation
statistics. For example, given a floor plan for school, we might use the
behavior models to estimate how long it will take students located in
particular classrooms to follow a robot to an exit. These estimates could then
be compared to a limited number of physical experiments.

Another potential way that these models could be used is to create different
demographic groups and estimate the amount of time a robot-guided evacuation
would require. For example, estimating how long it would take to evacuate able
bodied older adults from a nursing home or pregnant mothers could again inform
the design and use of the robots (e.g. speed of the robot, number of robots,
etc.). Of course, in order to evaluate these different demographic groups we
need to study participants of each type. Currently we have several older
adults that have participated, but we are not aware of any pregnant women that
have joined our study.

One final method for using these behavior models is to combine different
models to generate estimates of evacuation time for different scenarios. For
example, using a behavioral model for a shelter in place scenario, followed by
a behavior model for an evacuation, followed by a behavior model captured when
the robot redirects the person to a different exit. Although the individual
behavior models will have been generated by three different subjects, it may
be possible to produce a rough estimate of the evacuee’s behavior by simply
using different models for different events during the evacuation.

Experimentally, the resulting data could serve as a baseline estimate for in-
person experiments or it could govern the behavior of Non-Player Characters
(NPCs) in a simulation with a human subject. Hence, the human subject would be
surrounded by NPCs using these behavior models with the reactions of the human
subject being experimentally observed and recorded. Presented in an immersive
virtual reality simulation we believe that we will be able to produce
realistic environments with NPC behavior that reflects the behavior of true
human subjects.

##  III LIMITATIONS

We are making a number of important assumptions and there are certainly
limitations to our approach. One assumption is that the data generated by one
or many evacuees during a robot-guided evacuation will be correlated to and
predictive of other evacuees’ behavior during different emergencies. This
assumption only holds to a limited degree, however. Evacuee behavior after an
explosion and fire is unlikely to predict evacuee behavior during small, non-
threatening fire. Nevertheless, we believe that our research can provide a
type of foothold for representing how people evacuate in the presence of a
robot which can then be expanded to new situations and evacuee types.
Specifically, by capturing people’s behavior towards the robot in response to
an alarm, we record information about following speed, likelihood of
following, and response to robot guidance actions. This information is
necessary for designing evacuation robots and could be useful for related
applications, such as path clearing for evacuees [20].

Furthermore, we recognize that changes to the robot or the environment could
influence the resulting models and parameters. It will be necessary to
occasionally perform in-person experiments to evaluate if and how much the
behavior model has drifted from real world results. Simulation experiments may
inform robot and experimental design and focus the research on the most
promising hypotheses, yet cannot serve as a complete substitute for in-person
experiments.

Our long-term hope for this work is that a catalog of pedestrian behavior
models could be compiled allowing researchers to more accurately simulate the
actual social navigation behavior of individuals and groups during
emergencies. This catalog could serve as an important tool, perhaps allowing
HRI researchers to compare algorithms, systems, and behaviors for a fixed
context and reasonably accurate set of behavior models.

##  IV CONCLUSIONS

This paper has briefly introduced our ongoing research to create behavior
models of human evacuees while being guided to an exit by a robot during a
simulated emergency. The resulting behavior models will contribute to the
realism of our emergency evacuation simulations. This added realism will
assist in the design and use of emergency evacuation robots. Over time, we
hope to create a catalog of behavior models representing different evacuation
scenarios. This catalog could serve as a method to evaluate potential robot
social navigation algorithms applied to the emergency evacuation domain.

## References

  * [1] A. R. Wagner, “Robot-guided evacuation as a paradigm for human-robot interaction research,” _Frontiers in Robotics and AI_ , vol. 8, 2021.
  * [2] P. Trautman and A. Krause, “Unfreezing the robot: Navigation in dense, interacting crowds,” in _2010 IEEE/RSJ International Conference on Intelligent Robots and Systems_. IEEE, 2010, pp. 797–803.
  * [3] J. Leach, _Survival psychology_. Springer, 1994.
  * [4] E. Kuligowski, “Predicting human behavior during fires,” _Fire technology_ , vol. 49, no. 1, pp. 101–120, 2013.
  * [5] M. Pidd, F. De Silva, and R. Eglese, “A simulation model for emergency evacuation,” _European Journal of operational research_ , vol. 90, no. 3, pp. 413–419, 1996.
  * [6] E. D. Kuligowski, R. D. Peacock, B. L. Hoskins, _et al._ , _A review of building evacuation models_. US Department of Commerce, National Institute of Standards and Technology …, 2005.
  * [7] S. Sharma, H. Singh, and A. Prakash, “Multi-agent modeling and simulation of human behavior in aircraft evacuations,” _IEEE Transactions on aerospace and electronic systems_ , vol. 44, no. 4, pp. 1477–1488, 2008.
  * [8] G. Santos and B. E. Aguirre, “A critical review of emergency evacuation simulation models,” 2004.
  * [9] B. E. Aguirre, “Emergency evacuations, panic, and social psychology,” _Psychiatry: Interpersonal and Biological Processes_ , vol. 68, no. 2, pp. 121–129, 2005.
  * [10] X. Pan, _Computational modeling of human and social behaviors for emergency egress analysis_. Stanford University, 2006.
  * [11] X. Pan, C. S. Han, K. Dauber, and K. H. Law, “A multi-agent based framework for the simulation of human and social behaviors during emergency evacuations,” _Ai & Society_, vol. 22, no. 2, pp. 113–132, 2007.
  * [12] P. Robinette, W. Li, R. Allen, A. M. Howard, and A. R. Wagner, “Overtrust of robots in emergency evacuation scenarios,” in _2016 11th ACM/IEEE international conference on human-robot interaction (HRI)_. IEEE, 2016, pp. 101–108.
  * [13] E. Boukas, I. Kostavelis, A. Gasteratos, and G. C. Sirakoulis, “Robot guided crowd evacuation,” _IEEE Transactions on Automation Science and Engineering_ , vol. 12, no. 2, pp. 739–751, 2014.
  * [14] S. Zhang and Y. Guo, “Distributed multi-robot evacuation incorporating human behavior,” _Asian Journal of Control_ , vol. 17, no. 1, pp. 34–44, 2015.
  * [15] M. Nayyar and A. R. Wagner, “Effective robot evacuation strategies in emergencies,” in _2019 28th IEEE International Conference on Robot and Human Interactive Communication (RO-MAN)_. IEEE, 2019, pp. 1–6.
  * [16] P. Robinette, A. R. Wagner, and A. M. Howard, “Assessment of robot to human instruction conveyance modalities across virtual, remote and physical robot presence,” in _2016 25th IEEE International Symposium on Robot and Human Interactive Communication (RO-MAN)_. IEEE, 2016, pp. 1044–1050.
  * [17] ——, “Investigating human-robot trust in emergency scenarios: methodological lessons learned,” in _Robust Intelligence and Trust in Autonomous Systems_. Springer, 2016, pp. 143–166.
  * [18] A. R. Wagner, “Principles of evacuation robots,” in _Living with Robots_. Elsevier, 2020, pp. 153–164.
  * [19] H. C. Vorst, “Evacuation models and disaster psychology,” _Procedia Engineering_ , vol. 3, pp. 15–21, 2010.
  * [20] M. Nayyar and A. R. Wagner, “Aiding emergency evacuations using obstacle-aware path clearing,” in _2021 IEEE International Conference on Advanced Robotics and Its Social Impacts (ARSO)_. IEEE, 2021, pp. 7–14.