*Texas A&M University
**University of Sydney, CIREQ & CEBDA
June 2026
\[ Y_i = f(X_i) + V_i - U_i, \qquad U_i \ge 0 \]
parametric \(f(\cdot)\) interpretable but rigid
nonparametric \(f(\cdot)\) flexible but can violate production theory, aka shape constraints
Nonparametric SFM: Fan et al. (1996), Kumbhakar et al. (2007), Parmeter and Racine (2012), Ferrara and Vidoli (2017), Tran et al. (2023)
DNN for SFM: Wang (1996, 2003), Tsionas and Mamatzakis (2019), Tsionas (2022), Tsionas et al. (2023), Kutlu and Mao (2024)
Main idea: use deep neural networks, but build monotonicity and concavity into the architecture rather than hoping the fitted curve behaves itself.
\(\Rightarrow\) shape-awareshape-aware learninglearning
Flexible nonparametric approximation via hidden layers and activation functions.
Approximation theory (Hornik et al., 1989; Cybenko, 1989): single-layer perceptron (shaller NN) approximates any continuous function to a desired degree of accuracy.
Multi-layer networks (deep NN): greater flexibility and generalization especially for high-dimensions, interactions and nonlineary (Hornik, 1991; Poggio et al., 2017)
Layer 0: inputs \(\mathbf{h}^{(0)}=\mathbf{X}=(X_1,\ldots,X_K)^\top\) inputs of production
Layer \(L\): output \(Y = {\mathbf{h}^{(L)}}+ V-U\),
where \({\mathbf{h}^{(L)}}\) is defined recursively using an activation function \(\sigma^{(\ell)}\): \[ \mathbf{h}^{(\ell)} = \sigma^{(\ell)}\left(W^{(\ell)}\mathbf{h}^{(\ell-1)}+\mathbf{w}_0^{\ell}\right), \qquad \ell = 1,\ldots,L. \]
\[ \underbrace{ \begin{pmatrix} h_1^{(\ell)}\\ h_2^{(\ell)}\\ \vdots\\ h_{K^{(\ell)}}^{(\ell)}\\ \end{pmatrix}}_{\mathbf{h}^{(\ell)}} = \begin{pmatrix} \sigma^{(\ell)}(o_1^{(\ell)})\\ \sigma^{(\ell)}(o_2^{(\ell)})\\ \vdots\\ \sigma^{(\ell)}(o_{K^{(\ell)}}^{(\ell)})\\ \end{pmatrix}, %\equiv \sigma^{(\ell)}(\mathbf{o}^{(\ell)}), \text{ with }\textbf{activation function } \sigma^{(\ell)} \qquad\qquad %\text{Signal to layer $\ell$: } \underbrace{ \begin{pmatrix} o_1^{(\ell)}\\ o_2^{(\ell)}\\ \vdots\\ o_{K^{(\ell)}}^{(\ell)}\\ \end{pmatrix}}_{\textbf{o}^{(\ell)}} = \underbrace{ \begin{pmatrix} w^{(\ell)}_{1,1} & w^{(\ell)}_{1,2} & \ldots &w^{(\ell)}_{1,K^{(\ell-1)}}\\ w^{(\ell)}_{2,1} & w^{(\ell)}_{2,2} & \ldots &w^{(\ell)}_{2,K^{(\ell-1)}}\\ \vdots & \vdots & &\vdots\\ w^{(\ell)}_{K^{(\ell)},1} & w^{(\ell)}_{K^{(\ell)},2} & \ldots &w^{(\ell)}_{K^{(\ell)},K^{(\ell-1)}}\\ \end{pmatrix} }_{W^{(\ell)}} \underbrace{ \begin{pmatrix} h_1^{(\ell-1)}\\ h_2^{(\ell-1)}\\ \vdots\\ h_{K^{(\ell-1)}}^{(\ell-1)} \end{pmatrix}} _{\mathbf{h}^{(\ell-1)}} + \underbrace{ \begin{pmatrix} w_{01}^{(\ell)}\\ w_{02}^{(\ell)}\\ \vdots\\ w_{0K^{(\ell)}}^{(\ell)} \end{pmatrix}} _{\textbf{w}_0^{(\ell)}} \]
\[Y_i = \text{DNN}(\mathbf{X}_i;\mathbf{w}) + V_i-U_i, V_i \sim N(0,\sigma_v^2) \perp U_i \sim |N(0,\sigma_u^2)|, i=1, \ldots, n\]
and \(\mathbf{w}=\{W^{(\ell)},\mathbf{w}_0^{(\ell)}\}\) has all unknown weights \(W^{(\ell)}\) and biases \(\mathbf{w}_0^{(\ell)}\)
A production frontier should respect:
More input should not reduce feasible output:
\[ \frac{\partial f(X)}{\partial X_k} \ge 0 \]
for each input \(k\).
Marginal returns should not increase indefinitely:
\[ f(\lambda X + (1-\lambda)Z) \ge \lambda f(X)+(1-\lambda)f(Z) \]
for \(\lambda\in[0,1]\).
How to make DNN frontiers monotone and concave by construction?
Common activation functions:
\[ \sigma(t)=\max(0,t) \]
\[ \sigma(t)= \begin{cases} t, & t>0 \\ \alpha(e^t-1), & t\le 0 \end{cases} \]
Many other activation functions exist (see,e.g., LeCun et al., 1989; Nair and Hinton, 2010). All convex.
None produces DNNs that respect the shape constaints.
A flexible network can fit the data well and still produce economically nonsensical frontier shapes.
DGP (solid): Cobb-Douglas
Estimated (dashed): DNN with \(L=2\)
Point: DNN fit is neither monotone nor concave
A DNN is a valid production function if, for each hidden layer \(\ell\):
\[ W^{(\ell)} \ge 0, \qquad \mathbf{w}_0^{(\ell)} \ge 0. \]
Compositions of non-decreasing concave functions with non-negative affine maps remain monotone and concave.
\[ \text{DNN}(\mathbf{X};\mathbf{w}_0) = (\sigma^{(L)} \circ W^{(L)}\circ\cdots \sigma^{(1)} \circ W^{(1)})(\mathbf{X}) \]
Built-in architecture rather than after-estimation adjustements, rejection-sampling, etc
\[ \sigma(t)= \begin{cases} \alpha t, & t\ge 0,\quad 0<\alpha<1 \\ t, & t<0 \end{cases} \]
Piecewise linear and concave.
\[ \sigma(t)= \begin{cases} -\alpha(e^{-t}-1), & t\ge 0 \\ t, & t<0 \end{cases} \]
Smooth, concave, and non-decreasing.
These activations should be used when estimating DNN for production functions.
They also retain universal approximation properties of traditional DNNs – see paper.
They also have no vanishing gradient problem – see paper.
A longer list of activations in the paper, each having its own computational appeals
DGP (solid): Cobb-Douglas
Estimated (dashed): DNN with \(L=2\)
Point: DNN fit is both monotone and concave
\[\varepsilon_i= Y_i - DNN(\mathbf{X}_i;\mathbf{w}) \qquad \theta = (\mathbf{w}, \sigma_v^2, \sigma_u^2)\] \[ \text{Loss }L({\theta}) = \sum\limits_{i=1}^n \{-\ell_i({\theta}|\boldsymbol X_i,Y_i)\}, \text{ where }\ell_i({\theta}|\boldsymbol X_i,Y_i) = -\frac{1}{2}\log\left(\frac{\pi}{2}\right) - \frac{1}{2}\log\left(\sqrt{\sigma_v^2+\sigma_u^2}\right) + \log\Phi\left( -\frac{\sigma_u\varepsilon_i}{\sigma_v\sqrt{\sigma_v^2+\sigma_u^2}} \right) - \frac{\varepsilon_i^2}{{2(\sigma_v^2+\sigma_u^2})}. \]
Training by mini-batch stochastic gradient descent: \[ \theta_{t+1} \leftarrow \theta_t - a\frac{1}{|B_t|} \sum_{i\in B_t}\nabla_\theta L(\theta) \] with automatic differentiation and backpropagation, batch size \(|B_t|\), learning rate \(a\).
Inference using bootstrap.
Using Pytorch library in Python 3.
Note: unconstrained loss minimization problem with no extra computational tools beyond standard DNN.
Leave‐one‐out cross‐validation (LOOCV) to tune over a grid of candidate architectures varying in depth (\(L=1, 2, 3\)) and in width (\(2^2\), \(2^3\), \(2^4\) and \(2^5\) hidden nodes per layer) and over a grid of \(a\) between 0.1 and 0.01, with \(|B|=n\) when feasible.
DNNs are flexible, but not naturally interpretable. The paper uses SHAP values:
\[ DNN(x;\widehat w) = \phi_0 + \sum_{k=1}^K \phi_k(x). \]
For input \(j\) from a full set \([K]\):
\[ \phi_j(x)= \frac{1}{K}\sum_{S\subseteq[K]\setminus\{j\}} {K-1 \choose |S|}^{-1} \left[\nu_{S\cup\{j\}}(x)-\nu_S(x)\right]. \]
SHAP gives local input contributions, not marginal products or elasticities.
Shows how the prediction changes as inputs are included against a background distribution.
Note: new results on inference for SHAP:
“Shapley Values: Paired-Sampling Approximations” by M Mayer and M Wuethrich
“Profit and loss decomposition in continuous time and approximations” by G Junike, H Stier and M Christiansen
Three DGPs, two Cobb-Douglas, one Translog
Single relevant input:
\[ Y_i = \beta_0 X_{1i}^{\beta_1}\exp(V_i-U_i). \]
Includes irrelevant regressors:
\[ Y_i = \beta_0 X_{1i}^{\beta_1}X_{2i}^{0}X_{3i}^{0}X_{4i}^{0}\exp(V_i-U_i). \]
\[ Y_i = \beta_0 X_{1i}^{\beta_1} X_{2i}^{\beta_2} \exp\!\left[ \tfrac{1}{2}\beta_{11} (\ln X_{1i})^2 + \tfrac{1}{2}\beta_{22} (\ln X_{2i})^2 + \beta_{12} \ln X_{1i}\ln X_{2i} \right] \exp(V_i-U_i), \]
Compared estimators:
Linear SFM, GAM-SFM, BART-SFM, DNN-ReLU, DNN-ELU, DNN-CReLU, DNN-CELU.
| DGP | Best Parametric | Best Unconstrained DNN | Best Shape-Aware DNN |
|---|---|---|---|
| DGP1 (Cobb-Douglas) | Lin | DNN-ELU (RMSE=0.583) | DNN-CReLU (0.510) |
| DGP2 (+ irrelevant regressors) | Lin | DNN-ELU (1.025) | DNN-CReLU (0.394) |
| DGP3 (Translog) | Translog | DNN-ReLU (0.404) | DNN-CReLU (0.347) |
:::
Parametric models dominate only if true.
Shape-aware DNNs outperform unconstrained DNNs and other flexible estimators in almost any metric.
Better recovers inefficiency – lowest Bias\(_u\) and Bias\(_{TE}\).
More robust to irrelevant regressors.
More robust to misspecification risk – even a misspecified Translog performs worse than shape-aware DNNs.
Few instances of wrong-skewness.
Average runtime over 10 replications.
Trained in PyTorch using Adam optimization.
| Estimator | Runtime (sec) |
|---|---|
| Lin (misspecified) | 0.02 |
| Translog | 0.03 |
| GAM-SFM | 0.05 |
| BART-SFM | 78.00 |
| DNN-ReLU | 1.06 |
| DNN-ELU | 1.10 |
| DNN-CReLU | 1.20 |
| DNN-CELU | 5.52 |
:::
Accuracy gains come at modest computational cost
Data:
Models fitted: Linear SFM, BART-SFM, DNN-ELU, DNN-CELU
Linear SFM shows wrong-skew problem and estimates \(\widehat\sigma_u=0\), fails to produce efficiency scores in this example.
Flexible methods produce TE scores, only DNN-CELU produces monotone and concave partial dependence.
Efficiency estimates
Shape behavior
GitHub for this paper
with all Python codes and data
Doing SFA with Python: Adapting New Estimators to a High-Performance Computing Environment (with Yu Ma and Zheng Wei)
in Handbook on Research Methods and Applications in Productivity and Efficiency, Daraio, Cinzia (Ed.), 2026, Edward Edgar
iCEBDA26 in Xiamen, 4-7 Dec 2026
Keynotes: Mark Hallin (U Libre de Bruxelles), Peter Hansen (UNC Chapel Hill), Hashem Pesaran (U Combridge), Yixiao Sun (UC San Diego)