Enhancing Efficiency: Halton Draws in the Generalized True Random Effects Model

Abstract

This paper measures the impact of the number of Halton draws in excess of $\lceil \sqrt{n}\rceil$ on technical efficiency in the generalized true random effects (four-component) stochastic frontier model estimated by simulated maximum likelihood. A substantial set of Monte Carlo simulations demonstrates that increasing the number of Halton draws to $\lceil n^{3/4}\rceil$ ($\lceil n^{2/3}\rceil$) decreases the mean squared error of the total technical efficiency estimates by $6.1$ ($4.9$) percent. Furthermore, increasing the number of Halton draws either improves or has no detrimental impact on correlation, mean squared error, relative bias, and upward bias for persistent, transient, and total technical efficiency. An energy sector application is included, to demonstrate how these issues can arise in practice, and how increasing Halton draws can improve parameter and efficiency estimates in empirical work.

1. Introduction

The generalized true random effects (GTRE) model log likelihood contains an integral with no closed-form solution, which can be difficult to estimate directly.1 This integral can be estimated by simulation via the Butler and Moffitt (1982) formulation, with the ultimate goal of obtaining precise efficiency estimates. Technical efficiency estimates have been utilized by a wide swath of applied research since the seminal cross-sectional model by Aigner et al. (1977) and Meeusen and van Den Broeck (1977).2 The four-component panel data model (Colombi 2010; Colombi et al. 2011; Colombi et al. 2014; Kumbhakar et al. 2014; Filippini and Greene 2016), which builds upon the earlier cross-sectional form, allows for the parsing of persistent and transient inefficiency, adding depth to efficiency studies. This paper identifies an under-appreciated issue of integral estimation in the efficiency literature, and it presents practitioners with a practical guide on both the integral estimation itself and the selection of the number of Halton draws above $\lceil \sqrt{n}\rceil$ therein.3

By increasing the number of Halton draws in the GTRE model estimated by simulated maximum likelihood, above $\lceil \sqrt{n}\rceil$, the impact on the time-varying, persistent, and total technical efficiency values can be assessed. Greene (2003) deployed 200 Halton draws for a sample size of 158 ($\lceil n^{1.047}\rceil$) for the normal-gamma model, while Bernstein (2020) used 537 Halton draws for a sample size of 3553 ($\lceil n^{0.769}\rceil$) for the four-component model. Rather than relying on the ad hoc use of relatively large numbers ($>>\lceil \sqrt{n}\rceil$), this work gives the precise benefit to increasing the number of Halton draws, which is relevant, given the computational cost associated with increasing draws.

Our simulations indicate that increasing the number of Halton draws to $\lceil n^{3/4}\rceil$ ($\lceil n^{2/3}\rceil$) decreases the mean squared error of the total technical efficiency estimates by $6.1$ ($4.9$) percent. Similarly, increasing the number of Halton draws decreases the mean squared error of the transient technical efficiency estimates by $2.1$ ($1.7$) percent and decreases the mean squared error of the persistent technical efficiency estimates by $3.2$ ($1.5$) percent. These simulations also indicate either improvement or no negative (statistical) impact for correlation, mean squared error, relative bias, and upward bias4 for persistent, transient, and total technical efficiency (on average, other things equal).

The remainder of this paper proceeds as follows: Section 2 describes the foundational concepts for Halton sequences, the methods of Filippini and Greene (2016), and the results for the bias and noise of the simulated estimators. Section 3 builds on the framework of Badunenko and Kumbhakar (2016) and extends it by considering the number of Halton draws. Section 4 considers an empirical energy sector example based on Bernstein (2020) with additional techniques deployed in the estimation. Section 5 concludes and suggests future areas of exploration.

2. Methods

Numerical integral estimation is a pervasive problem in economics and elsewhere, encompassing numerous decision points and trade-offs. Although this work focuses on the GTRE in the efficiency literature, this section may be of interest to econometric practitioners of the mixed logit and discrete response models more broadly, or, generally, to researchers estimating an integral with no closed-form solution to the Butler and Moffitt (1982) formulation. In this section, a formulation and discussion of Halton sequences is provided, followed by their specific application to the four-component stochastic frontier model estimated by simulated maximum likelihood and to general theoretical properties for simulated estimators.

2.1. Halton Sequences

The Halton sequence, developed in Halton (1960), is a quasi-random sequence of numbers that approximates draws from a uniform density on the closed interval from zero to one.5 From the universality of the uniform, Halton draws can be mapped into any continuous distribution for numerical methods, including integral approximation.

Halton sequences ($h_1,\dots ,h_n$) are constructed by beginning with a prime number greater than or equal to 2. Given the starting value of 2, one simply begins with the sequence $\{h_1,h_2\}=\{0,1/2\}$ then adds $1/2^2=1/4$ to that, in order to obtain $\{h_3,h_4\}=\{1/4,3/4\}$. Then, one simply adds $1/2^3=1/8$ to the entire sequence $\{h_1,\dots ,h_4\}$ to obtain $\{h_5,\dots ,h_8\}=\{1/8,5/8,3/8,7/8\}$, and so on. The general case for a Halton sequence for prime k is

$$h_{t+1}=\{h_t,h_t+{\displaystyle \frac{1}{k^t}},h_t+{\displaystyle \frac{2}{k^t}},\dots ,h_t+{\displaystyle \frac{k-1}{k^t}}\},$$

(1)

where each subsequence is a cycle of length k Train (2009).

The advantage of Halton sampling is illustrated in Figure 1, for sample sizes of $N\in \{150, 300, 600, 1000\}$. The Halton draws are mapped to the Gaussian distribution and are generated from the R package randtoolbox, using the halton call. The R built-in package stats creates ordinary normal draws using the rnorm call. Visual inspection of Figure 1 illustrates that Halton draws better approximate the normal distribution in contrast to ordinary draws.

