**6. Conclusions**

Energy is one of the most important inputs in manufacturing industries. It is a scarce input that is expensive in both monetary and environmental terms. Hence, both policymakers and businesses should consider the efficient use of this input in the long-term.

This study uses the stochastic frontier approach to measure energy use efficiency in the US manufacturing during the time period 1958–2011 using the NBER-CES Manufacturing Industry Database. When panel data are available as in our case, we advocate using the latest or the 4-component SF model. We concentrate on the most and least energy-intensive manufacturing industries. More specifically, we first define energy intensity as the costs of energy in total economic activity. Then for each of five decades, we identify the top 10% and bottom 10% energy-intensive industries. We apply the 4-component stochastic frontier model that decomposes overall efficiency into the long-term or persistent and short-term efficiencies. Our main findings suggest that energy use efficiency in US manufacturing hit hard by the oil shock in the 1970s and it did not rebound until the 2000s. The major culprit of the low overall energy use efficiency was structural inefficiency, a finding that goes hand in hand with the "energy paradox" (see, e.g., [4]). It seems that one of the ways to mitigate low levels of energy use efficiency should be to do more research along the lines of [5,47,48] to promote, adopt, and establish energy-efficient technologies as the new benchmark.

**Author Contributions:** Conceptualization, O.B. and S.C.K.; Methodology, O.B. and S.C.K.; Data collection S.C.K, formal analysis, O.B.; Writing—original draft preparation, O.B., Writing—Review and Editing, O.B. and S.C.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** All authors declare that they have no conflict of interest.

#### **Appendix A**




**Table A1.** *Cont.*

#### **Appendix B**

Here we describe how to estimate the model in (11). To facilitate the discussion, rewrite

$$\log q\_{it} = r(\mathbf{X}\_{it}, trend; \boldsymbol{\omega}) + v\_{0i} - \boldsymbol{u}\_{0i} + \boldsymbol{v}\_{it} - \boldsymbol{u}\_{it} \tag{A1}$$

as

