Under active development. A stable version is expected in November 2026. Feedback is welcome by email or via GitHub issues.
5 Feed-Forward Neural Networks
5.1 Overview
Feed-forward neural networks are best understood as flexible nonlinear regression models. In classical econometrics, we often begin with a parametric specification such as a linear model, a polynomial expansion, or a spline basis. A feed-forward neural network keeps the same basic goal, estimating a conditional mean or conditional class probability, but replaces curated basis functions by a learned composition of affine transformations and nonlinear activation functions.
For econometricians, the main attraction is flexibility. When the relationship between two random variables is nonlinear, high-dimensional, or shaped by complicated interactions, a neural network offers a way to approximate that relationship without fully specifying it in advance. This is useful in applications such as the nonlinear asset-pricing setting studied by Gu, Kelly, and Xiu (2020) and in other prediction problems with large information sets:
Asset pricing: mapping many firm characteristics into expected returns
Macroeconomic forecasting: combining large predictor sets, survey variables, or text-based indicators
Credit risk: estimating default probabilities from borrower characteristics and payment histories
Household or firm behavior: approximating nonlinear production relationships or treatment-response interactions
5.2 Roadmap
We begin with the simplest building block, a neuron or unit, and then build up feed-forward networks as flexible nonlinear basis expansions.
We then discuss activation functions, the transformations that make successive layers nonlinear, and an approximation theorem that explains what sufficiently wide networks can represent.
Next, we connect neural-network estimation to loss minimization, maximum likelihood, and gradient descent, before introducing backpropagation as an efficient way to reuse chain-rule derivatives across layers.
We then show that the training objective is non-convex, connect the fitted network to familiar econometric objects, and explain why its individual weights are not identified.
We apply a network with one hidden layer to credit-card default probabilities, compare it with logistic regression on untouched test data, and then close with further applications and a practical training workflow.
5.3 Basic Architecture
A neural network is a composition of affine maps and nonlinear activations, organized in layers with weights that adjust during learning. We use the standard econometric convention that Y and \mathbf{X}=(X_1,\ldots,X_p)^\top denote the random outcome and random predictor vector, lowercase y and \mathbf{x} denote their realized values, and a hat denotes a fitted value. The cleanest entry point is a single neuron, also called a unit. For a realized predictor vector \mathbf{x}=(x_1,\ldots,x_p)^\top, let \mathbf{w}=(w_1,\ldots,w_p)^\top be a vector of adjustable coefficients called weights, and let b be an adjustable intercept called the bias. The unit first forms the affine score z=\mathbf{x}^\top\mathbf{w}+b and then applies an activation functiong to that score. Its scalar prediction for Y at input \mathbf{x} is
We use unit and neuron interchangeably. Historically, perceptron refers more narrowly to a threshold-based linear classifier; a sigmoid neuron is a smooth generalization rather than the original perceptron model. If g(z)=z is the identity function, the unit defines a linear predictor. For binary classification, a common choice is the sigmoid function
\sigma(z)=\frac{1}{1+e^{-z}},
which maps any real-valued score z monotonically into the interval (0,1). We can interpret \hat p=\sigma(z) as a fitted probability and classify an observation as Class 1 when \hat p>0.5. Because \sigma(0)=0.5, this classification rule has decision boundary z=\mathbf{x}^\top\mathbf{w}+b=0. Under symmetric 0–1 misclassification loss, the 0.5 cutoff is Bayes-optimal when the output equals the conditional class probability \mathbb{P}(Y=1\mid\mathbf{X}=\mathbf{x}). With an asymmetric loss—the normal case when the two errors have different economic costs, as in default or recession prediction—the optimal threshold moves away from 0.5. Here, Bayes-optimal means that the decision rule minimizes the expected loss under the conditional distribution of the outcome; for a general treatment of how the loss function determines the optimal decision, including such classification cutoffs, see Gneiting (2011).
Figure 5.1 illustrates this rule for two predictors using \mathbf{w}=(-0.2,1)^\top and b=-1.4. The dashed line collects the predictor combinations for which the score is zero.
Figure 5.1: A sigmoid neuron with weights \mathbf{w}=(-0.2,1)^\top and bias b=-1.4. The dashed line is the zero-score decision boundary x_2 = 1.4 + 0.2x_1, where the fitted probability equals 0.5. Class 1 points lie above the boundary and Class 0 points below it.
In Figure 5.1, every Class 1 point has \mathbf{x}^\top\mathbf{w}+b>0 and hence a fitted probability above 0.5, while every Class 0 point has \mathbf{x}^\top\mathbf{w}+b<0 and hence a fitted probability below 0.5. The sigmoid transforms the score nonlinearly, but the decision boundary of this single unit remains a straight line. A feed-forward network simply stacks many such units, so that the outputs of one layer become regressors for the next layer. That is why neural networks can be understood as learned basis expansions rather than as a fundamentally alien model class. In a polynomial or spline regression, the transformed regressors are chosen before estimation; in a neural network, the hidden activations a_j^{(l)} (the activation of unit j in layer l, defined formally below) are themselves learned from the data.
A network is described in part by its depthD. Under the convention used here, it has D-1 hidden layers followed by an output layer, which we index D. The output layer need not contain a single unit. With several output units, the network returns one number per unit. A multi-class classifier—with one output unit for each possible class—is an example developed later in this chapter. Thus the network output is generally a vector, which we denote generically in boldface by \hat{\mathbf y}(\mathbf{x}). Stacking the layer maps gives the composition
This composition is also what the name feed-forward refers to: information flows in one direction only, from the inputs through the hidden layers to the output. Each layer feeds its activations forward to the next layer, no activation is ever fed back into an earlier layer, and the graph of connections contains no cycles—so the network output is a function of the current input alone. The recurrent networks of a later chapter drop exactly this restriction by feeding activations back across time steps.
When the network has one output unit, \hat{\mathbf y}(\mathbf{x}) reduces to the scalar \hat y(\mathbf{x}) introduced above. The layers play distinct roles: the input layer holds the realized predictors \mathbf{x}, the hidden layers compute intermediate activations, and the output layer produces the prediction \hat{\mathbf y}(\mathbf{x}). Concretely, each layer l computes:
\mathbf{a}^{(l)} is the output of layer l (with \mathbf{a}^{(0)} = \mathbf{x})
H^{(l)} is the number of units (the width) of layer l, with H^{(0)}=p the input dimension
\mathbf{W}^{(l)} is the weight matrix for layer l, with dimensions H^{(l)} \times H^{(l-1)}. Its entry W_{hj}^{(l)} connects neuron j in layer l-1 to neuron h in layer l. The first subscript therefore identifies the receiving unit and the second identifies the sending unit: W_{hj}^{(l)} is the weight on the connection j\to h. The single-neuron vector \mathbf w introduced above is the transpose of one row of \mathbf W^{(1)}. More generally, the incoming weight vector for unit h in layer l is the transpose of row h of \mathbf W^{(l)}, and its j-th coordinate is W_{hj}^{(l)}. This destination-first ordering follows from the matrix product \mathbf{W}^{(l)}\mathbf{a}^{(l-1)}.
\mathbf{b}^{(l)} is the bias vector for layer l
g^{(l)} is the activation function for layer l, applied per vector entry
\boldsymbol{\theta} := \big(\mathbf{W}^{(1)}, \mathbf{b}^{(1)}, \dots, \mathbf{W}^{(D)}, \mathbf{b}^{(D)}\big) collects every weight matrix and bias vector of the network, stacked into a single parameter vector. Whenever we write f(\mathbf{x};\boldsymbol{\theta}), differentiate with respect to \boldsymbol{\theta}, or update \boldsymbol{\theta} during estimation, we refer to this stacked vector of all weights and biases.
Moreover, we define the pre-activation value in layer l: \mathbf{z}^{(l)} = \mathbf{W}^{(l)}\mathbf{a}^{(l-1)} + \mathbf{b}^{(l)}.
Figure 5.2 makes the layered vector map concrete in a network with three output units.
Figure 5.2: Example of a feed-forward network with three input coordinates, two hidden layers, and three output units. The left column holds the inputs, the middle columns the hidden-layer activations a_i^{(l)}, and the right column the coordinates \hat y_1, \hat y_2, and \hat y_3 of the vector output \hat{\mathbf y}(\mathbf{x}). In a three-class classifier, one output unit is associated with each class. Each arrow is a weighted connection.
The sigmoid neuron of Figure 5.1 is a single unit. Figure 5.2 shows what changes in a multilayer network: instead of specifying a basis expansion by hand, we let the model learn many intermediate transformed regressors a_j^{(l)} and then combine them again in later layers. The three right-hand nodes make the vector notation explicit. They receive the same final hidden activations but use output-unit-specific weights and biases; a scalar-output network is the special case with only one such node.
Single-Hidden-Layer Networks
Much of this chapter works with single-hidden-layer networks, as do Exercises 5.1 and 5.2. These networks have the form:
w_j^{(1)}, b_j^{(1)} are the input-to-hidden weights and biases
w_j^{(2)}, b^{(2)} are the hidden-to-output weights and bias
g(\cdot) is the activation function
This scalar-input, scalar-output notation is a componentwise specialization of the layer matrices above: w_j^{(1)}=W_{j1}^{(1)}, w_j^{(2)}=W_{1j}^{(2)}, and b^{(2)}=b_1^{(2)}. For a predictor vector \mathbf x, the incoming vector \mathbf w_j^{(1)} is the transpose of row j of \mathbf W^{(1)}.
This display is written for a scalar predictor x, which keeps the algebra readable and is all the exercises need. For a predictor vector \mathbf{x} \in \mathbb{R}^p, nothing changes except that the scalar product w_j^{(1)} x becomes the inner product \mathbf{w}_j^{(1)\top}\mathbf{x} with \mathbf{w}_j^{(1)} \in \mathbb{R}^p; the parameter count below is stated for that general case.
Figure 5.3 shows this architecture in its vector-input form for p=4 predictors and H=5 hidden units. Reading the diagram left to right retraces the formula: every input coordinate feeds every hidden unit through an input-to-hidden weight, each hidden unit applies the activation g to its weighted sum plus bias, and the single output node combines the five hidden activations once more into the prediction. In contrast to the multi-output network of Figure 5.2, the regression case needs only this one output unit.
Figure 5.3: A single-hidden-layer feed-forward network with four inputs, five hidden units, and one output. The left column holds the input coordinates x_1,\dots,x_4, the middle column the hidden activations a_1^{(1)},\dots,a_5^{(1)}, and the right node the scalar output \hat y. Each arrow is a weighted connection; each hidden unit and the output unit additionally carry a bias, which the diagram does not draw.
Why Width and Depth Matter Quantitatively A useful discipline is to count parameters before discussing “big” or “small” networks. If the input dimension is p, a single-hidden-layer network with H hidden units and one scalar output contains
(p+1)H + (H+1) = H(p+2)+1
parameters: pH input-to-hidden weights, H hidden biases, H hidden-to-output weights, and one output bias. For the small network of Figure 5.3, with p=4 and H=5, this already gives 5 \cdot 6 + 1 = 31 parameters—of which only the 20 input-to-hidden weights and 5 hidden-to-output weights are drawn as arrows. Even moderate changes in width can therefore increase estimation error variance quickly when p is already large.
More generally, under this chapter’s convention that depth D counts the parameterized layers and excludes the input layer, a fully connected network with hidden widths H^{(1)}, \dots, H^{(D-1)} and a scalar output (so H^{(D)} = 1) has
parameters. Each layer contributes its weight matrix plus one bias per unit. (For D = 2 this reduces to the single-hidden-layer count H(p+2)+1 above.)
Question for Reflection
Suppose a macro forecasting problem has p=40 predictors and you compare one hidden layer with H=10 against one hidden layer with H=100. How many parameters does each network have, and what does that imply for overfitting risk?
Suggested Answer
With one scalar output, the count is H(p+2)+1. For H=10, this gives 10(42)+1=421 parameters. For H=100, it gives 100(42)+1=4201 parameters. The larger network is much more flexible, but it also has roughly ten times as many free parameters, so it can fit idiosyncratic sample noise much more easily unless the sample is large and validation supports the extra complexity.
5.4 Activation Functions and Universal Approximation
The activation function controls both the shape a unit can represent and the gradient signal available during estimation. Six choices recur in this and the following network chapters:
The sigmoid, introduced in Section 5.3, is \sigma(z)=1/(1+e^{-z}). Its range is (0,1) and its derivative is \sigma'(z)=\sigma(z)(1-\sigma(z)), which makes it a natural output activation for a binary class probability.
The hyperbolic tangent (tanh) is \tanh(z)=(e^z-e^{-z})/(e^z+e^{-z}). It maps into (-1,1) and satisfies \tanh'(z)=1-\tanh^2(z). Its zero-centered output is useful in the recurrent architectures introduced in the recurrent neural networks chapter, although its derivative also becomes small when the unit saturates—that is, when the pre-activation is far enough from zero that the activation sits in one of its flat regions, here near \pm 1, where the derivative is close to zero.
The rectified linear unit (ReLU) is \text{ReLU}(z)=\max(0,z). Its range is [0,\infty) and its derivative is 1 for z>0 and 0 for z<0; at the kink z=0, software must choose a convention, usually 0, as the callout below explains. ReLU is useful in hidden layers, but it is a poor output activation for non-negative regressand: once an output pre-activation becomes negative, its zero derivative supplies no signal that can move the unit back into the positive region.
The softplus function, \text{Softplus}(z)=\ln(1+e^z), is a smooth approximation to ReLU with range (0,\infty). Its derivative is strictly positive, which makes it a safer output activation when a continuous prediction must remain positive.
For a K-class outcome, the softmax maps an output-score vector into probabilities:
where \delta_{kq}=1 if k=q and 0 otherwise is the Kronecker delta, distinct from the backpropagation error signal introduced in Section 5.6.
The linear, or identity, activation is \text{Identity}(z)=z. Its unrestricted range makes it the standard output activation for an unbounded continuous regression target.
Figure 5.4 shows the four most common non-linear activation functions.
Show the code
import numpy as npimport matplotlib.pyplot as pltz = np.linspace(-5, 5, 200)# The style is applied inside a context manager so that it does not leak# into the styling of any later figure in this chapter.with plt.style.context('seaborn-v0_8-whitegrid'): fig, axes = plt.subplots(2, 2, figsize=(10, 8))# Sigmoid sigmoid =1/ (1+ np.exp(-z)) axes[0, 0].plot(z, sigmoid, color='blue') axes[0, 0].set_title(r'Sigmoid: $\frac{1}{1+e^{-z}}$') axes[0, 0].set_xlabel('z') axes[0, 0].set_ylabel('g(z)')# Tanh tanh = np.tanh(z) axes[0, 1].plot(z, tanh, color='green') axes[0, 1].set_title(r'Tanh: $\tanh(z)$') axes[0, 1].set_xlabel('z') axes[0, 1].set_ylabel('g(z)')# ReLU relu = np.maximum(0, z) axes[1, 0].plot(z, relu, color='red') axes[1, 0].set_title(r'ReLU: $\max(0, z)$') axes[1, 0].set_xlabel('z') axes[1, 0].set_ylabel('g(z)')# Softplus: np.logaddexp(0, z) is the overflow-safe form of log(1 + exp(z)),# which would overflow for large z if written literally. softplus = np.logaddexp(0, z) axes[1, 1].plot(z, softplus, color='purple') axes[1, 1].set_title(r'Softplus: $\ln(1+e^z)$') axes[1, 1].set_xlabel('z') axes[1, 1].set_ylabel('g(z)')for ax in axes.flat: ax.spines['top'].set_visible(False) ax.spines['right'].set_visible(False) ax.axhline(y=0, color='gray', linestyle='--', lw=1) ax.axvline(x=0, color='gray', linestyle='--', lw=1) plt.tight_layout() plt.show()
Figure 5.4: The four most common activation functions—Sigmoid, Tanh, ReLU, and Softplus—plotted against the pre-activation input z (x-axis) and their output g(z) (y-axis), arranged in a 2×2 grid. Dashed reference lines mark zero on both axes.
Each function provides a different kind of nonlinearity, from the bounded “S-shape” of the sigmoid and tanh to the unbounded, one-sided linearity of ReLU and its smooth counterpart, Softplus.
Additional detail: ReLU at the kink
ReLU is not differentiable at z=0. Implementations therefore assign a conventional value to \operatorname{ReLU}'(0), usually 0(Goodfellow, Bengio, and Courville 2016). We return to what this convention means for network gradients after introducing backpropagation below.
Practical selection. For a feed-forward hidden layer, ReLU is a transparent baseline: away from its kink, positive pre-activations pass their gradient through unchanged instead of multiplying it by a derivative below one. This mitigates vanishing gradients—gradient signals shrinking toward zero through saturated activations—relative to sigmoid or tanh units (Glorot, Bordes, and Bengio 2011). Tanh remains relevant when the architecture specifically benefits from bounded, zero-centered hidden states, as in the recurrent networks of the coming chapters.
For the output layer, the support and interpretation of the prediction determine the activation:
Regression problems (unbounded): Linear activation
Regression (positive): Softplus
Binary classification: Sigmoid for probability interpretation
Multi-class Classification: Softmax
Why Activation Functions Matter Without nonlinear activation functions, a deep neural network collapses into a single linear model. With two layers of weight matrices \mathbf W^{(1)} and \mathbf W^{(2)} and zero biases for exposition, the output of a linear-activation network is again a matrix times the input:
Both ReLU and softplus have non-negative range. Using the activation descriptions above, why is softplus rather than ReLU the recommended output activation for a positive regression outcome such as firm sales or stock-market volatility?
Suggested Answer
Both activations keep the prediction non-negative, so the difference lies in the gradient signal. Once a ReLU output unit’s pre-activation turns negative, its output is zero and its derivative is zero, so no gradient can move the unit back into the positive region. Softplus has a strictly positive derivative everywhere, so the output unit keeps receiving a gradient signal however far its pre-activation drifts. Softplus is therefore the safer output activation for a positive outcome.
Universal Approximation Theorem
Theorem (Hornik 1991, Theorem 2): Suppose g is a nonconstant, bounded, continuous activation function, such as a sigmoid. For any continuous function f on a compact (i.e., bounded and closed) set K \subset \mathbb{R} and any \varepsilon > 0, there exist a positive integer H and parameters \{w_j^{(1)}, b_j^{(1)}, w_j^{(2)}\}_{j=1}^{H} and b^{(2)} such that the function
The theorem is an existence result: it says a suitable width H exists, but does not tell us how to pick H or how to find the weights.
The notation here matches the single-hidden-layer representation introduced above. The same result extends from scalar x to vector-valued \mathbf{x} \in \mathbb{R}^p, with K \subset \mathbb{R}^p compact and w_j^{(1)} x replaced by \mathbf{w}_j^{(1)\top}\mathbf{x}.
Attribution. This result is often called the “Cybenko–Hornik” theorem, but the two papers do not prove the same statement, and the difference is the hypothesis on g. Cybenko (1989) proves the earlier and weaker version: his activations must be continuous and sigmoidal, meaning g(z) \to 1 as z \to +\infty and g(z) \to 0 as z \to -\infty, and his domain is the unit cube. The bounded-nonconstant-continuous form stated above is Hornik (1991).
The activation condition matters. The statement quoted here is the classical bounded-activation theorem, with sigmoids as the leading example. Modern networks often use ReLU, which is unbounded and therefore not covered by this version of the theorem. However, the extension of the theorem that does cover ReLU is Leshno et al. (1993).
Practical Implications
The theorem provides existence, not construction. It says that some width H suffices, but it does not say how to choose H, how to find the weights, or how easy the resulting optimization will be. The width required in practice can be very large, and the underlying non-convex optimization remains hard whether or not the approximation theorem applies.
The guarantee is also local in an important sense: the approximation holds only on the compact set K, and the theorem says nothing about the network’s behavior outside K. This matters in econometric applications, where predictors such as returns or growth rates have unbounded support: for inputs beyond the region the compact set covers—say, a crisis-period observation larger than anything seen before—the approximation bound is silent, and the network’s extrapolation there is whatever the fitted weights and the activation function happen to produce.
With a fixed architecture, the network specifies a finite-dimensional model class: the family of candidate regression functions traced out as the parameter vector containing all weights and biases ranges over its admissible values. Universal approximation establishes that rich functions can be approximated; consistency of the estimated network additionally requires the complexity to grow slowly enough relative to the available data (White 1990; Chen 2007).
5.5 Loss Functions and Estimation
We follow the optimization chapter in writing L(\boldsymbol{\theta}) for the sample-average loss over N observations, and \ell(y_i, f(\mathbf{x}_i;\boldsymbol{\theta})) for observation i’s contribution, which that chapter writes as \ell_i(\boldsymbol{\theta}). Throughout, \boldsymbol{\theta} is the parameter vector defined in Section 5.3—all weight matrices and bias vectors of the network stacked into one vector—so a gradient \nabla_{\boldsymbol{\theta}} L collects the partial derivatives of the loss with respect to every individual weight and bias.
Mean Squared Error Loss (Regression)
For regression problems, we typically use the mean squared error:
The factor of \frac{1}{2} is sometimes included for convenience when taking derivatives. For a fixed network architecture, minimizing this criterion is nonlinear least squares: the fitted values depend nonlinearly on \boldsymbol{\theta}, but the estimator is still defined as a minimizer of a sum of squared residuals. We can employ other loss functions as well, to which we return when discussing distributional networks.
Cross-Entropy Loss (Classification)
For binary classification problems, we typically use the cross-entropy (log-likelihood) loss. Such a network has a single output unit, and we write z_i^{(D)} for that unit’s pre-activation value at input \mathbf{x}_i—the real-valued score the network produces before the output sigmoid is applied, in the notation of Section 5.3. Keep this separate from f(\mathbf{x}_i;\boldsymbol{\theta}), which denotes the network’s output throughout this chapter: in a classification network the output is the fitted probability \hat p_i, not the discrete outcome scale 0 or 1. If y_i \in \{0, 1\} are the true class labels and \hat{p}_i=\sigma(z_i^{(D)}) are the predicted probabilities, where \sigma is the sigmoid function defined above, then
For multi-class classification, the observed class first has to be represented as a vector.
Definition: One-hot Encoding
With K possible classes, a one-hot vector has length K, contains a 1 in the position of the observed class, and contains 0s elsewhere. For example, if the classes are recession, normal growth, and boom, the three encodings are [1, 0, 0], [0, 1, 0], and [0, 0, 1]. The encoding records class membership without imposing an ordering on the class labels.
Using this representation, the categorical cross-entropy loss with a softmax output is
where y_{ik} is the one-hot indicator for class k and \hat{p}_{ik} = \frac{\exp(z_{ik}^{(D)})}{\sum_{q=1}^K \exp(z_{iq}^{(D)})} is the corresponding softmax probability. Here z_{ik}^{(D)} is the class-k output-layer pre-activation for observation i. As in the binary case, the softmax is applied to the scores, and the network’s output is the probability vector (\hat p_{i1},\dots,\hat p_{iK}).
Connection to Information Theory
The cross-entropy loss is not just a convenient formula—it has a direct information-theoretic foundation. Recall from the Information Theory chapter that minimizing cross-entropy between the empirical data distribution and the model is equivalent to maximum likelihood estimation, which in turn minimizes the Kullback–Leibler (KL) divergence between the two distributions. In this sense, training a neural network with cross-entropy loss is maximum likelihood estimation of a Bernoulli (or multinomial) model whose success probability is parameterized by the network.
Two conditions from that chapter carry over, and it is worth being precise about them, because the loose version of this statement is false. First, the sample average L(\boldsymbol{\theta}) = \frac{1}{N}\sum_{i=1}^N \ell_i(\boldsymbol{\theta}) is 1/N times a negative log-likelihood only if the random outcomes Y_1,\ldots,Y_N are conditionally independent given the predictors, so that the joint likelihood factorizes into the product of the per-observation Bernoulli terms. For dependent data the corresponding object is the prediction-error decomposition, with the conditioning set \mathcal{F}_{t-1} in place of the predictors of observation i. Second, the per-observation terms constitute a log-likelihood only if the Bernoulli (or multinomial) distribution is itself the assumed conditional model of the outcome, with the network parameterizing its probabilities. Without that full conditional specification, minimizing cross-entropy remains a sensible estimation criterion—it is a proper scoring rule—but it is not maximum likelihood.
Question for Reflection
Why is the cross-entropy loss for binary classification a special case of the categorical cross-entropy loss with K = 2?
Suggested Answer
With K = 2, the one-hot outcome vector has two entries, for example (y_i, 1-y_i), and the softmax probability vector is (\hat p_i, 1-\hat p_i). Plugging these into the categorical cross-entropy loss
which is exactly the binary cross-entropy loss. So the binary case is not a different principle; it is the two-class version of the general categorical formula.
Gradient Descent Optimization
The optimization chapter developed gradient descent in detail and writes its iterates as \boldsymbol\theta^{(t)}. To keep parenthesized superscripts reserved for layer indices, we instead use square brackets here; with learning rate \eta, the update is:
Unlike the quadratic ordinary least squares (OLS) objective, a neural-network loss is generally non-convex. The returned stationary point can therefore depend on the initialization and learning rate. Figure 5.5 provides a visual reminder of that difference before we turn to the network-specific task: computing \nabla_{\boldsymbol\theta}L efficiently.
(b) Non-convex: multiple local minima and dependence on initialization
Figure 5.5: Gradient-descent paths for two scalar objectives. Left: six updates from \theta=2.5 with learning rate 0.3 approach the unique minimum of L(\theta)=\theta^2. Right: seven updates with learning rate 0.25 from each of \theta=-1 and \theta=0.5 approach different minima of L(\theta)=\theta^4-\theta^2+0.2\theta, showing dependence on initialization.
Convex function (left): For the plotted quadratic and the step size used in the figure, every initialization converges to the unique global minimum, with step sizes naturally shrinking as the gradient norm decreases. (Convexity alone does not guarantee this: for the quadratic plotted here the iteration is x_{k+1}=(1-2\alpha)x_k, which converges only for 0<\alpha<1. The optimization chapter gives the general condition.)
Non-convex function (right): Different initializations can fall into different local minima; one of these can be worse than the global minimum.
Neural networks: Their highly non-convex loss landscapes make optimization sensitive to initialization and learning rate choices.
Gradient Descent in the Simplest Possible Network Consider a single-hidden-layer network with one neuron, of the following form:
f(x; \boldsymbol{\theta}) = w^{(2)} g(w^{(1)} x + b^{(1)}) + b^{(2)}
The parenthesized superscripts continue to index layers. When we need to index gradient-descent iterations, we use square brackets, as in \boldsymbol\theta[t].
The gradient descent procedure for a neural network can be visualized as follows:
flowchart TD
%% Starting point
START["θ[0] (Initial Parameters)"] --> FORWARD["Forward Pass: Compute f(x; θ[t])"]
%% Forward computation
FORWARD --> LOSS["Compute Loss L(θ[t])"]
%% Gradient computation
LOSS --> GRAD["Compute Gradients ∇L(θ[t])"]
%% Parameter update
GRAD --> UPDATE["Update: θ[t+1] = θ[t] - η∇L(θ[t])"]
%% Convergence check
UPDATE --> CHECK{"Converged?<br/>|∇L(θ[t+1])| < ε"}
CHECK -->|No| FORWARD
CHECK -->|Yes| OPTIMAL["θ̂ (Returned Parameters)"]
%% Styling
classDef startend fill:#e8f5e8,stroke:#333,stroke-width:2px
classDef process fill:#e1f5fe,stroke:#333,stroke-width:2px
classDef decision fill:#fff3e0,stroke:#333,stroke-width:2px
classDef update fill:#f3e5f5,stroke:#333,stroke-width:2px
class START,OPTIMAL startend
class FORWARD,LOSS,GRAD process
class CHECK decision
class UPDATE update
Figure 5.6: Flowchart of the gradient descent algorithm. Rectangular nodes represent computation steps (forward pass, loss evaluation, gradient computation, parameter update) and the diamond node is the convergence check. The loop continues until the gradient norm falls below a threshold \varepsilon, at which point the current parameters are returned. In a non-convex problem, a small gradient certifies only an approximate stationary point, not a global minimum.
Reading Figure 5.6 from top to bottom separates the two kinds of computation used in training: the forward pass produces predictions and a loss, and the gradient computation supplies the direction for the parameter update. The arrow back to the forward-pass node is the iterative optimization loop. The next section shows that the gradient computation is itself organized as a sweep through the network in the reverse direction—the backward pass of backpropagation, which is also why the prediction step carries the name forward pass. The terminal node should be read as “the algorithm stopped,” not “a global optimum was found,” because a small gradient in a non-convex network certifies only approximate stationarity.
5.6 Backpropagation
We now examine how the chain rule computes the gradient with respect to all parameters in a wide multilayer network with the objects defined in Section 5.3. Computing the gradient naively would require a separate chain-rule calculation for each weight and bias, which is inefficient. Backpropagation solves this problem by reusing intermediate quantities from the forward pass to compute all gradients in one backward sweep.
The algorithm was popularized by Rumelhart, Hinton, and Williams (1986). For a weight from unit j in layer l-1 to unit h in layer l, we want
differentiating the finite sum lets us derive the gradient for one observation and then average the resulting gradients; for example, the weight gradient is
For the one-observation derivation below, fix i and suppress its index: \ell, a_h^{(l)}, and z_h^{(l)} stand for \ell_i, a_{ih}^{(l)}, and z_{ih}^{(l)}, respectively. The layer index l remains explicit, and h and j label connected units. We treat gradients with respect to vectors as column vectors.
Error Signal Definition For each layer l, define the vector of error signals
This error signal measures how the loss changes when the pre-activation of unit h is perturbed, holding the other pre-activations in the layer fixed. Writing the pre-activation defined in Section 5.3 entry by entry gives
The weight W_{hj}^{(l)} affects z_h^{(l)} but no other pre-activation in layer l. Moreover, the activations a_r^{(l-1)} are computed from the input and the parameters of layers 1,\dots,l-1, so they do not vary with W_{hj}^{(l)}. Treating the vector \mathbf z^{(l)} as an intermediate quantity, the multivariate chain rule therefore gives
The (h,j)-th entry of the outer product is \delta_h^{(l)}a_j^{(l-1)}, which recovers the scalar weight-gradient formula.
Although the loss may depend on z_h^{(l)} through many downstream paths, all those paths are summarized by \delta_h^{(l)}. Once the error signals for layer l are available, the stored activation a_j^{(l-1)} supplies the remaining factor: each weight derivative requires one multiplication, while each bias derivative is obtained directly from the error signal.
Everything so far has therefore reduced the gradient computation to the error signals: if we knew every \delta_h^{(l)}, all weight and bias gradients would follow at the cost of one multiplication each. What we have not yet done is compute the error signals themselves. That is the task of the backward pass, and it is where the recursion runs in the reverse direction: the output layer’s signal follows directly from the loss, and each hidden layer’s signal is then obtained from the signals of the layer after it.
Backward Pass
For an element-wise output activation, the output-layer signal is
This formula assumes that a_h^{(D)} depends only on its own pre-activation, which holds for the element-wise activations introduced earlier. The softmax output is the notable exception; the foldable remark below shows how its error signal accounts for the normalization across coordinates.
Additional Detail: The Output Signal Under Softmax
Softmax violates the element-wise assumption because, through the normalization constraint across outputs indexed by k, every output coordinate depends on every output-layer pre-activation. Its error signal therefore sums over the full Jacobian of the output map,
With cross-entropy as the loss, the sum simplifies to \delta_k^{(D)} = \hat{p}_k-y_k: fitted class probability minus the one-hot indicator. Exercise 5.3 derives the analogous residual form for a sigmoid output and Bernoulli cross-entropy.
How can we reuse this information to calculate \boldsymbol{\delta}^{(l)} for hidden layersl = D-1,D-2,\ldots,1? Begin with the affine map into the next layer:
The next-layer error signal is \boldsymbol{\delta}^{(l+1)}=\nabla_{\mathbf z^{(l+1)}}\ell. The next-layer pre-activation \mathbf z^{(l+1)} is an affine function of \mathbf a^{(l)}, with Jacobian \mathbf W^{(l+1)}. Therefore, the chain rule for column-vector gradients gives
The transpose reverses the forward linear map: \mathbf W^{(l+1)} sends activations from layer l to layer l+1, whereas \mathbf W^{(l+1)\top} carries the gradient from layer l+1 back to layer l. To see how this operation determines the order of the subscripts, consider its h-th coordinate:
Thus W_{qh}^{(l+1)} appears because it is the (h,q)-th entry of the transposed matrix. Under the destination-first convention, it is the weight on the forward connection from unit h in layer l to unit q in layer l+1.
Finally, applying the chain rule through the element-wise activation at layer l gives the backpropagation recursion
\boldsymbol{\delta}^{(l)}
=
\underbrace{\left(
\mathbf W^{(l+1)\top}\boldsymbol{\delta}^{(l+1)}
\right)}_{\text{reused from layer }l+1}
\odot
\underbrace{g'^{(l)}(\mathbf z^{(l)})}_{\text{computed once at layer }l}.
Equivalently, its h-th coordinate is
\delta_h^{(l)}
=
\underbrace{\left(
\sum_{q=1}^{H^{(l+1)}}
W_{qh}^{(l+1)}\delta_q^{(l+1)}
\right)}_{\text{reused from layer }l+1}
\underbrace{g'^{(l)}(z_h^{(l)})}_{\text{computed once at layer }l}.
ReLU kinks and the computed gradient
Suppose the loss and all other activations are differentiable, so ReLU is the network’s only nonsmooth component. If no ReLU pre-activation for any training observation equals zero, the sample objective is differentiable at the current parameter value and the backpropagation recursion returns its ordinary gradient. If some pre-activation does equal zero, ReLU’s derivative does not exist. The composite objective may or may not be differentiable at that parameter value, but standard implementations still apply the chain rule using their assigned value for \operatorname{ReLU}'(0).
Question for Reflection
Suppose you ignored the recursion and instead computed each \frac{\partial \ell}{\partial W_{hj}^{(l)}} separately, applying the chain rule from the loss all the way down to layer l each time. Which quantities would you recompute, and which does the backward pass compute only once?
Suggested Answer
Computing the gradient for a weight in layer l from scratch requires propagating the derivative through every layer between the output and layer l—that is, re-deriving the entire chain \delta^{(D)}, \delta^{(D-1)}, \dots, \delta^{(l)}. Repeating this separately for parameters recomputes substantial portions of the same chain.
The recursion avoids this because \boldsymbol{\delta}^{(l)} is computed once from \boldsymbol{\delta}^{(l+1)}, and everything above layer l is already summarized in that one vector. The forward pass stored \mathbf{z}^{(l)} and \mathbf{a}^{(l-1)}. At layer l, the backward pass evaluates g'^{(l)} once at the stored \mathbf{z}^{(l)} and performs one matrix-vector product with \mathbf{W}^{(l+1)\top}. The gradients for all weights in the layer then follow from \boldsymbol{\delta}^{(l)} and \mathbf{a}^{(l-1)} without further propagation. Reverse-mode backpropagation therefore costs a small constant multiple of the forward pass. Both passes scale with the network’s active connections and arithmetic operations; the advantage is reuse, not independence from the number of parameters.
Full Gradient Update Step
In your optimization algorithm, you are currently at parameter value \boldsymbol{\theta}[t]. The two passes below run for one observation at a time; the resulting per-observation gradients are then averaged over the mini-batch before the parameter update, exactly as in the optimization chapter.
Part 1: Perform a Forward Pass of Your Data
For an observation with predictors \mathbf x, compute the network at the current parameters. Store the pre-activations and activations for each layer:
\mathbf{a}^{(l-1)}
\mathbf{z}^{(l)}
for l=1,\dots,D, together with the final output \mathbf{a}^{(D)}. The backward pass will evaluate each activation derivative from the stored \mathbf z^{(l)}.
Part 2: Reverse Iteration of Gradients
Compute the output-layer error signal \boldsymbol{\delta}^{(D)}, and immediately convert it into the output-layer gradients:
2.1 Calculate the error signal for layer l by combining the \boldsymbol{\delta}^{(l+1)} just produced in the previous iteration with the activation derivative at layer l:
(We write \mathcal{H} for the Hessian to avoid a clash with the layer width H.)
A twice-differentiable function on an open convex domain is convex if and only if its Hessian is positive semidefinite (PSD) at every point of that domain. Neural-network loss functions are nonlinear in the parameters (because of the activation function g), and the loss is in general not a convex function of \boldsymbol \theta. We establish this with a symmetry argument, stated for the one-unit network f(x) = w^{(2)}\, g(w^{(1)} x + b^{(1)}) + b^{(2)} of Section 5.5. The same symmetry then yields two consequences that are worth keeping apart, because they concern different things: non-convexity of the loss surface, an optimization statement about how hard the minimum is to find, and non-uniqueness of the minimizers, an identification statement about what the estimated parameters mean, developed formally under “Non-identification of \boldsymbol\theta” below.
The mirror symmetry. For the sigmoid, g(-z) = 1 - g(z), so the parameter vectors
implement exactly the same functionf(x) = w^{(2)}\, g(w^{(1)} x + b^{(1)}) + b^{(2)} and therefore attain exactly the same loss, L(\boldsymbol\theta) = L(\boldsymbol\theta_{\text{mir}}).
Consequence 1: the loss is not convex. Fix any parameter vector \boldsymbol\theta = (w^{(1)}, b^{(1)}, w^{(2)}, b^{(2)}) and pair it with its own mirror image \boldsymbol\theta_{\text{mir}}. Averaging the two coordinate by coordinate, the first three entries cancel and the midpoint is
Because the midpoint’s hidden-to-output weight is zero, it represents a constant function with value b^{(2)} + w^{(2)}/2. If L were convex, the defining inequality applied to this pair would give
where the final equality uses L(\boldsymbol\theta) = L(\boldsymbol\theta_{\text{mir}}) from the mirror symmetry. Since \boldsymbol\theta was arbitrary, convexity would force every network fit to be weakly beaten by a constant fit—namely, by the constant b^{(2)} + w^{(2)}/2 built from its own parameters. In particular, no parameter vector could attain strictly lower loss than every constant fit. Whenever the data are such that some network fit is strictly better than every constant fit—the typical case when the outcome varies systematically with the predictor—this is a contradiction, so L cannot be convex. Note that the argument never assumes a minimizer exists, which matters because on an unbounded parameter space the infimum need not be attained.
Consequence 2: minimizers are not unique. Whenever a minimizer does exist, it is not unique: if the optimal fit is non-constant then any minimizer has w^{(2)} \neq 0, and its mirror image is a different parameter vector with the same loss; and if the optimal fit is constant, then every (w^{(1)}, b^{(1)}, 0, c) with the same c attains it, which is a continuum. Unlike Consequence 1, this is not a statement about the difficulty of optimization—it is a statement about what the minimizing parameters mean. We return to it in the discussion of identification below.
A Hessian Diagnostic
The symmetry argument establishes non-convexity without calculating the Hessian. A second route is to find a parameter point at which the Hessian is not PSD. For a loss with several parameters, the sign of one cross-partial derivative is not enough: off-diagonal entries of a PSD matrix may be positive or negative. A negative determinant of a symmetric 2\times 2 principal Hessian block, by contrast, shows that the block is indefinite. Here, a principal block means the submatrix obtained by retaining the same parameter coordinates in the rows and columns. The full Hessian must then also have a direction of negative curvature. Exercise 5.2 applies this diagnostic to the interaction between the input-to-hidden weight w^{(1)} and the hidden-to-output weight w^{(2)}.
Connection to Econometric Concepts
From an econometric perspective it pays to keep four objects apart, because the claims made in the rest of this chapter attach to different ones and are not interchangeable.
where \Theta denotes the set of admissible parameter vectors. At an arbitrary \boldsymbol{\theta} the network is a parameterized candidate function and nothing more. Since \boldsymbol{\theta} is a free argument rather than a function of the sample, f(\cdot\,;\boldsymbol{\theta}) is not an estimator of anything.
The estimator is defined by the rule that lets the sample choose the parameter value. Write \hat{\boldsymbol{\theta}} for the value that this rule selects from the minimizers of the criterion below, treating any tie-breaking convention as part of the rule; because it is a function of the sample, it is a random object, and so is the fitted regression function
So \widehat m is an estimator of m in the ordinary sense, and the one curve it returns on the dataset in front of us is the corresponding estimate; claims about sampling variability are claims about the rule, not about that single curve. (We continue to write \hat y and \hat y_i for the network’s output at whatever parameter value is currently held—during optimization, that is not \hat{\boldsymbol{\theta}}—and reserve \widehat m for the fitted estimator.)
Correct specification means m \in \mathcal{M}_\Theta: the target is exactly representable by some parameter value in the class. It does not mean \widehat m = m. Even under correct specification the two differ in any finite sample, because \hat{\boldsymbol{\theta}} is computed from finite data and estimation error remains even though it may vanish in the limit N\rightarrow \infty. Under misspecification, m \notin \mathcal{M}_\Theta, and an approximation error is present on top of the estimation error.
The estimation rule selects a minimizer of the sample loss:
where \ell is the loss appropriate to the prediction problem, such as squared loss or cross-entropy. This is the same sample-average criterion introduced above, now written in terms of the full network f.
Non-identification of \boldsymbol{\theta}. The \arg\min above should be read as a set, not a point. Consequence 2 of the mirror symmetry already produced two parameter vectors with the same fitted function and loss. A layer with H>1 units adds a second, more general source of multiplicity: permuting the units together with their incoming and outgoing weights leaves the network’s output unchanged, generating up to H! equivalent parameter vectors (distinct except at parameter values with symmetries of their own). The map from parameters to functions is therefore many-to-one, and \boldsymbol{\theta} is not identified, even in population and with infinite data.
Three conclusions for econometric practice follow. First, an individual estimated weight has no structural interpretation, and a standard error attached to it is not interpretable the way an OLS coefficient’s standard error is. Second, the object that is invariant to these symmetries—the fitted prediction function \widehat m—is the economically relevant one, and evaluating it out of sample is how the rest of the chapter proceeds. (Invariance does not by itself make \widehat m unique: distinct functions can attain the same minimum sample loss, and pinning down the population regression function as a unique risk minimizer requires the usual uniqueness and support conditions.) Third, good out-of-sample prediction is evidence that \widehat m approximates the conditional mean under the distribution that generated the data—not evidence about what would happen to Y if a predictor were intervened on. Establishing such a causal interpretation requires a research design and identifying assumptions that connect variation in \mathbf X to interventions; the non-identification of individual network weights is a separate parameterization issue.
Preprocessing neural-network inputs
Before fitting a network, a useful default is to standardize continuous predictors with nonzero training-sample (equivalently called estimation sample) variance: subtract each predictor’s training-sample mean and divide by its training-sample standard deviation. Large differences in input scales can produce unequal first-layer gradient magnitudes and make a common learning rate difficult to choose. Standardization reduces this source of imbalance; the optimization chapter explains the step-size problem.
Standardization is a practical default rather than a mathematical requirement. With unrestricted first-layer weights and biases, an invertible affine rescaling can be absorbed into those parameters, preserving the functions the network can represent.
The preceding sections separate the fitted probability function from the particular weights used to represent it. We now make that distinction concrete with a feed-forward neural network (FNN) in a cross-sectional credit-risk application: the object we evaluate is each client’s predicted probability of default, not an individual network coefficient.
The Taiwan credit-card default data contain 30,000 clients (Yeh 2009). For each client, the predictors combine a granted credit limit, age and coded demographic characteristics with six months of repayment statuses, bill amounts, and payment amounts. The repayment histories run from April through September 2005, and the binary outcome y_i records whether client i defaulted in the following month. The original study uses the data to compare several statistical and machine-learning methods, including a neural network, with an explicit emphasis on probability forecasts rather than only hard classifications (Yeh and Lien 2009).
The default share is 6{,}636/30{,}000=0.2212. A classifier that always predicts no default therefore attains 77.88% accuracy without using a single predictor. That number is a useful warning: accuracy at a fixed threshold can look respectable even when a model has learned nothing about how risk varies across clients.
Data design. We remove the row identifier and treat the documented and undocumented category codes as labels rather than quantities with meaningful numerical distances. The continuous variables are standardized, while the demographic and monthly repayment-status variables are one-hot encoded with one reference category omitted. Every transformation is estimated on the data available for fitting and then applied unchanged to held-out observations.
We divide the observations randomly and stratify by the outcome, meaning that the default share is kept approximately constant across a 60% training sample, a 20% validation sample, and a 20% test sample. The validation sample chooses among four architectures that follow a simple pyramid rule: each successive hidden layer is narrower as the network approaches its scalar output. The candidates have hidden-layer widths 8, 16\to8, 32\to16\to8, and 64\to32\to16\to8. The first is the smallest pyramid, contracting directly from eight hidden units to the output. We use no weight penalty, so architecture is the only tuned dimension. All candidates use the same rule that stops training when the training loss no longer improves by a stated tolerance. The test sample remains untouched until the architecture has been chosen. The comparison model is an unpenalized logistic regression fitted to the same transformed predictors. After selection, both models are refitted on the combined training and validation observations, which we call the development sample, before the final test evaluation.
This randomized design has a narrow interpretation. All clients share the same calendar period, so the test results estimate within-cohort cross-sectional performance under an exchangeability approximation. They do not estimate how a model trained in 2005 would perform on borrowers from a later credit regime. The data contain no later target period with which to conduct that out-of-time evaluation.
The data also contain sex, education, and marital-status codes. Including these historical fields in a benchmark prediction exercise is not an endorsement of their use in an operational lending decision. The score comparison below answers whether the stated procedures predict this sample’s outcome; it does not establish that the predictors are appropriate, that the resulting decisions would be fair, or that the model should be deployed.
We evaluate the probability forecasts using the cross-entropy loss defined in Section 5.5 and the Brier score
The Brier score is the mean squared error of a binary probability forecast and is a strictly proper scoring rule: in population, its expected value is uniquely minimized by reporting the true conditional event probability (Brier 1950; Gneiting and Raftery 2007). Lower cross-entropy and Brier scores are better. We also report accuracy after thresholding the probabilities at 0.5, but this last column evaluates the resulting decisions rather than the probability forecasts themselves.
Table 5.1: Validation log loss for four unpenalized pyramidal FNN architectures fitted to the Taiwan credit-card default training sample, each using the fixed random seed 42. Arrows separate successive hidden layers. Lower values are better.
Validation log loss
Architecture
8
0.4382
16 → 8
0.4608
32 → 16 → 8
0.5565
64 → 32 → 16 → 8
1.0123
On this split and with the single fixed initialization used for each candidate, validation log loss selects the one-hidden-layer architecture with eight units. Adding unpenalized layers makes validation performance worse for these fitted candidates, but the comparison is illustrative rather than an initialization-robust architecture ranking. The result illustrates why the pyramid rule defines a candidate set, not a guarantee that the deepest member should be used.
Table 5.2: Probability-forecast performance on the 6,000-client test sample from the Taiwan credit-card default data. The constant benchmark reports the development-sample default share for every client. Logistic regression and the selected FNN are refitted on the combined training and validation samples. Lower log loss and Brier score are better; higher accuracy under the 0.5 decision threshold is better.
Log loss
Brier score
Accuracy at 0.5
Model
Constant probability
0.5284
0.1723
0.7788
Logistic regression
0.4431
0.1390
0.8170
FNN
0.4369
0.1376
0.8170
On the test sample, the selected network lowers log loss from 0.4431 for logistic regression to 0.4369 and lowers the Brier score from 0.1390 to 0.1376. The gain is modest, while both procedures classify 81.70% correctly under the 0.5 threshold. There is no contradiction: the first two criteria evaluate the entire probability forecast, while accuracy discards probability magnitudes after applying a decision threshold.
A grouped calibration diagnostic compares predicted probabilities with observed event frequencies. For each model separately, we order its 6,000 test-set predicted probabilities from smallest to largest and split the clients into ten groups of roughly 600: the first group collects the clients with the lowest predicted default probabilities, the last group those with the highest. The groups are only approximately equal-sized because clients with tied predicted probabilities are kept in the same group. Because the grouping uses each model’s own predictions, the two models’ groups generally contain different clients. A point’s horizontal coordinate is a group’s mean predicted probability, and its vertical coordinate is that group’s realized default frequency. A perfectly calibrated set of grouped forecasts would lie on the 45-degree line, although agreement in ten groups would not by itself establish calibration conditional on the full predictor vector.
Figure 5.7: Grouped calibration of logistic-regression and FNN default probabilities on the 6,000-client test sample. Each point summarizes one model-specific probability decile: the horizontal coordinate is the decile’s mean predicted default probability and the vertical coordinate is its observed default frequency. The dashed 45-degree line marks equality of the two quantities. Deciles are formed separately for each model and contain equal numbers of clients up to ties.
Figure 5.7 shows that both models track the diagonal reasonably closely at low and intermediate predicted risks. Both also overpredict in their highest-risk decile: the mean forecast is about 0.73, while the realized frequency is about 0.67 for logistic regression and 0.68 for the network. The plot therefore adds information that the average scores suppress, but it remains a descriptive test-set diagnostic rather than a proof of conditional calibration.
The comparison is intentionally limited. The gradient-boosting chapter later compares this selected network with a single tree, a random forest, and gradient boosting on the same test clients.
Question for Reflection
On the fixed test sample, why can the neural network have lower cross-entropy and Brier scores than logistic regression while having exactly the same accuracy under the 0.5 threshold?
Suggested Answer
Cross-entropy and the Brier score use every reported probability and reward probabilities that move toward the realized outcomes. Accuracy first converts each probability into a class label at 0.5 and then ignores its magnitude. The network can therefore improve many probability estimates without changing the total number classified correctly.
5.9 Applications in Economics and Finance
The credit-default illustration is one cross-sectional application. Feed-forward networks also appear in asset-pricing and macroeconomic prediction problems with larger or time-indexed information sets.
Asset pricing and portfolio management. Predicting stock returns from a large set of firm characteristics and macroeconomic variables is the flagship application. The traditional route specifies linear factor models with pre-specified factors; feed-forward networks instead capture nonlinear interactions between characteristics and let the data suggest the functional form. Gu, Kelly, and Xiu (2020) show that flexible nonlinear specifications can outperform linear benchmarks for individual-stock return prediction.
Macroeconomic forecasting and nowcasting. For predicting gross domestic product (GDP) growth, inflation, or unemployment from high-dimensional datasets, the traditional tools are vector autoregressions (VARs) and dynamic factor models. Feed-forward networks handle many predictors and nonlinear interactions naturally, especially when the problem is cast as one-step-ahead prediction with a large information set. Richardson, van Florenstein Mulder, and Vehbi (2021) nowcast New Zealand GDP growth in real time from roughly 600 domestic and international macroeconomic and financial series, using genuine data vintages so that each nowcast uses only information available at the time. They find that machine-learning algorithms—among them a small feed-forward network with a single sigmoid hidden layer—improve on both an autoregressive benchmark and a dynamic factor model. The real-time design matters econometrically: macroeconomic data are revised, so evaluating on final-vintage data would credit the model with information no forecaster had.
5.10 Practical Training Workflow
The complexity of the optimization problem, together with large modern datasets, explains why practitioners rarely rely on plain batch gradient descent alone.
Batch gradient descent uses the entire training set for each update. Mini-batch gradient descent estimates the gradient from a subset, while stochastic gradient descent uses one randomly sampled observation. Modern training usually uses mini-batches, often with an adaptive optimizer such as Adam, RMSprop, or AdaGrad. The optimization chapter compares these update rules and explains their tuning parameters.
Step-by-step Guide for Coming Up With Good Network Architecture
1. Data Preprocessing
Split into training, validation, and test sets
For time-series applications, respect the information set available at the forecast date rather than shuffling observations at random
Standardize continuous predictors using training-sample statistics and apply the same transformation to validation and test data; see the input preprocessing discussion.
2. Architecture Design
Start with one or two hidden layers, which are often enough for tabular econometric data
A pyramidal candidate makes successive hidden layers narrower toward the output, but its depth still has to be justified by validation performance
Increase width or depth only if validation performance justifies the extra complexity. Architecture is itself a tuning dimension: Christensen, Siggaard, and Veliyev (2023) compare networks of one to four hidden layers in a realized-variance forecasting exercise, and find that the ranking across depths varies while all of them beat the heterogeneous autoregressive (HAR) benchmark. Their broader finding is a useful corrective in the other direction too—even a shallow network, implemented with minimal tuning, was competitive, so extra depth is not where the gains in that application came from. Gu, Kelly, and Xiu (2020) report the sharper version of the same lesson for return prediction: performance there deteriorates beyond roughly three hidden layers
Choose an appropriate output activation for the outcome type, using the guidelines in Section 5.4
3. Training Process
Use Adam or another mini-batch optimizer as a practical baseline; the optimization chapter explains these update rules and their tuning parameters
Monitor both training and validation loss
Treat the learning rate and batch size as tuning choices rather than fixed defaults
4. Iteration and Refinement
Overfitting: If training loss keeps falling while validation loss begins to rise, simplify the architecture or stop at the best validation iteration
Underfitting: If both training and validation losses remain poor relative to a simple benchmark, increase model capacity or continue training
5.11 Summary
Key Takeaways
Feed-forward networks compose affine transformations and nonlinear activations into flexible prediction functions.
The output activation and loss determine whether a network predicts a mean, a binary probability, or a vector of class probabilities.
Backpropagation reuses chain-rule calculations to obtain gradients for all network parameters.
Probability forecasts and thresholded classifications answer different questions and require matching evaluation criteria.
Common Pitfalls
Universal approximation does not choose an architecture or guarantee successful optimization and generalization.
Individual network weights are not identified structural coefficients or causal effects.
Treat category codes as labels unless their numerical ordering has a documented meaning.
A randomized same-period client split does not test forecasting performance in later credit regimes.
5.12 Exercises
Exercise 5.1: ReLU Networks as Adaptive Piecewise Linear Regressions
where the knots satisfy c_1 < c_2 < \cdots < c_K. Adopt the conventions c_0 = -\infty and c_{K+1} = +\infty, so that the intervals (c_r, c_{r+1}) for r = 0, 1, \dots, K cover the real line except for the knots c_1,\dots,c_K.
Show that f(x) is piecewise linear. More precisely, show that on each open interval (c_r,c_{r+1}), r = 0,1,\dots,K, its slope is
\alpha_1+\sum_{j=1}^r \gamma_j,
with the empty sum for r=0 equal to zero, so that the slope on (-\infty,c_1) is \alpha_1.
Show that this function can be written as a one-hidden-layer neural network with ReLU activation using at most K+2 hidden units, that is, in the form
f(x)=b^{(2)}+\sum_{m=1}^{K+2} w_m^{(2)}\,\mathrm{ReLU}(w_m^{(1)}x+b_m^{(1)}),
by choosing suitable values for b^{(2)}, w_m^{(1)}, b_m^{(1)}, and w_m^{(2)}. Here m indexes hidden units, whereas j continues to index the K knots. (The count is “at most” rather than “exactly” because units with w_m^{(2)}=0 contribute nothing and may be dropped, which happens for the two linear units when \alpha_1=0.)
Exam level. The exercise gives a formal representation result: every linear spline in the stated fixed-knot hinge form can be represented exactly by a one-hidden-layer ReLU network.
Hint for Part 2
Use \mathrm{ReLU}(x-c_j)=(x-c_j)_+. Also note that
The convention c_0=-\infty makes r=0 the leftmost piece: before the first knot no ReLU term is active, the sum is empty, and the slope on (-\infty,c_1) is simply \alpha_1. The convention c_{K+1}=+\infty makes r=K the rightmost piece, on which all K terms are active. Therefore f is piecewise linear, with slope changes only at the knot locations.
Part 2: Network Representation
We construct the network one part of f(x) at a time. Write the activation of hidden unit m as
The construction uses two hidden units for the linear term and one hidden unit for each of the K knots, hence at most K+2 hidden units.
Exercise 5.2: Gradient Descent and Non-Convexity in a One-Unit Network
Consider the problem of estimating a neural network for a regression problem with economic interpretation. You have data \{(x_i, y_i)\}_{i=1}^N where x_i represents log income and y_i represents log consumption.
You want to estimate a single-hidden-layer network:
f(x; \boldsymbol{\theta}) = w^{(2)} g(w^{(1)} x + b^{(1)}) + b^{(2)}
where g(z) = \frac{1}{1 + e^{-z}} is the sigmoid function and \boldsymbol{\theta} = (w^{(1)}, b^{(1)}, w^{(2)}, b^{(2)}), ordered by layer to match Section 5.7. Throughout this exercise, parenthesized superscripts index layers, the subscript i indexes observations, and square brackets such as [t] index gradient-descent iterations.
For every observation i=1,\ldots,N, define the pre-activation, hidden activation, fitted value, and residual by
flowchart LR
%% Network nodes
X["Input<br/>x"] -->|"w⁽¹⁾, b⁽¹⁾"| H["Hidden unit<br/>a"]
H -->|"w⁽²⁾, b⁽²⁾"| Y["Output<br/>ŷ"]
%% Styling
classDef input fill:#e1f5fe
classDef hiddenlayer fill:none,stroke:#333,stroke-width:2px
classDef output fill:#e8f5e8
class X input
class H hiddenlayer
class Y output
Figure 5.8: The single-hidden-layer network of Exercise 5.2. Log income x enters one sigmoid hidden unit a = g(w^{(1)} x + b^{(1)}), whose output feeds the linear output \hat{y} = w^{(2)} a + b^{(2)}. The arrow labels identify the parameters acting on each connection.
Derive the sample gradients \frac{\partial L}{\partial w^{(2)}}, \frac{\partial L}{\partial b^{(2)}}, \frac{\partial L}{\partial w^{(1)}}, and \frac{\partial L}{\partial b^{(1)}} for
L(\boldsymbol{\theta}) = \frac{1}{2N} \sum_{i=1}^N (y_i - f(x_i; \boldsymbol{\theta}))^2.
Report the chain-rule factor \frac{\partial \hat y_i}{\partial \theta_j} separately for each parameter before assembling the gradient, and identify the activation-derivative factor that appears in both hidden-layer gradients but neither output-layer gradient. Then state the gradient descent update for the full parameter vector \boldsymbol{\theta} in one line, with learning rate \eta>0. (The gradients themselves are stated in Section 5.5; derive them without looking, since the remaining parts build directly on these expressions.)
Establish non-convexity directly from the Hessian. Consider the one-observation sample N=1 with x_1=y_1=1, and evaluate the loss at
\boldsymbol{\theta}_0=(w^{(1)},b^{(1)},w^{(2)},b^{(2)})=(0,0,0,0).
Starting from the gradients in Part 1, compute the 2\times 2 principal Hessian block corresponding to (w^{(1)},w^{(2)}),
\mathcal H_{\{w^{(1)},w^{(2)}\}}(\boldsymbol{\theta}_0)
=
\begin{pmatrix}
\dfrac{\partial^2L}{\partial (w^{(1)})^2} & \dfrac{\partial^2L}{\partial w^{(1)}\partial w^{(2)}}\\[4pt]
\dfrac{\partial^2L}{\partial w^{(2)}\partial w^{(1)}} & \dfrac{\partial^2L}{\partial (w^{(2)})^2}
\end{pmatrix}_{\boldsymbol{\theta}=\boldsymbol{\theta}_0}.
Show that its determinant is negative, and explain why this proves that the full loss L(\boldsymbol{\theta}) is non-convex. Finally, explain why a negative cross-partial derivative by itself would not be sufficient.
Using g'(z)=g(z)(1-g(z)), show that
0<g'(z)\le \frac{1}{4}
for all z. To quantify saturation, suppose that for some M\ge 0,
|z_i|\ge M
\qquad\text{for every }i=1,\ldots,N.
Show that
g'(z_i)
\le
\varepsilon_M
:=
\frac{e^{-M}}{(1+e^{-M})^2}
\qquad\text{for every }i=1,\ldots,N.
Verify that M=0 recovers the global bound g'(0)=1/4, and show that \varepsilon_M\to0 as M\to\infty.
Combine the bound from Part 3 with the gradients from Part 1 to deduce
\left|\frac{\partial L}{\partial w^{(1)}}\right|
\le
\frac{|w^{(2)}|\varepsilon_M}{N}\sum_{i=1}^N |r_i x_i|,
\qquad
\left|\frac{\partial L}{\partial b^{(1)}}\right|
\le
\frac{|w^{(2)}|\varepsilon_M}{N}\sum_{i=1}^N |r_i|.
Explain why, for a fixed learning rate and when the other factors in the bounds remain bounded, stronger saturation can make the gradient-descent updates to w^{(1)} and b^{(1)} arbitrarily small. Finally, explain why saturation alone does not imply that the output-layer updates to w^{(2)} and b^{(2)} must also be small.
Exam level. The exercise derives all gradients of the one-unit network, uses a principal Hessian block to prove non-convexity, bounds the sigmoid derivative under saturation, and then shows how that bound attenuates the hidden-layer updates without necessarily attenuating the output-layer updates.
Hint for Part 1
Use the chain rule with the residual notation defined above:
\frac{\partial L}{\partial \theta_j}
=
\frac{1}{N}\sum_{i=1}^N r_i\frac{\partial r_i}{\partial \theta_j},
\qquad
\frac{\partial r_i}{\partial \theta_j}
=
-\frac{\partial\hat y_i}{\partial\theta_j}.
Remember that g'(z) = g(z)(1 - g(z))
Hint for Part 2
At \boldsymbol{\theta}_0, first calculate a_1, r_1, and \partial a_1/\partial w^{(1)}. Differentiate the Part 1 gradients once more to obtain the three distinct entries of the symmetric block.
A symmetric 2\times 2 matrix with a negative determinant has one positive and one negative eigenvalue. To see what this implies for the full Hessian, extend a direction in the (w^{(1)},w^{(2)}) coordinates by setting its b^{(1)} and b^{(2)} coordinates equal to zero.
Hint for Part 3
Rewrite g'(z) so that z enters only through e^z+e^{-z}, and study how this sum depends on |z|.
Solution
Part 1: Gradients and Updates
Using the observation-level notation defined in the exercise, the loss is
where the last two use \partial a_i/\partial z_i = g'(z_i)=a_i(1-a_i) together with \partial z_i/\partial w^{(1)} = x_i and \partial z_i/\partial b^{(1)} = 1. Substituting these factors gives the gradients.
The distinguishing activation-derivative factor is the sigmoid derivative
g'(z_i)=a_i(1-a_i),
which enters both hidden-layer gradients and neither output-layer gradient. The sigmoid derivative is absent at the output layer because w^{(2)} and b^{(2)} enter the fitted value linearly once the hidden activation a_i is given.
Part 2: A Hessian Witness of Non-Convexity
For N=1, x_1=y_1=1, and \boldsymbol{\theta}_0=(0,0,0,0), the hidden activation, fitted value, and residual are
A symmetric 2\times2 matrix with a negative determinant has one positive and one negative eigenvalue. If \mathbf v=(v_1,v_2)^\top is a corresponding direction of negative curvature in the (w^{(1)},w^{(2)}) coordinates, then \widetilde{\mathbf v}=(v_1,0,v_2,0)^\top is a direction in the full (w^{(1)},b^{(1)},w^{(2)},b^{(2)}) parameter space with a negative Hessian quadratic form. The full Hessian is therefore not PSD at \boldsymbol{\theta}_0. Because the Hessian of a twice-differentiable convex function must be PSD everywhere, L(\boldsymbol{\theta}) is non-convex.
The determinant is essential to the argument. A negative cross-partial derivative is only a negative off-diagonal Hessian entry, and PSD matrices may have negative off-diagonal entries. The negative determinant shows that this particular block is indefinite and hence supplies the required direction of negative curvature.
Part 3: Bounding the Sigmoid Derivative
The sigmoid derivative is
g'(z)=g(z)(1-g(z)),
and g(z)\in(0,1) for all z, so g'(z)>0. For u\in(0,1), the quadratic u(1-u) is maximized at u=1/2, where it equals 1/4. Therefore
This expression is symmetric around zero and decreases as |z| increases, because e^z+e^{-z} is symmetric and increasing in |z|. Consequently, if |z_i|\ge M for every i=1,\ldots,N, then
g'(z_i)
\le
\frac{1}{2+e^M+e^{-M}}
=
\frac{e^{-M}}{(1+e^{-M})^2}
=
\varepsilon_M
\qquad\text{for every }i=1,\ldots,N.
Setting M=0 gives \varepsilon_0=1/(2+1+1)=1/4, the global bound from the first display. Moreover, the numerator e^{-M} tends to zero while the denominator tends to one, so
\lim_{M\to\infty}\varepsilon_M=0.
Part 4: Gradient Bounds Under Saturation
Substituting the bound from Part 3 into the hidden-layer gradients from Part 1 and applying the triangle inequality gives
These inequalities condition on the current output weight, residuals, inputs, and data set. If the factors multiplying \varepsilon_M remain bounded and the learning rate \eta is fixed, increasing M makes the gradient-descent updates -\eta\,\partial L/\partial w^{(1)} and -\eta\,\partial L/\partial b^{(1)} arbitrarily small. This is the precise sense in which saturation can slow learning of the hidden-layer parameters.
The same conclusion does not follow for the output-layer parameters because their gradients are
which contain no sigmoid-derivative factor. Saturation therefore imposes no generic attenuation through g'(z_i) on either output-layer gradient. The residuals r_i still depend on the saturated activations, so these gradients are not literally unaffected, but nothing in their form forces them toward zero. The w^{(2)} gradient contains the activation a_i itself rather than g'(z_i): if a_i\approx1 for every i, it approximately equals the b^{(2)} gradient, whereas if a_i\approx0 for every i and the residuals remain bounded, the factor a_i attenuates it. Saturation may therefore suppress the hidden-layer updates even while the output-layer parameters continue to move.
Exercise 5.3: Cross-Entropy, Likelihood, and the Sigmoid Output Layer
Consider the observed sample \{(x_i,y_i)\}_{i=1}^N, where x_i is the network input and y_i\in\{0,1\} is the observed realization of the binary random outcome Y_i. Let \boldsymbol{\theta} collect all network weights and biases, and write z_i(\boldsymbol{\theta}) for the network’s output-layer pre-activation at input x_i; that is, the real-valued score before the sigmoid. The model-implied success probability is
which is the quantity that Section 5.5 writes as \hat p_i, with z_i(\boldsymbol{\theta}) playing the role of z_i^{(D)} there; we write \pi_i(\boldsymbol{\theta}) here to keep its dependence on the parameter vector explicit. Below, we work at a generic \boldsymbol{\theta} and suppress its argument, writing z_i and \pi_i.
Write x_{1:N}=(x_1,\ldots,x_N) for the full collection of predictors. Conditioning on x_{1:N} or on x_i below abbreviates conditioning on the event that the corresponding random predictors take these observed values. Assume throughout that, conditional on x_{1:N}, the random outcomes Y_1,\ldots,Y_N are independent and the conditional model satisfies
\mathbb{P}_{\boldsymbol{\theta}}(Y_i=1\mid x_{1:N})
=
\mathbb{P}_{\boldsymbol{\theta}}(Y_i=1\mid x_i)
=
\pi_i
\qquad\text{for every }i=1,\ldots,N.
Show that the negative Bernoulli log-likelihood contribution of observation i is
\ell_i(\boldsymbol{\theta})
=
-\left[y_i \log \pi_i + (1-y_i)\log(1-\pi_i)\right],
and explain where the conditional independence assumption is used when you combine these contributions into the sample-average loss L(\boldsymbol{\theta})=\frac{1}{N}\sum_{i=1}^N \ell_i(\boldsymbol{\theta}).
Viewing L as a function of the intermediate pre-activations z_{1:N}=(z_1,\ldots,z_N), show that
\frac{\partial \ell_i}{\partial z_i}=\pi_i-y_i,
and conclude that \frac{\partial L}{\partial z_i}=\frac{\pi_i-y_i}{N}.
Keep the same sigmoid output layer but replace cross-entropy by the observation contribution to the mean squared error (MSE) criterion,
\ell_i^{\mathrm{MSE}}
=
\frac12(y_i-\pi_i)^2.
Derive \frac{\partial \ell_i^{\mathrm{MSE}}}{\partial z_i}. Then compare the two derivatives for a confidently wrong prediction, y_i=1 with \pi_i\to0. Explain why cross-entropy avoids an additional attenuation from the output sigmoid in this case. Does this comparison establish that gradients for all earlier-layer parameters remain bounded away from zero? Give a reason for your answer.
Exam level. The exercise links binary classification in neural networks back to likelihood theory, derives the output-layer error signal used in backpropagation, and then shows exactly which source of gradient attenuation cross-entropy removes—and which sources may remain in earlier layers. Part 3 reuses the saturation mechanism of Exercise 5.2.
Hint for Part 2
Use the chain rule
\frac{\partial \ell_i}{\partial z_i}
=
\frac{\partial \ell_i}{\partial \pi_i}\frac{\partial \pi_i}{\partial z_i}.
Remember that \sigma'(z_i)=\pi_i(1-\pi_i).
Hint for Part 3
The chain rule runs through \pi_i in exactly the same way as in Part 2; only \partial \ell_i/\partial \pi_i changes.
Once you have both derivatives, evaluate each in the limit \pi_i\to0 with y_i=1 and compare the two limits.
For the last question, write out the full chain rule from \ell_i to a hidden-layer parameter and check which of its factors the comparison above has actually examined.
Solution
Part 1: From Bernoulli Likelihood to Cross-Entropy
Under the maintained conditional model, the probability mass assigned to the observed outcome y_i is \pi_i^{y_i}(1-\pi_i)^{1-y_i}. This observation-level step uses the restriction that conditioning on the other predictors does not change observation i’s success probability; it does not yet use independence across outcomes. Taking logs and negating gives
Conditional independence licenses the next step. It lets the joint likelihood conditional on x_{1:N} factorize into a product of the individual Bernoulli terms,
\prod_{i=1}^N \pi_i^{y_i}(1-\pi_i)^{1-y_i},
so that its negative log is the sum \sum_{i=1}^N \ell_i(\boldsymbol{\theta}), and the sample-average loss
is exactly 1/N times the negative log-likelihood—and is exactly the binary cross-entropy loss displayed in the chapter. Without the full conditional specification—both the Bernoulli probabilities and their factorization—the sum can still be used as a prediction loss or a separable estimation criterion, but it is not generally the joint negative log-likelihood. This matters for binary outcomes such as recession indicators or the sign of an asset return, which may be serially dependent; a likelihood can instead be built through the prediction-error decomposition, conditioning each term on the relevant history.
Part 2: Derivative with Respect to the Pre-Activation
Hence \begin{align*}
\frac{\partial \ell_i}{\partial z_i}
&=
\left(-\frac{y_i}{\pi_i}+\frac{1-y_i}{1-\pi_i}\right)\pi_i(1-\pi_i) \\
&=
-y_i(1-\pi_i)+(1-y_i)\pi_i \\
&=
\pi_i-y_i.
\end{align*} Viewing L as a function of z_{1:N}, the pre-activation z_i enters only through \ell_i. It follows that
The two output-layer error signals therefore differ by exactly the factor \pi_i(1-\pi_i)=\sigma'(z_i), which cancels from the cross-entropy derivative but remains in the squared-error derivative.
Now take the confidently wrong case y_i=1 with \pi_i\to0, equivalently z_i\to-\infty:
Thus the cross-entropy error signal with respect to the output pre-activation remains of order one, whereas squared error introduces an additional sigmoid-derivative factor that drives this signal to zero. Cross-entropy therefore avoids this particular source of gradient attenuation for confidently wrong predictions.
The comparison does not establish that gradients for all earlier-layer parameters remain bounded away from zero. For any network parameter \theta_j, the observation-level chain rule continues with
Even under cross-entropy, the upstream factor \partial z_i/\partial\theta_j may be small because earlier hidden activations are saturated or because many small derivatives are multiplied across layers. Cross-entropy removes the extra attenuation from the output sigmoid; it does not eliminate every source of vanishing gradients in the network.
Chen, Xiaohong. 2007. “Large Sample Sieve Estimation of Semi-Nonparametric Models.” In Handbook of Econometrics, edited by James J. Heckman and Edward E. Leamer, 6B:5549–5632. Elsevier. https://doi.org/10.1016/S1573-4412(07)06076-X.
Christensen, Kim, Mathias Siggaard, and Bezirgen Veliyev. 2023. “A machine learning approach to volatility forecasting.”Journal of Financial Econometrics 21 (5): 1680–1727. https://doi.org/10.1093/jjfinec/nbac020.
Cybenko, G. 1989. “Approximation by Superpositions of a Sigmoidal Function.”Mathematics of Control, Signals, and Systems 2 (4): 303–14. https://doi.org/10.1007/BF02551274.
Glorot, Xavier, Antoine Bordes, and Yoshua Bengio. 2011. “Deep Sparse Rectifier Neural Networks.” In Proceedings of the 14th International Conference on Artificial Intelligence and Statistics, 15:315–23. Proceedings of Machine Learning Research. PMLR.
Gneiting, Tilmann. 2011. “Making and Evaluating Point Forecasts.”Journal of the American Statistical Association 106 (494): 746–62. https://doi.org/10.1198/jasa.2011.r10138.
Gneiting, Tilmann, and Adrian E. Raftery. 2007. “Strictly proper scoring rules, prediction, and estimation.”Journal of the American Statistical Association 102 (477): 359–78. https://doi.org/10.1198/016214506000001437.
Leshno, Moshe, Vladimir Ya. Lin, Allan Pinkus, and Shimon Schocken. 1993. “Multilayer Feedforward Networks with a Nonpolynomial Activation Function Can Approximate Any Function.”Neural Networks 6 (6): 861–67. https://doi.org/10.1016/S0893-6080(05)80131-5.
Richardson, Adam, Thomas van Florenstein Mulder, and Tuğrul Vehbi. 2021. “Nowcasting GDP using machine-learning algorithms: A real-time assessment.”International Journal of Forecasting 37 (2): 941–48. https://doi.org/10.1016/j.ijforecast.2020.10.005.
Rumelhart, David E., Geoffrey E. Hinton, and Ronald J. Williams. 1986. “Learning Representations by Back-Propagating Errors.”Nature 323: 533–36. https://doi.org/10.1038/323533a0.
Yeh, I-Cheng, and Che-hui Lien. 2009. “The Comparisons of Data Mining Techniques for the Predictive Accuracy of Probability of Default of Credit Card Clients.”Expert Systems with Applications 36 (2): 2473–80. https://doi.org/10.1016/j.eswa.2007.12.020.