In order to reduce correlation between/among Halton draws, the usual practice is to discard the first several draws in excess of the largest prime number used to generate the sequence(s) and to scramble the sequence(s) order, especially when moving beyond the one-dimensional (single-sequence) case (Train 2009). Additionally, for problems requiring more than one sequence, each sequence is generated with a different prime to reduce correlation. For the two-dimensional Halton draws in the simulation in Section 3, antithetic variates (see Section 4 for discussion) were not utilized, no draws were discarded or scrambled, and the 7th and 8th primes (17,19) were selected as the seeds for $r_i$ and $h_i$, respectively. However, for the empirical exercise in Section 4.1, the first and second primes were selected (2,3) and antithetic variates, discarding, and scrambling were all explored.

In a maximum simulated likelihood (MSL) context, there are different views in the literature as to whether draws should differ from one observation to the next (McFadden 1989; Train 2009) or if draws should be the same from observation to observation. McFadden (1989) suggests using several independent draws for each observation, to reduce “simulation chatter” and increase efficiency.6 Hence, the MSL in the McFadden (1989) approach is evaluated at more points and, therefore, may converge differently than MSL using the same draws, making the two approaches difficult to compare directly while holding draws fixed. A benefit to using the same draws for each observation is that the results are invariant to the order selection of the draws, which is not the case when the draws differ for some or all observations. In this context (and generally), more MSL evaluations improve the simulated log likelihood but increase computational intensity and model run times. Therefore, to economize on time and computational resources, this work utilized the same Halton draws for each observation.

2.2. The Generalized True Random Effects Model

The generalized true random effects (GTRE) model is a four-component model of the form

$$Y_{it}=e^{{\alpha }_0}\left(\prod _{j=1}^K{X}_{j,it}^{{\alpha }_j}\right)e^{r_i-h_i+v_{it}-u_{it}}.$$

(2)

The usual panel residuals are distributed normally with variance ${\sigma }_v^2$, where $v_{it}\sim N(0,{\sigma }_v^2)$ and firm heterogeneity $r_i\sim N(0,{\sigma }_r^2)$. The persistent component of inefficiency, $h_i$, is distributed half-normal with pre-truncation variance ${\sigma }_h^2$, and the time-varying component, $u_{it}$, is distributed half-normal, with pre-truncation variance ${\sigma }_u^2$. After taking the natural logarithm of both sides, where $log\left(X\right)=x$, Equation (2) becomes

$$y_{it}={\alpha }_0+\sum _{j=1}^K{\alpha }_j{x}_{j,it}+r_i-h_i+v_{it}-u_{it}={\alpha }_0+\sum _{j=1}^K{\alpha }_j{x}_{j,it}+{ϵ}_{it}.$$

(3)

The vectorized version of Equation (3) is estimated by simulated log likelihood (see Filippini and Greene (2016)), as follows:

$$\begin{array}{ccc}\hfill ln{L}_s(\alpha ,\lambda ,\sigma ,{\sigma }_r,{\sigma }_h)& & =\hfill \\ \hfill \sum _{i=1}^{N}ln{\displaystyle \frac{1}{R}}\sum _{r=1}^R\left\{\prod _{t=1}^T\right[{\displaystyle \frac{2}{\sigma }}& & \varphi \left((y_{it}-{\alpha }^{\prime }{x}_{it}-{\sigma }_r{r}_{ir}-{\sigma }_h|h_{ir}\left|\right)/\sigma \right)·\hfill \\ & & \Phi \left(-(y_{it}-{\alpha }^{\prime }{x}_{it}-{\sigma }_r{r}_{ir}-{\sigma }_h|h_{ir}\left|\right)\lambda /\sigma \right)\left]\right\}.\hfill \end{array}$$

(4)

In Equation (4), $\varphi \left(z\right)$ and $\Phi \left(z\right)$ describe the normal PDF and CDF, respectively; $\lambda ={\sigma }_u/{\sigma }_v$ and $\sigma ={({\sigma }_v^2+{\sigma }_u^2)}^{1/2}$. The normal half-normal density becomes a conditional density with the inclusion of $r_i$ and $h_i$; thus, to obtain the unconditional density, integration is required. The summation over R replaces an integral that eliminates $r_i$ and $h_i$ via MSL estimation and the law of large numbers. Hence, these plant-specific effects are simulated by Halton sequences, which greatly improves upon the ordinary normal draws, as illustrated in Figure 1.7 MSL estimation of the GTRE is a consistent, asymptotically normal, and efficient computational method for integral estimation, with a crucial decision point being the number of draws.

The technical efficiency estimation for the GTRE, as implemented in Filippini and Greene (2016), follows Colombi (2010); Kumbhakar et al. (2014), where

$$\mathbb{E}(exp(t·{(h_i,u_{i1},\dots ,u_{iT})}^{\prime })|{ϵ}_i)=\left[{\displaystyle \frac{{\overline{\Phi }}_{T+1}(R{\widehat{ϵ}}_i+\Lambda t^{\prime },\Lambda )}{{\overline{\Phi }}_{T+1}(R{\widehat{ϵ}}_i,\Lambda )}}\right]exp(tR{\widehat{ϵ}}_i+{\displaystyle \frac{1}{2}}t\Lambda t^{\prime }),$$

(5)

for $i=1,\dots ,N$, $t=(-1,0,\dots ,0),(0,-1,0,\dots ,0),\cdots ,(0,\dots ,0,-1)$, and $R=\Lambda A^{\prime }{\Sigma }^{-1}$. $\Lambda ={[V^{-1}+A^{\prime }{\Sigma }^{-1}A]}^{-1}$, $A=-\left[1_t\right|I_T]$ and $\Sigma ={\widehat{\sigma }}_v^2{I}_T+{\widehat{\sigma }}_r^2{1}_T{1}_T^{\prime }$. Additionally,

$$V=\left(\begin{array}{cc}{\widehat{\sigma }}_h^2& 0_T^{\prime }\\ 0_T& {\widehat{\sigma }}_u^2{I}_T\end{array}\right)$$

and ${\overline{\Phi }}_{T+1}(X,{\Sigma }_X)$ is the joint probability that a $T+1$-order normal vector is in the non-negative orthant, given mean X and variance ${\Sigma }_X$. The computation of Equation (5) uses multivariate normal integration.8