$$\log q\_{it} = r(\mathbf{X}\_{it'} \, trend; \boldsymbol{\omega}) + \epsilon\_{0i} + \epsilon\_{it'}$$

where *it* = *vit* − *uit* and 0*<sup>i</sup>* = *v*0*<sup>i</sup>* − *u*0*<sup>i</sup>* decompose the error term into two 'composed error' terms (both of which contain a two-sided and a one-sided error terms). Assume the most general case where all four components are heteroskedastic

$$
\sigma\_{u\_{it}}^2 = \exp\left(\mathbf{z}\_{u\_{it}}\boldsymbol{\Psi}\_u\right), \quad i = 1, \cdots, n, \quad t = 1, \cdots, \tau,\\
T\_{i\prime} \tag{A2}
$$

$$
\sigma\_{\mathbf{u}\mathbf{q}}^2 = \exp\left(\mathbf{z}\_{\mathbf{u}\mathbf{q}}\boldsymbol{\Psi}\_{\mathbf{u}\mathbf{0}}\right), \quad \mathbf{i} = \mathbf{1}, \dots, \mathbf{n}, \tag{A3}
$$

$$
\sigma\_{v\_{\rm vir}}^2 = \exp\left(z\_{v\_{\rm vir}}\Psi\_v\right), \quad i = 1, \cdots, n, \quad t = 1, \cdots, \tau,\\
T\_{\rm ir} \tag{A4}
$$

$$
\sigma\_{v0i}^2 = \exp\left(\mathbf{z}\_{v0i}\boldsymbol{\Psi}\_{v0}\right), \quad i = 1, \cdots, n,\tag{A5}
$$

where *<sup>z</sup>uit* are the determinants of transient inefficiency, *<sup>z</sup>u*0*<sup>i</sup>* are the determinants of persistent inefficiency, and *<sup>z</sup>vit* and *<sup>z</sup>v*0*<sup>i</sup>* define the heteroskedasticity functions of the noise and random effects. The homoskedastic error component is easily derived from (A2–A5) by setting the vector of determinants to a constant. For example if *vit* is homoskedastic, *<sup>z</sup>vit* is a vector of ones of length <sup>∑</sup>*<sup>n</sup> <sup>i</sup>*=<sup>1</sup> *Ti*. The conditional density of *<sup>i</sup>* = (*i*1,..., *iTi* ) is given by

$$f\left(\boldsymbol{\epsilon}\_{i}|\boldsymbol{\epsilon}\_{0i}\right) = \prod\_{t=1}^{T\_{i}} \frac{2}{\sigma\_{it}} \boldsymbol{\phi}\left(\frac{\boldsymbol{\epsilon}\_{it}}{\sigma\_{it}}\right) \boldsymbol{\Phi}\left(\frac{\boldsymbol{\epsilon}\_{it}\boldsymbol{\lambda}\_{it}}{\sigma\_{it}}\right),$$

where *<sup>σ</sup>it* <sup>=</sup> [exp (*zuitψu*) <sup>+</sup> exp (*zvitψv*)]1/2 and *<sup>λ</sup>it* <sup>=</sup> [exp (*zuitψu*)/ exp (*zvitψv*)]1/2.

Integrate 0*<sup>i</sup>* (the distribution of which we know) out to get the unconditional density of *<sup>i</sup>*

$$f\left(\boldsymbol{\varepsilon}\_{i}\right) = \int\_{-\infty}^{\infty} \left[\prod\_{t=1}^{T\_{i}} \frac{2}{\sigma\_{it}} \boldsymbol{\Phi}\left(\frac{\boldsymbol{\varepsilon}\_{it}}{\sigma\_{it}}\right) \boldsymbol{\Phi}\left(\frac{\boldsymbol{\varepsilon}\_{it}\boldsymbol{\lambda}\_{it}}{\sigma\_{it}}\right)\right] \times \frac{2}{\sigma\_{0i}} \boldsymbol{\Phi}\left(\frac{\boldsymbol{\varepsilon}\_{0i}}{\sigma\_{0i}}\right) \boldsymbol{\Phi}\left(\frac{\boldsymbol{\varepsilon}\_{0i}\boldsymbol{\lambda}\_{0i}}{\sigma\_{0i}}\right) d\boldsymbol{\varepsilon}\_{0i} \boldsymbol{\lambda}\_{i}$$

where *<sup>σ</sup>*0*<sup>i</sup>* = [exp (*zu*0*<sup>i</sup> <sup>ψ</sup>u*0) <sup>+</sup> exp (*zv*0*<sup>i</sup> ψv*0)] 1/2 and *<sup>λ</sup>*0*<sup>i</sup>* = [exp (*zu*0*<sup>i</sup> <sup>ψ</sup>u*0)/ exp (*zv*0*<sup>i</sup> ψv*0)] 1/2. The log-likelihood function for the *i*-th observation of model (A1) is therefore given by

log *Li* (*β*, *ψu*0, *ψv*0, *ψu*, *ψv*) = log ⎡ ⎢ ⎢ ⎣ <sup>+</sup><sup>∞</sup> −∞ ⎛ ⎜⎜⎝ *Ti* ∏*t*=1 ⎧ ⎪⎪⎨ ⎪⎪⎩ 2 *σit φ rit* − 0*<sup>i</sup> <sup>σ</sup>it* ×Φ (*rit* − 0*i*)*λit <sup>σ</sup>it* ⎫ ⎪⎪⎬ ⎪⎪⎭ ⎞ ⎟⎟⎠ 2 *σ*0*i φ* 0*<sup>i</sup> σ*0*i* Φ 0*iλ*0*<sup>i</sup> σ*0*i* d0*<sup>i</sup>* ⎤ ⎥ ⎥ ⎦ <sup>=</sup> log <sup>+</sup><sup>∞</sup> −∞ *Ti* ∏*t*=1 2 *σit φ it σit* Φ *itλit <sup>σ</sup>it* × 2 *σ*0*i φ* 0*<sup>i</sup> σ*0*i* Φ 0*iλ*0*<sup>i</sup> σ*0*i* d0*<sup>i</sup>* , (A6)

where *it* <sup>=</sup> *rit* <sup>−</sup> (*v*0*<sup>i</sup>* <sup>+</sup> *<sup>u</sup>*0*i*) and *rit* <sup>=</sup> log *qit* <sup>−</sup> *<sup>r</sup>*(*Xit*, *trend*; *<sup>ω</sup>*). We rely on the Monte-Carlo integration as a method to approximate the integral in (A6). For estimation purposes, we write 0*<sup>i</sup>* = [exp (*zu*0*<sup>i</sup> ψu*0)] 1/2*Vi* + [exp (*zv*0*<sup>i</sup> ψv*0)] 1/2|*Ui*|, where both *Vi* and *Ui* are standard normal random variables. The resulting simulated log-likelihood function for the *i*-th observation is

log *L<sup>S</sup> <sup>i</sup>* (*β*, *ψu*0, *ψv*0, *ψu*, *ψv*) = log ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 1 *R R* ∑ *r*=1 ⎛ ⎜⎜⎜⎜⎝ *Ti* ∏*t*=1 ⎧ ⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎩ 2 *σit φ rit* <sup>−</sup> ([exp (*zu*0*iψu*0)] 1/2*Vir* + [exp (*zv*0*iψv*0)] 1/2|*Uir*|) *<sup>σ</sup>it* ×Φ [*rit* <sup>−</sup> ([exp (*zu*0*iψu*0)] 1/2*Vir* + [exp (*zv*0*iψv*0)] 1/2|*Uir*|)]*<sup>λ</sup> <sup>σ</sup>it* ⎫ ⎪⎪⎪⎪⎬ ⎪⎪⎪⎪⎭ ⎞ ⎟⎟⎟⎟⎠ ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ <sup>=</sup> log 1 *R R* ∑ *r*=1 *Ti* ∏*t*=1 2 *σ φ itr σ* Φ *itrλ σ* , (A7)

where *Vir* and *Uir* are *R* random deviates from the standard normal distribution, and *itr* = *rit* − ([exp (*zu*0*<sup>i</sup> ψu*0)] 1/2*Vir* + [exp (*zv*0*<sup>i</sup> ψv*0)] 1/2|*Uir*|). *<sup>R</sup>* is the number of draws for approximating the log-likelihood function. The full log-likelihood is the sum of panel-*i* specific log-likelihoods given in (A7).

We use the results of [39] to calculate persistent and time-varying cost efficiencies. Using the moment generating function of the closed skew normal distribution, the conditional means of *u*0*i*, *ui*1, ··· , *uiTi* are given by:

$$E\left(\exp\left\{t'u\_{i}\right\}|r\_{i}\right) = \frac{\overline{\Phi}\_{T\_{i}+1}\left(\mathcal{R}\_{i}r\_{i} + \mathcal{A}\_{i}t, \mathcal{A}\_{i}\right)}{\overline{\Phi}\_{T\_{i}+1}\left(\mathcal{R}\_{i}r\_{i}, \mathcal{A}\_{i}\right)} \times \exp\left(t'\mathcal{R}\_{i}r\_{i} + 0.5t'\mathcal{A}\_{i}t\right),\tag{A8}$$

where *r<sup>i</sup>* = ! *ri*1,...,*riTi* " , *<sup>A</sup>* <sup>=</sup> <sup>−</sup>[**1***Ti ITi* ], **<sup>1</sup>***Ti* is the column vector of length *Ti* and *<sup>I</sup>Ti* is the identity matrix of dimension *Ti*, the diagonal elements of *<sup>V</sup><sup>i</sup>* are [exp (*zu*0*<sup>i</sup> <sup>ψ</sup>u*0) exp (*zuitψu*)], **<sup>Σ</sup>***<sup>i</sup>* <sup>=</sup> exp (*zvitψv*)*ITi* <sup>+</sup> exp (*zv*0*<sup>i</sup> ψv*0)**1***Ti* **1** *Ti* , **<sup>Λ</sup>***<sup>i</sup>* <sup>=</sup> *<sup>V</sup><sup>i</sup>* <sup>−</sup> *<sup>V</sup>iA* (**Σ***<sup>i</sup>* <sup>+</sup> *AViA* ) <sup>−</sup><sup>1</sup> *AV<sup>i</sup>* <sup>=</sup> *V*−1 *<sup>i</sup>* <sup>+</sup> *<sup>A</sup>* **Σ**−<sup>1</sup> *<sup>i</sup> <sup>A</sup>* −<sup>1</sup> , *R<sup>i</sup>* = *ViA* (**Σ***<sup>i</sup>* + *AViA* ) <sup>−</sup><sup>1</sup> <sup>=</sup> **<sup>Λ</sup>***iA* **Σ**−<sup>1</sup> *<sup>i</sup>* , *φ<sup>q</sup>* (*x*, *μ*, **Ω**) is the density function of a *q*-dimensional normal variable with expected value *μ* and variance **Ω** and Φ*<sup>q</sup>* (*μ*, **Ω**) is the probability that a *q*-variate normal variable of expected value *<sup>μ</sup>* and variance **<sup>Ω</sup>** belongs to the positive orthant., *<sup>u</sup><sup>i</sup>* = (*u*0*i*, *ui*1, ... , *uiTi* ) , and <sup>−</sup>*<sup>t</sup>* is a row of the identity matrix of dimension (*Ti* <sup>+</sup> <sup>1</sup>). If <sup>−</sup>*<sup>t</sup>* is the *<sup>τ</sup>*-th row, Equation (A8) provides the conditional expected value of the *<sup>τ</sup>*-th component of the cost efficiency vector exp (−*ui*). In particular, for *τ* = 1, we get the conditional expected value of the persistent technical efficiency.