2.3. Bias and Noise of Simulated Estimators

To better understand the theoretical properties of the MSL, begin by setting $\theta$ to the parameter vector so that $l_s\left(\theta \right)=\partial ln{L}_s\left(\theta \right)/\partial \theta$. Next, let ${\theta }^{*}$ be the parameter values that maximize $ln{L}_s$, as $l_s\left({\theta }^{*}\right)=0$. The simulated gradient can be decomposed by adding and subtracting the true gradient at ${\theta }^{*}$ and the expectation of the simulated gradient at ${\theta }^{*}$ to arrive at the following:

$$l_s\left({\theta }^{*}\right)=l\left({\theta }^{*}\right)+\underset{=\mathrm{simulation}\mathrm{bias}}{\underbrace{[{\mathbb{E}}_r{l}_s\left({\theta }^{*}\right)-l\left({\theta }^{*}\right)]}}+\underset{=\mathrm{simulation}\mathrm{noise}}{\underbrace{[l_s\left({\theta }^{*}\right)-{\mathbb{E}}_r{l}_s\left({\theta }^{*}\right)]}}$$

(6)

Equation (6) (Train 2009) shows that the simulated gradient can be decomposed into simulation bias and simulation noise. The simulation bias normalized by the sample size can be shown to equal $Z\sqrt{n}{R}^{-1}$, where Z is constant. Thus, the simulation bias asymptotically requires that R increases faster than $\lceil \sqrt{n}\rceil$. For a finite sample, the simulation bias decreases with R. The simulation noise can be shown to be asymptotically normal with mean zero and variance $S{\left(nR\right)}^{-1}$, where S is a constant (Train 2009). Hence, for a given finite sample, simulation noise continues to decrease as R rises. Thus, the choice of R impacts the parameter estimates in practice, as is illustrated by the Taylor expansion:

$$\widehat{\theta }={\theta }^{*}-({\displaystyle \frac{\partial }{\partial \theta }}{l}_s\left({\theta }^{*}\right){)}^{-1}·l_s\left({\theta }^{*}\right)$$

(7)

Equation (7) (Train 2009) can be shown to be normally distributed with the same distribution as the maximum likelihood, yet the finite sample properties are clearly impacted by bias and noise in Equation (6). While the impact of this bias and noise is known for the parameters themselves, the impact on the efficiency estimates of the GTRE is the novel focus of this work.

3. Simulation Results

The simulation framework followed Badunenko and Kumbhakar (2016). A Cobb–Douglas production frontier with two inputs, $X_1$ and $X_2$, $Y=e^{{\alpha }_0}{X}_1^{{\alpha }_1}{X}_2^{1-{\alpha }_1}$, was deployed. The production frontier was constant returns to scale with ${\alpha }_0=0.3$ and ${\alpha }_1=0.4$. The covariates, $X_1$ and $X_2$, were truncated exponential distributions with truncation $log\left(2\right)$ and $log\left(10\right)$, respectively.

The $2^4=16$ different scenarios in

Table 1. Monte Carlo $\sigma$-combinations, as in Table 1 in Badunenko and Kumbhakar (2016).

Scenario	${\mathit{\sigma }}_{\mathit{h}}$	${\mathit{\sigma }}_{\mathit{u}}$	${\mathit{\sigma }}_{\mathit{r}}$	${\mathit{\sigma }}_{\mathit{v}}$	${\mathit{\lambda }}_0={\mathit{\sigma }}_{\mathit{h}}/{\mathit{\sigma }}_{\mathit{r}}$	$\mathit{\lambda }={\mathit{\sigma }}_{\mathit{u}}/{\mathit{\sigma }}_{\mathit{v}}$	$\mathit{\Lambda }={\mathit{\sigma }}_{\mathit{h}}/{\mathit{\sigma }}_{\mathit{u}}$
1	$0.04$	$0.04$	$0.04$	$0.04$	1	1	1
2	$0.04$	$0.2$	$0.04$	$0.04$	1	5	$0.2$
3	$0.2$	$0.04$	$0.04$	$0.04$	5	1	5
4	$0.2$	$0.2$	$0.04$	$0.04$	5	5	1
5	$0.04$	$0.04$	$0.04$	$0.2$	1	$0.2$	1
6	$0.04$	$0.2$	$0.04$	$0.2$	1	1	$0.2$
7	$0.2$	$0.04$	$0.04$	$0.2$	5	$0.2$	5
8	$0.2$	$0.2$	$0.04$	$0.2$	5	1	1
9	$0.04$	$0.04$	$0.2$	$0.04$	$0.2$	1	1
10	$0.04$	$0.2$	$0.2$	$0.04$	$0.2$	5	$0.2$
11	$0.2$	$0.04$	$0.2$	$0.04$	1	1	5
12	$0.2$	$0.2$	$0.2$	$0.04$	1	5	1
13	$0.04$	$0.04$	$0.2$	$0.2$	$0.2$	$0.2$	1
14	$0.04$	$0.2$	$0.2$	$0.2$	$0.2$	1	$0.2$
15	$0.2$	$0.04$	$0.2$	$0.2$	1	$0.2$	5
16	$0.2$	$0.2$	$0.2$	$0.2$	1	1	1

(see also Table 1 in Badunenko and Kumbhakar (2016)) were utilized for the variances and pre-truncation variances in $\{{\sigma }_r^2,{\sigma }_h^2,{\sigma }_u^2,{\sigma }_v^2\}$. The $\sigma$-values were either $0.04$ or $0.20$, which made the $\lambda ,{\lambda }_0$, and $\Lambda$ values either $0.2$, 1, or 5. Each scenario had 1000 Monte Carlo trials, with $N\in \{50,100,500\}$ and $t\in \{3,6,10\}$, and a maximum of 1000 function evaluations for each individual trial maximization.9 Letting $n=N·t$ and $\lceil \rceil$ be the ceiling operator, then the length of Halton values we considered were the following: ${\mathcal{H}}_{1/2}\left(n\right)=\lceil n^{1/2}\rceil$, ${\mathcal{H}}_{13/24}\left(n\right)=\lceil n^{13/24}\rceil$, ${\mathcal{H}}_{7/12}\left(n\right)=\lceil n^{7/12}\rceil$, ${\mathcal{H}}_{2/3}\left(n\right)=\lceil n^{2/3}\rceil$, and ${\mathcal{H}}_{3/4}\left(n\right)=\lceil n^{3/4}\rceil$.

The importance of this research with respect to run-times is illustrated in Figure 2. The log median second run-times were plotted against the nine combinations of $(t,N)$ for variance scenario 5, which represented ${\sigma }_v=0.2$, ${\sigma }_u={\sigma }_r={\sigma }_h=0.04$. For small sample sizes, the time difference for practical applications was very small, but for larger sample sizes the difference was enormous. In scenario 5, for example, when $N=50$ and $t=3$, median ${\mathcal{H}}_{1/2}$ took nearly 6 s and ${\mathcal{H}}_{3/4}$ took about 11 s, resulting in over 83 percent more time to run ${\mathcal{H}}_{3/4}$. However, for scenario 5, when $N=500$ and $t=10$, median ${\mathcal{H}}_{1/2}$ took 131 s and ${\mathcal{H}}_{3/4}$ took 1231, implying that ${\mathcal{H}}_{3/4}$ resulted in over 841 percent more run-time.10 In practical applications when there is not a known DGP and the panel might be unbalanced or there are outliers in the data these run-times tend to grow substantially, as was the case in Section 4.

The small-sample nature of when more Halton draws were the most beneficial is highlighted in Figure 3. Figure 3 presents the median MSE for total technical efficiency by different Halton spaces. The MSE levels were most separated for small sample sizes but tended to decrease for the largest of sample sizes, as the theory suggests.

Table 2. Internal meta-regressions: Three sources of log(mse).

	Persistent	Transient	Total
N = 100, t = 3	0.265 ***	0.430 ***	0.259 ***
	(0.006)	(0.005)	(0.005)
N = 100, t = 6	0.106 ***	0.154 ***	0.104 ***
	(0.006)	(0.005)	(0.005)
N = 50, t = 10	0.148 ***	0.097 ***	0.142 ***
	(0.006)	(0.005)	(0.005)
N = 50, t = 3	0.353 ***	0.589 ***	0.401 ***
	(0.006)	(0.005)	(0.005)
N = 50, t = 6	0.244 ***	0.275 ***	0.245 ***
	(0.006)	(0.005)	(0.005)
N = 500, t = 10	−0.167 ***	−0.156 ***	−0.198 ***
	(0.005)	(0.004)	(0.005)
N = 500, t = 3	0.054 ***	0.171 ***	−0.013 **
	(0.006)	(0.005)	(0.005)
N = 500, t = 6	−0.131 ***	−0.048 ***	−0.161 ***
	(0.005)	(0.004)	(0.005)
${\mathcal{H}}_{13/24}$	−0.030 ***	−0.011 ***	−0.040 ***
	(0.004)	(0.004)	(0.004)
${\mathcal{H}}_{7/12}$	0.003	−0.011 ***	−0.020 ***
	(0.005)	(0.004)	(0.004)
${\mathcal{H}}_{2/3}$	−0.016 ***	−0.017 ***	−0.050 ***
	(0.004)	(0.004)	(0.004)
${\mathcal{H}}_{3/4}$	−0.033 ***	−0.021 ***	−0.063 ***
	(0.004)	(0.004)	(0.004)
$\Lambda =1$	1.167 ***	−1.257 ***	−0.070 ***
	(0.004)	(0.003)	(0.003)
$\Lambda =5$	2.208 ***	−2.563 ***	0.305 ***
	(0.005)	(0.003)	(0.004)
${\lambda }_0=1$	−0.316 ***	0.532 ***	−0.043 ***
	(0.004)	(0.003)	(0.003)
${\lambda }_0=5$	−1.025 ***	1.154 ***	−0.653 ***
	(0.005)	(0.003)	(0.004)
$\lambda =1$	0.358 ***	−0.694 ***	0.003
	(0.004)	(0.003)	(0.003)
$\lambda =5$	0.845 ***	−1.406 ***	−0.349 ***
	(0.004)	(0.003)	(0.004)
Constant	−6.782 ***	−4.741 ***	−4.905 ***
	(0.007)	(0.005)	(0.006)
Observations	720,000	720,000	720,000
R2	0.197	0.337	0.118
Adjusted R2	0.197	0.337	0.118
Residual Std. Error (df = 719981)	1.200	0.955	1.040
F Statistic (df = 18; 719981)	9827 ***	20,372 ***	5342 ***
Base case average ${\mathcal{H}}_1$	0.0095	0.0049	0.0117

Notes: Robust standard errors clustered at the replication level in parenthesis. Significance levels are the following: * $p<$ 0.1; ** $p<$ 0.05; *** $p<$ 0.01. The base case in these regressions is ${\mathcal{H}}_{1/2}$, $\lambda ={\lambda }_0=\Lambda =0.2$, and N = 100, t = 10.

demonstrates that increasing the number of Halton draws to ${\mathcal{H}}_{3/4}$ (${\mathcal{H}}_{2/3}$) decreased the mean squared error (MSE) of the total technical efficiency estimates by $6.1$ ($4.9$) percent.11 Similarly, increasing the number of Halton draws decreased the mean squared error of the transient technical efficiency estimates by $2.1$ ($1.7$) percent, and increasing the number of Halton draws to ${\mathcal{H}}_{3/4}$ decreased the mean squared error of the persistent technical efficiency estimates by $3.2$ percent. The results for relative bias, upward bias, and correlation all showed a very similar pattern to MSE and are, therefore, included in the appendix, for the interested reader.

Table 2 and

Table A1. Internal meta-regressions: Three sources of relative bias.

	Persistent	Transient	Total
N = 100, t = 3	0.014 ***	0.003 ***	0.017 ***
	(0.0004)	(0.0002)	(0.0005)
N = 100, t = 6	0.003 ***	0.001 ***	0.003 ***
	(0.0004)	(0.0002)	(0.0004)
N = 50, t = 10	−0.006 ***	0.003 ***	−0.003 ***
	(0.0004)	(0.0002)	(0.0004)
N = 50, t = 3	0.005 ***	0.006 ***	0.011 ***
	(0.0004)	(0.0003)	(0.001)
N = 50, t = 6	−0.003 ***	0.003 ***	0.0004
	(0.0004)	(0.0002)	(0.0005)
N = 500, t = 10	0.001 ***	−0.001 ***	0.001 *
	(0.0003)	(0.0002)	(0.0004)
N = 500, t = 3	0.015 ***	−0.0001	0.015 ***
	(0.0003)	(0.0002)	(0.0004)
N = 500, t = 6	−0.0003	−0.001 ***	−0.001 ***
	(0.0003)	(0.0002)	(0.0004)
${\mathcal{H}}_{13/24}$	−0.006 ***	−0.00004	−0.006 ***
	(0.0003)	(0.0002)	(0.0004)
${\mathcal{H}}_{7/12}$	−0.005 ***	−0.0002	−0.006 ***
	(0.0003)	(0.0002)	(0.0004)
${\mathcal{H}}_{2/3}$	−0.013 ***	−0.0003 *	−0.014 ***
	(0.0003)	(0.0002)	(0.0003)
${\mathcal{H}}_{3/4}$	−0.015 ***	−0.0002	−0.016 ***
	(0.0003)	(0.0002)	(0.0003)
$\Lambda =1$	0.022 ***	−0.008 ***	0.015 ***
	(0.0002)	(0.0002)	(0.0003)
$\Lambda =5$	0.033 ***	−0.019 ***	0.014 ***
	(0.0004)	(0.0002)	(0.0004)
${\lambda }_0=1$	0.064 ***	0.0004 ***	0.065 ***
	(0.0002)	(0.0001)	(0.0003)
${\lambda }_0=5$	0.042 ***	0.001 ***	0.042 ***
	(0.0003)	(0.0002)	(0.0004)
$\lambda =1$	0.004 ***	0.034 ***	0.039 ***
	(0.0002)	(0.0002)	(0.0003)
$\lambda =5$	0.010 ***	0.022 ***	0.032 ***
	(0.0003)	(0.0002)	(0.0004)
Constant	−0.054 ***	−0.009 ***	−0.064 ***
	(0.0004)	(0.0003)	(0.001)
Observations	720,000	720,000	720,000
R2	0.155	0.117	0.116
Adjusted R2	0.155	0.117	0.116
Residual Std. Error (df = 719981)	0.076	0.047	0.092
F Statistic (df = 18; 719981)	7347 ***	5287 ***	5229 ***
Base case average ${\mathcal{H}}_1$	0.016	0.0069	0.0226

Notes: Robust standard errors clustered at the replication level in parenthesis. Significance levels are the following: * $p<$ 0.1; ** $p<$ 0.05; *** $p<$ 0.01. The base case in these regressions was ${\mathcal{H}}_{1/2}$, $\lambda ={\lambda }_0=\Lambda =0.2$, and N = 100, t = 10.

,

Table A2. Internal meta-regressions: Three sources of upward bias.

	Persistent	Transient	Total
N = 100, t = 3	0.067 ***	0.032 ***	0.045 ***
	(0.002)	(0.001)	(0.001)
N = 100, t = 6	0.014 ***	0.009 ***	0.009 ***
	(0.002)	(0.001)	(0.001)
N = 50, t = 10	−0.017 ***	0.023 ***	−0.004 ***
	(0.002)	(0.001)	(0.001)
N = 50, t = 3	0.033 ***	0.064 ***	0.038 ***
	(0.002)	(0.001)	(0.001)
N = 50, t = 6	0.007 ***	0.034 ***	0.008 ***
	(0.002)	(0.001)	(0.001)
N = 500, t = 10	0.012 ***	−0.016 ***	−0.0004
	(0.002)	(0.001)	(0.001)
N = 500, t = 3	0.037 ***	−0.004 ***	0.030 ***
	(0.002)	(0.001)	(0.001)
N = 500, t = 6	−0.004 **	−0.014 ***	−0.008 ***
	(0.002)	(0.001)	(0.001)
${\mathcal{H}}_{13/24}$	−0.033 ***	0.001	−0.020 ***
	(0.001)	(0.001)	(0.001)
${\mathcal{H}}_{7/12}$	−0.025 ***	0.0003	−0.018 ***
	(0.001)	(0.001)	(0.001)
${\mathcal{H}}_{2/3}$	−0.064 ***	0.0002	−0.044 ***
	(0.001)	(0.001)	(0.001)
${\mathcal{H}}_{3/4}$	−0.072 ***	0.001	−0.048 ***
	(0.001)	(0.001)	(0.001)
$\Lambda =1$	0.012 ***	0.031 ***	0.052 ***
	(0.001)	(0.001)	(0.001)
$\Lambda =5$	−0.038 ***	0.042 ***	0.039 ***
	(0.002)	(0.001)	(0.001)
${\lambda }_0=1$	0.220 ***	−0.037 ***	0.161 ***
	(0.001)	(0.001)	(0.001)
${\lambda }_0=5$	0.103 ***	−0.076 ***	0.068 ***
	(0.001)	(0.001)	(0.001)
$\lambda =1$	−0.016 ***	0.031 ***	0.100 ***
	(0.001)	(0.001)	(0.001)
$\lambda =5$	−0.038 ***	0.022 ***	0.099 ***
	(0.001)	(0.001)	(0.001)
Constant	0.436 ***	0.497 ***	0.313 ***
	(0.002)	(0.002)	(0.002)
Observations	720,000	720,000	720,000
R2	0.080	0.016	0.086
Adjusted R2	0.080	0.016	0.086
Residual Std. Error (df = 719981)	0.334	0.276	0.282
F Statistic (df = 18; 719981)	3460 ***	636 ***	3742 ***
Base case average ${\mathcal{H}}_1$	0.5674	0.5214	0.5337

Notes: Robust standard errors clustered at the replication level in parenthesis. Significance levels are the following: * $p<$ 0.1; ** $p<$ 0.05; *** p < 0.01. The base case in these regressions was ${\mathcal{H}}_{1/2}$, $\lambda ={\lambda }_0=\Lambda =0.2$, and N = 100, t = 10.

and

Table A3. Internal meta-regressions: Three sources of correlation.

	Persistent	Transient	Total
N = 100, t = 3	−0.084 ***	−0.052 ***	−0.038 ***
	(0.001)	(0.0003)	(0.001)
N = 100, t = 6	−0.030 ***	−0.015 ***	−0.017 ***
	(0.001)	(0.0002)	(0.001)
N = 50, t = 10	0.004 ***	0.002 ***	−0.010 ***
	(0.001)	(0.0002)	(0.001)
N = 50, t = 3	−0.066 ***	−0.066 ***	−0.046 ***
	(0.001)	(0.0004)	(0.001)
N = 50, t = 6	−0.033 ***	−0.015 ***	−0.025 ***
	(0.001)	(0.0002)	(0.001)
N = 500, t = 10	−0.0003	0.0001	0.021 ***
	(0.001)	(0.0001)	(0.001)
N = 500, t = 3	−0.073 ***	−0.044 ***	−0.021 ***
	(0.001)	(0.0002)	(0.001)
N = 500, t = 6	−0.007 ***	−0.011 ***	0.014 ***
	(0.001)	(0.0001)	(0.001)
${\mathcal{H}}_{13/24}$	0.015 ***	0.002 ***	0.006 ***
	(0.001)	(0.0002)	(0.0004)
${\mathcal{H}}_{7/12}$	0.007 ***	0.001 ***	0.005 ***
	(0.001)	(0.0002)	(0.0004)
${\mathcal{H}}_{2/3}$	0.022 ***	0.001 ***	0.010 ***
	(0.001)	(0.0002)	(0.0004)
${\mathcal{H}}_{3/4}$	0.029 ***	0.002 ***	0.014 ***
	(0.001)	(0.0002)	(0.0004)
$\Lambda =1$	0.137 ***	−0.049 ***	−0.176 ***
	(0.001)	(0.0002)	(0.0003)
$\Lambda =5$	0.295 ***	−0.082 ***	−0.032 ***
	(0.001)	(0.0003)	(0.001)
${\lambda }_0=1$	0.172 ***	0.046 ***	0.130 ***
	(0.0005)	(0.0002)	(0.0003)
${\lambda }_0=5$	0.557 ***	0.079 ***	0.444 ***
	(0.001)	(0.0002)	(0.0005)
$\lambda =1$	0.097 ***	0.335 ***	0.167 ***
	(0.001)	(0.0002)	(0.0004)
$\lambda =5$	0.190 ***	0.728 ***	0.483 ***
	(0.001)	(0.0002)	(0.0004)
Constant	−0.065 ***	0.143 ***	0.256 ***
	(0.001)	(0.0003)	(0.001)
Observations	720,000	720,000	720,000
R2	0.696	0.965	0.828
Adjusted R2	0.696	0.965	0.828
Residual Std. Error (df = 719981)	0.177	0.052	0.116
F Statistic (df = 18; 719981)	91,686 ***	1,097,208 ***	192,541 ***
Base case average ${\mathcal{H}}_1$	0.3666	0.4681	0.5274

Notes: Robust standard errors clustered at the replication level in parenthesis. Significance levels are the following: * $p<$ 0.1; ** $p<$ 0.05; *** $p<$ 0.01. The base case in these regressions was ${\mathcal{H}}_{1/2}$, $\lambda ={\lambda }_0=\Lambda =0.2$, and N = 100, t = 10.

present meta-regressions of the $16\times 3\times 3\times 1000\times 5=\mathrm{720,000}$ simulations. This design is similar to that of Andor et al. (2023), which advocates for internal meta-regressions in large-scale Monte Carlo simulations.

4. Empirical Exercise

An empirical exercise was considered, to illustrate how issues surrounding Halton draws might manifest in practice. This section also introduces some best practices for Halton draws, including discarding, scrambling, the inverse error function, and antithetic draws (abbreviated as “original and enhancements” or simply “enhanced”) into the original simulation estimation procedure/code (abbreviated as “original”). Discarding was performed by removing the first 1000 draws for both $r_i$ and $h_i$. Scrambling was introduced by taking 9999 samples of the Halton draw for $r_i$ and then selecting the one with the smallest absolute value of correlation with the $h_i$.12 The square root of two times the inverse error function ($\sqrt{2}{\mathit{erf}}^{-1}(·)$) was utilized for the $h_i$ rather than taking the absolute value of the standard normal draws to enhance precision, as $h_i$ was distributed half-normal. Finally, antithetic draws of the form $-r_i$ were considered by evaluating the likelihood at $-r_i$ for each i and then taking a simple average with the likelihood evaluated at the usual $r_i$.13

4.1. Production Function Application

A Cobb–Douglas production function from Bernstein (2020) was utilized on a panel of United States natural gas plants. The data came from the FERC Form 1: Electric Utility Annual Report. In order to demonstrate the importance of having a sufficiently large number of draws, the following model was estimated:

$$y_{it}={\alpha }_0+\sum _{j=1}^4{\alpha }_j{x}_{j,it}+\sum _{m=1}^3{\beta }_m{c}_{m,it}+\underset{={ϵ}_{it}}{\underbrace{r_i-h_i+v_{it}-u_{it}}},$$

(8)

with Halton spaces utilized for powers from $\lceil n^{50/100}\rceil$ (baseline $\lceil \sqrt{n}\rceil$) to $\lceil n^{105/100}\rceil$ with all $5/100$ unit increments in between, as well as the setting of $\lceil n^{0.769}\rceil$ ($\lceil n^{0.77}\rceil$, henceforth) from Bernstein (2020).14 In this example, $y_{it},x_{j,it}$, and $c_{m,it}$ were all in natural logarithms, with the upper case variables in levels. Furthermore, $Y_{it}$ was the plant output in kWh, the production factors $X_{j,it}$ included plant capacity MW (h) ($X_1=\mathrm{K}$), number of employees ($X_2=\mathrm{L}$), natural gas in Mcf ($X_3=\mathrm{NG}$), and barrels of oil ($X_4=\mathrm{Oil}$). Three control variables $c_{m,it}$ were also included.15 The sample size was 3553, which was an unbalanced panel with 448 identified plants and time periods ranging from 2 to 23, with an average of $7.9$ reporting years per entity. Summary statistics for the data are provided in

Table 3. Summary statistics for natural gas electric generating plants from 1994–2016.

	Y	K	L	Oil	NG
units	kWh/103	MW (h)	employees	barrels/103	Mcf/103
mean	483,178	432	26	76	6072
sd	576,553	421	28	619	10,965
min	1018	6	1	0	2
max	2,145,656	2764	214	24,331	154,002
Observations $=3553$ (unbalanced panel)

Notes: Y, Oil, and NG are each divided by 10³. MW stands for megawatts, kWh for kilowatt-hours, h for hour, and Mcf for thousands of cubic feet.

:

For each model run, the bobyqa optimizer of the minqa package was utilized with 1,000,000 maximum function evaluations and, subsequently, the optim call in the stats package calculated the numerical Hessian for computing standard errors.16 For the simulations, the optim call was utilized with a maximum of 1000 function evaluations.

In

Table 4. Selected coefficients, TE averages, and run times.

x	$\lceil {\mathit{n}}^{\mathit{x}/100}\rceil$	$\mathit{\lambda }$	$\mathit{\sigma }$	${\mathit{\sigma }}_{\mathit{r}}$	${\mathit{\sigma }}_{\mathit{h}}$	TE-p	TE-tr	Hours	log ℓ
Original
50	60	1.38 ***	0.57 ***	0.81 ***	0.63 ***	0.63	0.73	0.57	−2707.63
55	90	1.39 ***	0.57 ***	0.82 ***	0.49 ***	0.69	0.73	0.63	−2717.54
60	136	1.40 ***	0.57 ***	0.72 ***	0.47 ***	0.70	0.73	0.77	−2690.55
65	204	1.39 ***	0.57 ***	0.74 ***	0.46 ***	0.71	0.73	1.02	−2688.93
70	306	1.39 ***	0.57 ***	0.81 ***	0.00	1.00	0.73	1.18	−2688.94
75	461	1.39 ***	0.57 ***	0.81 ***	0.00	1.00	0.73	2.11	−2690.57
77	537	1.41 ***	0.57 ***	0.74 ***	0.29 ***	0.80	0.73	2.51	−2673.41
80	693	1.41 ***	0.57 ***	0.73 ***	0.34 *	0.77	0.73	4.07	−2674.24
85	1043	1.40 ***	0.57 ***	0.73 ***	0.24 ***	0.83	0.73	5.81	−2667.48
90	1569	1.40 ***	0.57 ***	0.73 ***	0.24 ***	0.83	0.73	7.20	−2668.59
95	2361	1.40 ***	0.57 ***	0.67 ***	0.39 ***	0.75	0.73	10.79	−2655.32
100	3553	1.40 ***	0.57 ***	0.67 ***	0.41 ***	0.74	0.73	21.48	−2656.65
105	5348	1.40 ***	0.57 ***	0.67 ***	0.40 ***	0.74	0.73	20.91	−2657.81
Original and Enhancements
50	60	1.28 ***	0.56 ***	0.71 ***	0.32 ***	0.78	0.74	0.90	−2661.90
55	90	1.36 ***	0.57 ***	0.73 ***	0.12 **	0.91	0.73	0.79	−2656.36
60	136	1.39 ***	0.57 ***	0.74 ***	0.00	1.00	0.73	0.90	−2662.65
65	204	1.40 ***	0.57 ***	0.73 ***	0.08	0.94	0.73	1.58	−2660.55
70	306	1.40 ***	0.57 ***	0.73 ***	0.07	0.95	0.73	2.32	−2661.22
75	461	1.40 ***	0.57 ***	0.74 ***	0.00	1.00	0.73	3.72	−2663.09
77	537	1.39 ***	0.57 ***	0.72 ***	0.20 *	0.86	0.73	5.04	−2661.04
80	693	1.40 ***	0.57 ***	0.73 ***	0.11 *	0.92	0.73	6.18	−2663.44
85	1043	1.40 ***	0.57 ***	0.73 ***	0.14 **	0.90	0.73	6.73	−2664.02
90	1569	1.40 ***	0.57 ***	0.72 ***	0.14 *	0.90	0.73	12.43	−2657.93
95	2361	1.41 ***	0.57 ***	0.70 ***	0.21 **	0.85	0.73	16.31	−2657.61
100	3553	1.40 ***	0.57 ***	0.68 ***	0.30 ***	0.80	0.73	27.60	−2654.04
105	5348	1.40 ***	0.57 ***	0.70	0.00	1.00	0.73	38.64	−2655.33

Notes: Significance levels are given by 0.01 ***, 0.05 **, and 0.10 * from the Student’s t-Distribution on 3553 − 12 = 3541 degrees of freedom; x refers to the x in $\lceil n^{x/100}\rceil$. The largest log likelihood (log ℓ) values are highlighted in each set of simulations. “Enhancements” refers to discards, antithetic draws, $\sqrt{2}{\mathit{erf}}^{-1}(·)$, and shuffling. Appendix Table A4 provides further testing of modifications made to the original case. The largest log likelihood (log ℓ) values are highlighted in pink in each set of simulations.

, the efficiency parameters are presented for each set of Halton draws ($\lceil n^{x/100}\rceil$) along with the average persistent (TE-p) and transient (TE-tr) technical efficiencies and the number of hours each model took to run for both the original and enhanced settings. The parameters $\sigma$ and $\lambda$ were estimated extremely well from model to model (as were the $\alpha$ and $\beta$ parameters, excluding the intercept). The consistent estimation of $\sigma$ and $\lambda$ from run to run allowed for nearly identical estimates of the transient technical efficiency, with a minimum correlation of $0.964$ for the original and $0.991$ for the enhanced in the 13-by-13 Pearson correlation matrices; ${\sigma }_r$ varied up to 22 percent for the original and 9 percent for the enhanced from the smallest to the largest case, while ${\sigma }_h$ was zero twice in the original and zero thrice in the enhanced, making the estimation of the persistent TE more variable from run to run. Interestingly, ${\sigma }_h$ had a much tighter range of 0.00–0.32 in the enhanced model relative to 0.00–0.63 in the original, ${\sigma }_r$ had a tighter range of 0.68–0.74 in the enhanced model relative to 0.67–0.82 in the original, $\sigma$ was nearly the same for all runs, and $\lambda$ only differed somewhat in the two smallest Halton settings for the enhanced estimation.

In Figure 4, six parameter plots are provided, in order to illustrate issues of measuring certain parameters, with solid parameter lines and 95 percent confidence bands in dotted lines.17 The figures highlight the strong ability of all the models to tightly estimate the coefficient for the log of NG, and, generally, all the $\alpha$s and $\beta$s look similar.18 The models’ difficulty in estimating the intercept is very striking for lower values of Halton draws, and this illustrates the importance of precisely estimating ${\sigma }_h$ and ${\sigma }_r$ for the constant term, as the composed error term does not have zero mean. In every original plot, the slope of the parameter appears to level off towards zero, indicating that $\lceil n^{95/100}\rceil$ sufficed. For the enhanced estimation, ${\sigma }_r,{\sigma }_h$ and the intercept do not appear to level off, as a zero for ${\sigma }_h$ repeatedly appears, likely due to a local maximum of the log likelihood function that is very close to nonzero values for ${\sigma }_h$.19 Consequently, in this setup, it is recommended to utilize the estimates that maximize the log likelihood. The original case of $\lceil n^{95/100}\rceil$ and $\lceil n^{100/100}\rceil$ for the enhanced each maximized the log likelihood, as highlighted in Table 4 (this is denoted as $x\left[\max \right]$, so that for the original $x\left[\max \right]=95$, and $x\left[\max \right]=100$ for the enhanced).

Considering the case of $\lceil n^{x\left[\max \right]/100}\rceil$ Halton draws to the sample size as the true (oracle) efficiency values for the original setting, we were able to compute the MSE of the transient TE as of the $\lceil \sqrt{n}\rceil$ case of $0.00038$, an upward bias of $0.35$, a relative bias of $-0.002$, and a correlation of $0.98$. As for the persistent TE, there was a larger divergence generally as MSE was $0.014$, upward bias was $0.00$, relative bias was $-0.16$, and correlation was $0.97$. Similarly, for the enhanced setting where $x\left[\max \right]=100$, the MSE of the transient TE was $0.000098$, the upward bias was $0.97$, the relative bias was $0.013$, and the correlation was $1.00$. As for the persistent TE, the MSE was $0.00017$, the upward bias was $0.00$, the relative bias was $-0.016$, and the correlation was $1.00$. This result is strongly suggestive of the superior performance of the enhanced estimation methodologies in recovering similar root n results to the $x\left[\max \right]$ case.

Large discrepancies in TE estimates could be problematic for European energy efficiency regulators, as are general discrepancies writ large. Andor et al. (2019) notes that regulators prefer to overestimate firm efficiency measures because if “the efficiency of such a firm is underestimated, the consequence is a cost saving requirement, which is too strict. In the long run this could lead to insolvency through no fault of the firm”. As is shown herein for the enhanced case, the upward bias of the transient TE was $0.97$ for the $\lceil \sqrt{n}\rceil$ case, which implied that the majority of the $\lceil \sqrt{n}\rceil$ estimates were likely overestimated. As for the persistent TE for both the original and enhanced cases, all of the $\lceil \sqrt{n}\rceil$ estimates were underestimates. Hence, from the Andor et al. (2019) regulatory perspective, some combination approach may be desirable to increase the likelihood of overestimated TE measures.

Figure 5 highlights the differences between the three technical efficiency estimates from the case of $\lceil n^{x\left[\mathit{max}\right]/100}\rceil$ Halton draws on the y-axis and $\lceil \sqrt{n}\rceil$ on the x-axis. The transient TE estimates being very highly correlated, they also tended to stay on the $y=x$ line as they had very low relative bias amounting to just $0.20$ percent, on average, for the original and $1.27$ percent for the enhanced. For the original, a half dozen or so estimates appeared far away from the $y=x$ line, which is not surprising in a sample of 3553 observations, yet these were eliminated when the enhancements were utilized. For the persistent TE, it is clear why upward bias was zero, as full separation of the data mass from the 45 degree line is noticeable, albeit less pronounced for the enhanced case. Overall, the enhanced estimation shows much tighter efficiency estimates across all three types of TE from the root n to the “oracle” settings relative to the original estimation.

5. Conclusions

This paper measured the impact of the number of Halton draws in excess of $\lceil \sqrt{n}\rceil$ on three technical efficiency measures in the four-component stochastic frontier setting. Increasing the number of Halton draws to $\lceil n^{3/4}\rceil$ ($\lceil n^{2/3}\rceil$) decreased the mean squared error of the total technical efficiency estimates by $6.1$ ($4.9$) percent, and either improved ($\uparrow /\downarrow$) or had no damaging impact on correlation (↑), mean squared error (↓), relative bias (↓), and upward bias (↓) for the three technical efficiency measures. These findings suggest that efficiency metrics are improved by increasing the number of Halton draws well beyond what is required for consistency, especially for small samples. As a result, practitioners should implement best practices for simulating integrals such as discarding and scrambling, and they should balance the time it takes them to simulate with the precision they require for their analysis, as illustrated in the empirical example.

Future work might examine higher dimensions for Halton sequences such as the mixed logit model or the performance of the (deterministic) quadrature for integral approximations in the present setting. Another avenue of exploration is the presence and absence of determinants of inefficiency and the impact of scrambling, discarding, and antithetic draws on the simulation results herein. Furthermore, examining more draws per observation and utilizing other quasi-random sequences, such as the Sobol sequences or Hammersley sets, in lieu of Halton sequences would provide new lines of inquiry.