Under active development. A stable version is expected in November 2026. Feedback is welcome by email or via GitHub issues.

12 Gradient Boosting

12.1 Overview

Gradient boosting is a sequential tree-based method that builds an additive predictor by repeatedly fitting new trees to the mistakes of the current ensemble. Random forests are designed primarily to stabilize high-variance trees through averaging, whereas boosting often emphasizes bias reduction through successive corrections. In practice, both methods can change both bias and variance. The modern gradient-boosting formulation is due to J. H. Friedman (2001); adaptive boosting (AdaBoost) is an important predecessor and special case in the broader boosting literature (Freund and Schapire 1997).

For econometricians, boosting is a flexible predictive tool for tabular data. It can learn nonlinearities, interactions, threshold effects, and asymmetric predictive relationships without the researcher specifying them. It is well suited to settings where the signal is spread across many variables: representing such structure with one regularized tree may require many splits and produce an unstable approximation.

At the same time, boosting has several interacting tuning choices. Learning rate, tree depth, early stopping, and validation design are all consequential. In time-series and forecasting applications, careless validation can create large but spurious performance gains. A common credit-default application near the end of the chapter compares boosting with a single tree, a random forest, logistic regression, and a feed-forward neural network under the same data split and probability scores.

12.2 Roadmap

1. We begin with the additive structure of boosting and contrast it with random forests.
2. We then derive the squared-error version, where each new tree fits residuals.
3. Next we generalize to arbitrary differentiable losses, with each new tree fitted to the negative loss gradients—the pseudo-residuals.
4. We then discuss shrinkage (scaling each tree’s contribution by a learning rate), tree depth, subsampling, and early stopping.
5. We then revisit credit-card default prediction to compare the tree-based methods with logistic regression and a feed-forward neural network.
6. Finally, we summarize how boosting should be interpreted and validated in econometric applications.

12.3 Additive Stagewise Modeling

Following the standard notation in the boosting literature, this chapter writes F_m for the ensemble prediction function after iteration m. This chapter-local notation is distinct from the use of F for a predictive cumulative distribution function in earlier chapters. In regression, F_m(x) is a point prediction; in binary classification below, it is a real-valued score.

Boosting constructs the predictor

\hat F_M(x) = \hat F_0(x) + \nu \sum_{m=1}^M \hat h_m(x),

where:

• \hat F_0 is an initial simple predictor
• \hat h_m is the tree added at iteration m
• \nu \in (0,1] is the learning rate
• M is the number of boosting iterations

Each tree here plays the role of a weak learner (equivalently, a base learner): a deliberately low-complexity base model—in this chapter a shallow regression tree. The name comes from classification boosting, where a weak learner performs only slightly better than random guessing. Gradient boosting imposes no such requirement on any tree in isolation. Each tree’s job is instead to move the current ensemble in a descent direction of the loss, so it is fitted to the current residuals or pseudo-residuals rather than directly to the original target.

Random Forests vs Boosting

The two methods use the same basic building block, the decision tree, but in different ways:

1. Random forests grow many trees largely independently and average them, primarily to stabilize high-variance tree predictions.
2. Boosting grows trees sequentially, with each new tree correcting the current ensemble; this often reduces bias, but can also change variance.
3. Random forests do not have a learning-rate–iteration-count tradeoff; boosting requires these two choices to be tuned jointly.

Method	Ensemble structure	Main regularization or stabilization lever	Main risk
Single tree	No ensemble	Tree depth, leaf size, pruning	High variance and discontinuity
Bagging	Parallel bootstrap trees	Tree size; B controls finite-ensemble simulation error	Trees may remain highly correlated
Random forest	Parallel bootstrap trees with feature subsampling	max_features and leaf size; B controls finite-ensemble simulation error	Out-of-bag (OOB) error can mislead under dependence
Gradient boosting	Sequential trees fit to residuals or pseudo-residuals	Learning rate, number of trees, tree depth, early stopping	Sensitive tuning and validation leakage

12.4 Squared-Error Boosting

Start with regression and squared loss

L(F) = \frac{1}{2}\sum_{i=1}^N (y_i - F(x_i))^2.

The factor \frac{1}{2} is algebraically convenient and does not affect the minimizer.

Two conventions in this display are worth flagging, since they differ from earlier chapters. First, L here is the total training loss, summed and not averaged over observations—the convention of the boosting literature, and the one that makes the pseudo-residuals below come out as y_i - F(x_i) without a trailing 1/N. The optimization chapter instead wrote L(\boldsymbol\theta) = N^{-1}\sum_{i=1}^N \ell_i(\boldsymbol\theta) for the sample average; the two differ by the constant factor N, which changes neither the minimizer nor the direction of the gradient. Second, the per-observation loss is written outcome first, \ell(y_i, F(x_i)), as in the feed-forward networks chapter.

Step 0: Initial Model

The best constant predictor minimizes

\frac{1}{2}\sum_{i=1}^N (y_i - c)^2.

Differentiating with respect to c gives

\sum_{i=1}^N (c - y_i) = 0 \quad \Longrightarrow \quad \hat F_0(x) \equiv \bar y.

So the initial model is simply the sample mean.

Step 1: Residual Fitting

At iteration m, define the residuals of the current ensemble as

r_{im} = y_i - \hat F_{m-1}(x_i).

Under squared loss, these are exactly the negative gradients of the objective with respect to the current fitted values, since differentiating L(F) gives -\,\partial L/\partial F(x_i) = y_i - F(x_i), evaluated at F = \hat F_{m-1}. Boosting therefore fits the next tree \hat h_m(x) to the residuals:

\hat h_m \approx \arg\min_h \sum_{i=1}^N (r_{im} - h(x_i))^2.

The update is then

\hat F_m(x) = \hat F_{m-1}(x) + \nu \hat h_m(x).

If \nu = 1, the correction is taken in full. If \nu < 1, the update is shrunken.

For an econometrician, this is shrunken forward-stagewise regression. If the candidate weak learners were single regressors rather than trees, each stage would run a low-complexity ordinary least squares regression of the current residuals on the candidate regressors, add only a fraction \nu of the selected fitted value, and then recompute the residuals. This is the componentwise L_2-boosting connection developed by Bühlmann and Yu (2003). A shallow tree changes the function class from one regressor to a piecewise-constant rule, but the stagewise residual-regression logic is the same.

Squared-Error Boosting Algorithm

For the regression-tree case, the algorithm can be summarized as follows.

1. Initialize \hat F_0(x)=\bar y and residuals r_i=y_i-\bar y.
2. For m=1,\ldots,M:
   • fit a shallow tree \hat h_m with maximum depth d to the current residuals
   • update \hat F_m(x)=\hat F_{m-1}(x)+\nu \hat h_m(x)
   • update residuals r_i=y_i-\hat F_m(x_i)
3. Return the boosted predictor \hat F_M(x).

The maximum depth d bounds the interaction order available along a root-to-leaf path, while \nu controls how cautiously each correction is added.

Figure 12.1 follows the fitted function from the initial constant through two later boosting stages.

import numpy as np
import matplotlib.pyplot as plt
from sklearn.ensemble import GradientBoostingRegressor

rng = np.random.default_rng(202)
x = np.sort(rng.uniform(-3, 3, 200))
y = np.sin(1.3 * x) + 0.35 * x + 0.8 * (x > 0.8) + rng.normal(scale=0.35, size=x.shape[0])

gbr = GradientBoostingRegressor(
    n_estimators=60,
    learning_rate=0.1,
    max_depth=2,
    min_samples_leaf=8,
    random_state=202,
)
gbr.fit(x.reshape(-1, 1), y)

x_grid = np.linspace(-3.1, 3.1, 500)
staged = list(gbr.staged_predict(x_grid.reshape(-1, 1)))
f0 = np.full_like(x_grid, np.mean(y))

fig, axes = plt.subplots(1, 3, figsize=(13, 4.5), sharey=True)

for ax, pred, title in zip(
    axes,
    [f0, staged[2], staged[39]],
    ["Initial constant model", "After 3 trees", "After 40 trees"],
):
    ax.scatter(x, y, s=15, color="0.4", alpha=0.65)
    ax.plot(x_grid, pred, linewidth=2.4, color="C3")
    ax.set_title(title)
    ax.set_xlabel("Predictor")
    ax.grid(True, alpha=0.3)

axes[0].set_ylabel("Outcome")
plt.tight_layout()
plt.show()

Figure 12.1: Training observations (gray points) and the boosted fit (red line) at three stages: the initial constant model (left), after 3 trees (center), and after 40 trees (right).

The ensemble starts as a constant. With each iteration, the fit becomes more refined because the next tree is trained on what the current ensemble still gets wrong.

Why Small Trees Are Used

Boosting usually uses shallow trees rather than deep ones. A decision stump is a depth-1 tree with one split, so it captures a single threshold effect. A tree of depth 2 or 3 can capture low-order interactions. This is deliberate regularization:

• shallow trees make each correction simple
• many simple corrections can approximate a complicated function
• restricting tree depth reduces the risk of overfitting each stage

This stagewise logic is one reason boosting performs well on tabular data when the regression function has low- to moderate-order interactions.

12.5 General Gradient Boosting

Gradient boosting extends beyond squared loss by replacing the residual with the negative gradient of a differentiable loss. Suppose the objective is

\sum_{i=1}^N \ell(y_i, F(x_i)),

where \ell is differentiable with respect to the fitted value F(x_i).

Before the stagewise updates begin, choose the loss-optimal constant predictor

\hat F_0(x)\equiv\hat\gamma_0, \qquad \hat\gamma_0 \in \arg\min_{\gamma\in\mathbb{R}} \sum_{i=1}^N \ell(y_i,\gamma).

For the losses used in this chapter, this constant is the sample mean under squared loss, an empirical \tau-quantile under quantile (pinball) loss, a Huber location estimate under Huber loss, and the sample log-odds under Bernoulli log-loss when both classes occur. The membership symbol \in allows for losses, such as pinball loss, whose sample minimizer need not be unique.

At iteration m, the pseudo-residual is

r_{im} = - \left. \frac{\partial \ell(y_i, F(x_i))}{\partial F(x_i)} \right|_{F=\hat F_{m-1}}.

The next tree is fitted to these pseudo-residuals, and the ensemble is updated in the direction of steepest descent in function space. That is why the method is called gradient boosting.

Additional Detail: Pinball and Huber Losses

The differentiability assumption can be relaxed. Two useful examples show how the loss changes both the statistical target and the pseudo-residuals.

Pinball loss for conditional quantiles. Let u_i=y_i-F(x_i) be the current forecast error and choose a quantile level \tau\in(0,1). The pinball loss is

\rho_\tau(u_i) = \begin{cases} \tau u_i, & u_i\geq 0,\\ (\tau-1)u_i, & u_i<0. \end{cases}

The loss penalizes underprediction at rate \tau and overprediction at rate 1-\tau. In population—for a random draw (X,Y) whose realizations the training pairs (x_i,y_i) represent—minimizing the conditional expected value of \rho_\tau(Y-F(x)) targets a conditional \tau-quantile. If the conditional quantile is not unique, any q satisfying \mathbb{P}(Y<q\mid X=x)\leq\tau\leq\mathbb{P}(Y\leq q\mid X=x) is a valid minimizer. When \tau=0.5, the loss is proportional to absolute error and targets a conditional median. The two linear pieces meet at u_i=0, where the ordinary derivative does not exist. Away from that kink, the negative derivative with respect to F(x_i) is

r_{im} = \begin{cases} \tau, & y_i>\hat F_{m-1}(x_i),\\ \tau-1, & y_i<\hat F_{m-1}(x_i). \end{cases}

At equality, the loss subgradient with respect to F(x_i) can be any value in [-\tau,1-\tau], so the corresponding pseudo-residual—its negative—can be any value in [\tau-1,\tau]. A subgradient is a generalized slope for a non-differentiable point. Implementations choose one valid pseudo-residual and fit the next tree to these asymmetric targets.

Huber loss for robust location prediction. For a threshold \delta>0, the Huber loss is quadratic for small errors and linear for large ones:

H_\delta(u_i) = \begin{cases} \tfrac{1}{2}u_i^2, & |u_i|\leq\delta,\\ \delta\bigl(|u_i|-\tfrac{1}{2}\delta\bigr), & |u_i|>\delta. \end{cases}

Its negative loss gradient with respect to F(x_i) is the clipped residual \max\{-\delta,\min(u_i,\delta)\}. Large forecast errors therefore cannot generate arbitrarily large pseudo-residuals. The threshold \delta is measured on the outcome scale: smaller values create more clipping and a more median-like Huber location target, while larger values approach squared loss and a more mean-like target. It should therefore be fixed or selected in a scale-aware, validation-safe way. The general lesson is that gradient boosting is not tied to squared, Bernoulli, or exponential loss: the chosen loss and its tuning constants determine the population target—the minimizer of the population risk \mathbb{E}[\ell(Y, F(X))] over the random draw (X,Y)—as well as the pseudo-residuals and the loss-optimal leaf values.

To see the gradient structure precisely, think of F as a vector (F(x_1), \ldots, F(x_N))^\top \in \mathbb{R}^N, where each coordinate is the fitted value at one training observation. The total loss is a function of this vector, L(F) = \sum_{i=1}^N \ell(y_i, F(x_i))—the same training objective introduced for squared loss in Section 12.4, now written for a general \ell—and the i-th coordinate of its negative gradient at \hat F_{m-1} is exactly r_{im}. Boosting cannot step directly in this direction, because the result would be N numbers with no way to generalize to new observations. Instead, each tree \hat{h}_m approximates the negative gradient within a structured, tree-shaped function class. Fitting a tree to the pseudo-residuals is therefore fitting a generalizable approximation to the steepest-descent direction.

Refining the leaf values. Fitting a least-squares tree to the pseudo-residuals serves to determine the split structure—the partition of the predictor space into leaves R_{jm}, j=1,\dots,J_m (the j-th of the J_m leaves of the tree fitted at iteration m)—but the least-squares leaf means need not be the loss-optimal constants on the scale of \ell. The loss-optimal leaf values are obtained by a separate one-dimensional minimization within each leaf,

\gamma_{jm}^{\star} \in \arg\min_{\gamma} \sum_{i:\, x_i \in R_{jm}} \ell\big(y_i,\, \hat F_{m-1}(x_i) + \gamma\big),

after which the update becomes \hat F_m = \hat F_{m-1} + \nu \sum_{j=1}^{J_m} \gamma_{jm}^{\star} \mathbf{1}\{x \in R_{jm}\} rather than \hat F_{m-1} + \nu \hat h_m. For squared error, \gamma_{jm}^{\star} is the mean pseudo-residual in the leaf, so the two updates coincide. For pinball loss, the minimizer can be set-valued; an implementation selects one empirical quantile according to its tie convention. The distinction between the tree’s least-squares leaf means and the loss-optimal values matters for general losses such as the Bernoulli loss introduced next. This per-leaf line search is the TreeBoost refinement proposed by J. H. Friedman (2001).

Classification with Bernoulli Log-Loss

For binary outcomes, let F(x_i)\in\mathbb{R} denote the ensemble’s unrestricted score. The sigmoid function converts this score into an event probability:

p_i = \sigma(F(x_i)) = \frac{1}{1 + e^{-F(x_i)}}.

For any probability p_i\in(0,1), the odds of the event are p_i/(1-p_i). The log-odds, also called the logit, are the logarithm of those odds. Because the logit is the inverse of the sigmoid,

F(x_i) = \log\left(\frac{p_i}{1-p_i}\right).

Thus F(x_i) is measured on the log-odds scale rather than the probability scale: it can take any real value, while the sigmoid maps it into a probability between zero and one. Adding one unit to F(x_i) multiplies the event odds by e; it does not add a fixed amount to the probability. The Bernoulli log-loss is

\ell(y_i, F(x_i)) = - y_i \log p_i - (1-y_i)\log(1-p_i).

Differentiating gives

\frac{\partial \ell}{\partial F(x_i)} = p_i - y_i.

At boosting iteration m this derivative is evaluated at the current score, so write p_{i,m-1} = \sigma(\hat F_{m-1}(x_i)) for the fitted probability that observation i carries into iteration m. The probability changes at every iteration, and the subscript keeps that visible. The negative gradient is then

r_{im} = y_i - p_{i,m-1}.

Classification boosting therefore retains a residual-like target, y_i-p_{i,m-1}, formed from the binary outcome and the current probability. The ensemble score F, however, lives on the log-odds scale. The Newton leaf calculation below uses Bernoulli curvature to translate the probability-scale discrepancy into an update of that log-odds score.

The Bernoulli case is exactly where the leaf-value refinement above matters: the tree fitted to the pseudo-residuals r_{im} = y_i - p_{i,m-1} fixes only the splits, and because Bernoulli log-loss is not squared error, the exact leaf constants \gamma_{jm}^{\star} are set by the per-leaf line search rather than by the mean residual in the leaf. That minimization has no closed form in general. Implementations commonly approximate \gamma_{jm}^{\star} with one Newton step from \gamma=0,

\widetilde\gamma_{jm} = \frac{\sum_{i:\, x_i \in R_{jm}} (y_i - p_{i,m-1})}{\sum_{i:\, x_i \in R_{jm}} p_{i,m-1}(1-p_{i,m-1})},

using the gradient p_{i,m-1} - y_i and the curvature p_{i,m-1}(1-p_{i,m-1}) evaluated at the probabilities carried in from iteration m-1. The tilde distinguishes this one-step approximation from the exact minimizer \gamma_{jm}^{\star}.

The scaling becomes important when the current probabilities are extreme. The curvature contribution p_{i,m-1}(1-p_{i,m-1}) is close to zero when p_{i,m-1} is close to zero or one. If such an extreme prediction is wrong, the corresponding residual y_i-p_{i,m-1} is instead close to 1 or -1. A leaf containing this confidently misclassified observation can therefore have a non-negligible numerator but little total curvature in the denominator, producing a large correction on the log-odds scale.

The denominator is also the observed Fisher curvature for a one-parameter Bernoulli log-odds model within the leaf. Thus the Newton leaf update is a local Fisher-scoring step: it divides the score by the amount of statistical curvature in that leaf. The natural-gradient section of Chapter 13 generalizes the same information-scaling idea to a vector of distribution parameters.

Connection to Cross-Entropy

For classification, boosting with Bernoulli log-loss is minimizing cross-entropy, i.e., the negative log-likelihood of a Bernoulli model. The equivalence links boosting directly to the maximum-likelihood and Kullback–Leibler discussion: fitting the boosted classifier means reducing the discrepancy between observed class outcomes and predicted event probabilities.

AdaBoost as Exponential-Loss Boosting

AdaBoost predates gradient boosting, but exponential loss places the two algorithms in a common framework. Freund and Schapire (1997) introduced AdaBoost. J. Friedman, Hastie, and Tibshirani (2000) later connected its observation-weighting rule to exponential-loss boosting and interpreted the resulting additive model as additive logistic regression.

For binary classification with y_i \in \{-1, +1\} for i=1,\ldots,N, define

\ell(y_i, F(x_i)) = e^{-y_i F(x_i)}.

The negative gradient at the current fit is

r_{im} = -\left.\frac{\partial \ell(y_i, F(x_i))}{\partial F(x_i)}\right|_{F=\hat F_{m-1}} = y_i \, e^{-y_i \hat F_{m-1}(x_i)}.

Weighted-classification equivalence. Define the observation weight

w_{im}=e^{-y_i\hat F_{m-1}(x_i)},

so that r_{im}=y_iw_{im}. Suppose the weak learner is a classifier with h(x_i)\in\{-1,+1\}. Because h(x_i)^2=1, its least-squares criterion satisfies

\sum_{i=1}^N\bigl(r_{im}-h(x_i)\bigr)^2 = C_m-2\sum_{i=1}^N w_{im}y_i h(x_i), \qquad C_m=\sum_{i=1}^N\bigl(w_{im}^2+1\bigr).

The constant C_m does not depend on h. Moreover,

\mathbf{1}\{y_i\neq h(x_i)\} = \frac{1-y_i h(x_i)}{2}.

Minimizing the least-squares criterion is therefore equivalent to minimizing weighted misclassification error. This equivalence relies on the restriction h(x_i)\in\{-1,+1\}; it does not apply to the real-valued regression-tree leaves used elsewhere in this chapter.

The signed margin of observation i is y_i\hat F_{m-1}(x_i). A negative margin indicates misclassification and implies w_{im}>1, so the next classifier places more weight on the current mistakes. This is the AdaBoost observation-reweighting step. The full AdaBoost algorithm also chooses

\alpha_m = \frac{1}{2}\log\!\left(\frac{1-\mathrm{err}_m}{\mathrm{err}_m}\right),

where \mathrm{err}_m is the weak classifier’s weighted misclassification rate. The correspondence above concerns the observation weights, not this step-size choice.

What the Bernoulli residual bound controls. Exponential loss assigns an unbounded weight w_{im} to an observation whose negative margin keeps growing. That observation can eventually dominate the first-order fitting criterion. Under Bernoulli log-loss, by contrast, the pseudo-residual satisfies |y_i-p_{i,m-1}|<1. The bound prevents this exponentially increasing first-order contribution at any one iteration.

The bound does not guarantee that the final ensemble is insensitive to one observation. Split selection is discrete: one observation can change which of two similar candidate splits is selected and thereby affect every case routed through that node. The observation also re-enters the fitting problem at each iteration, so its bounded contributions can influence several successive trees.

The pseudo-residual bound also does not bound the approximate Newton leaf value \widetilde\gamma_{jm}. In the Newton update derived above, the denominator is the leaf’s total curvature,

\sum_{i:\,x_i\in R_{jm}}p_{i,m-1}(1-p_{i,m-1}).

The fitting procedure does not guarantee that this sum is bounded away from zero. For a confidently misclassified observation, the numerator contribution approaches 1 or -1, while its curvature contribution approaches zero. The numerator can therefore remain non-negligible when the denominator is small, producing a large leaf update.

Different safeguards address different channels of influence. Minimum leaf sizes discourage splits based on very few observations; a lower bound on total leaf curvature rules out nearly flat leaves; an update cap bounds |\widetilde\gamma_{jm}| directly; and penalties on leaf values shrink the update. Log-loss prevents the unbounded per-observation first-order weight generated by exponential loss, but it does not make the fitted ensemble insensitive to individual observations.

12.6 Regularization in Boosting

Boosting can keep refining the fit indefinitely. That is why regularization is essential, not optional.

Learning Rate

The learning rate \nu shrinks each update:

\hat F_m(x) = \hat F_{m-1}(x) + \nu \hat h_m(x).

Smaller \nu means:

• slower learning
• more trees needed
• often better generalization, because the fit evolves more cautiously

Scaling each update by \nu<1 is called shrinkage. The interaction between \nu and the optimal number of trees M is roughly reciprocal: a smaller \nu requires more iterations to reach the same training risk, so halving \nu approximately doubles the optimal M chosen by validation. This reciprocal relation is an empirical regularity reported by J. H. Friedman (2001), not a theorem. Zhang and Yu (2005) supply theoretical support for treating the iteration count as a regularization device: for boosting over the linear span of a base class, with independent and identically distributed data and particular early-stopping rules, they establish numerical convergence and consistency under regularity conditions. This is a result about a specific class of stopping strategies, not a blanket guarantee that any validation-based stopping rule applied to boosted trees on dependent data inherits the same property. The product \nu \cdot M behaves like a complexity budget, which is the sense in which shrinkage acts as a regularization device.

Number of Trees

The number of iterations M controls how long the algorithm keeps refining the fit.

• too few trees can underfit
• too many trees can eventually overfit
• the optimal M is usually chosen by validation or early stopping—that is, by stopping boosting at the iteration where validation error is minimized (see Section 12.7)

Tree Depth

Depth controls the maximum interaction order available along a root-to-leaf path.

• depth 1: additive threshold effects only
• depth 2: two-way interactions can appear, but need not
• deeper trees: higher-order local interactions are permitted, but repeated variables or additive structure can keep the realized interaction order lower

Subsampling

Some boosting implementations fit each new tree on a random subsample of the data rather than the full sample. This is called stochastic gradient boosting (J. H. Friedman 2002).

Subsampling can reduce variance and speed computation by making each stage cheaper and less tied to one full-sample pseudo-residual fit. It introduces additional randomness into the stagewise updates, so any predictive benefit must be checked over repeated seeds and with the appropriate validation design.

Monotonicity Constraints

Sometimes the econometrician has credible prior knowledge about the sign of a predictive relationship even if the rest of the conditional mean is complicated. For example, one might regard it as implausible that a higher debt-service burden should lower predicted default risk, or that a tighter financial-conditions index should lower predicted recession risk, all else equal. A monotonicity constraint imposes exactly such a shape restriction on the fitted function.

In a tree-ensemble model, this means restricting the fitted prediction to move only in one direction as a chosen predictor increases, holding the other predictors fixed. The idea is not unique to boosting, but it is especially prominent in modern boosted-tree implementations because it fits naturally into the broader regularization toolkit. This does not make the model structural or causal. It is better understood as a regularization device that uses economic theory to rule out locally implausible wiggles in the prediction surface.

When the sign restriction is well founded, monotonicity constraints can improve credibility, stabilize partial-effect summaries, and prevent the ensemble from fitting noise in directions that contradict economic reasoning. But they are not free. If the true conditional relationship is non-monotone, changes sign across regimes, or is only monotone after conditioning on variables the model does not include, the constraint can introduce systematic misspecification.

How the constraint is enforced. A monotonicity constraint is a hard restriction on the fitted function, not a penalty added to the loss. Implementations restrict candidate splits and leaf values so that the resulting ensemble respects the requested ordering, but the precise bookkeeping is library-specific. The scikit-learn documentation linked below states the supported constraint and illustrates the resulting global shape restriction; other libraries may use different tree-building algorithms to obtain the same functional property.

In scikit-learn. The histogram-based gradient-boosting estimator HistGradientBoostingRegressor accepts a monotonic_cst argument with one integer per feature: +1 for increasing, -1 for decreasing, and 0 for no constraint. For three predictors, with the first constrained to have a monotonically increasing effect:

from sklearn.ensemble import HistGradientBoostingRegressor

model = HistGradientBoostingRegressor(
    monotonic_cst=[1, 0, 0],
    max_iter=500,
    learning_rate=0.05,
)

See the scikit-learn example Monotonic Constraints for a visual comparison of constrained and unconstrained fits and a full working interface example.

Implementation Note: Extreme Gradient Boosting (XGBoost)

Extreme gradient boosting (XGBoost) (Chen and Guestrin 2016) is an optimized implementation family for gradient-boosted trees. Relative to the basic algorithm in this chapter, it adds engineering and regularization features such as efficient split search, \ell_1 and \ell_2 penalties, sparse or missing-value handling, monotone constraints, and scalable training. The econometric interpretation remains the same: it is still a tuned predictive tree ensemble, so validation design and the information set are central.

Figure 12.2 plots training and validation mean squared error (MSE) along the boosting path for two learning rates. Because the held-out sample is inspected at every iteration, it is a tuning sample rather than a final test sample.

import numpy as np
import matplotlib.pyplot as plt
from sklearn.ensemble import GradientBoostingRegressor

rng = np.random.default_rng(77)
n_train = 280
n_validation = 280

X_train = rng.uniform(-3, 3, size=(n_train, 1))
X_validation = rng.uniform(-3, 3, size=(n_validation, 1))

def g(x):
    z = x[:, 0]
    return np.sin(1.5 * z) + 0.4 * z + 0.7 * (z > 1.0)

y_train = g(X_train) + rng.normal(scale=0.35, size=n_train)
y_validation = g(X_validation) + rng.normal(scale=0.35, size=n_validation)

configs = [(0.05, "C0"), (0.2, "C3")]

fig, axes = plt.subplots(1, 2, figsize=(12, 4.8), sharey=True)

for lr, color in configs:
    model = GradientBoostingRegressor(
        n_estimators=180,
        learning_rate=lr,
        max_depth=2,
        min_samples_leaf=6,
        random_state=77,
    )
    model.fit(X_train, y_train)
    initial_prediction = np.mean(y_train)
    train_mse = [np.mean((y_train - initial_prediction) ** 2)]
    validation_mse = [np.mean((y_validation - initial_prediction) ** 2)]
    for pred_train, pred_validation in zip(
        model.staged_predict(X_train), model.staged_predict(X_validation)
    ):
        train_mse.append(np.mean((y_train - pred_train) ** 2))
        validation_mse.append(np.mean((y_validation - pred_validation) ** 2))

    tree_counts = np.arange(len(train_mse))
    ax = axes[0] if lr == 0.05 else axes[1]
    ax.plot(tree_counts, train_mse, label="Training MSE", color="C0", linewidth=2)
    ax.plot(tree_counts, validation_mse, label="Validation MSE", color="C3", linewidth=2)
    ax.set_title(f"Learning rate = {lr}")
    ax.set_xlabel("Number of trees")
    ax.grid(True, alpha=0.3)
    ax.legend(frameon=False)

axes[0].set_ylabel("Mean squared error")
plt.tight_layout()
plt.show()

Figure 12.2: Training MSE (blue) and independent validation-sample MSE (red) from the initial constant model at zero trees through 180 boosting iterations. The validation path is available for choosing the iteration count, so it is not a final test evaluation. The left panel uses a learning rate of 0.05; the right panel uses 0.2.

The smaller learning rate often gives a flatter, more forgiving validation-error profile. With the larger learning rate, validation error initially falls more quickly but can begin rising earlier.

Question for Reflection

If two boosted models have similar validation error, what evidence from the validation curve would justify the small-learning-rate, many-tree model, and what evidence would justify the larger-learning-rate, fewer-tree model?

Suggested Answer

A small learning rate with many trees is better supported when its validation curve is flatter and overfitting appears later along the boosting path. A larger learning rate with fewer trees is better supported when the two models have genuinely similar time-ordered validation error, but the larger-learning-rate model reaches that performance with fewer boosting steps and lower computational cost. The comparison should be based on validation over the relevant forecast origins, not on in-sample fit.

Real macroeconomic data can be less forgiving. Figure 12.3 shows training and validation MSE curves for annualized gross domestic product (GDP) growth in quarter t+1 from the quarterly Federal Reserve Economic Data collection (FRED-QD) developed by McCracken and Ng (2020); the data appendix documents the repository copy. The five predictors are dated quarter t: GDP growth, inflation, the term spread, the log volatility index (VIX), and the unemployment rate. The evaluation uses a fixed-estimation-window, time-ordered 75/25 holdout: each model is fitted once on the first 75% of quarters and evaluated, without re-estimation, on the remaining 25%.¹ Training MSE falls steadily in both panels. Neither validation curve improves on the stage-zero constant forecast: both attain their global minimum at zero trees, although they fluctuate locally thereafter. The learning rate governs how fast the deterioration accumulates—gradually at \nu = 0.02, faster and more erratically at \nu = 0.2.

The immediate deterioration is worth confronting rather than hiding. The lesson of early stopping is not that boosting always improves for a while and then degrades; it is that the iteration count is a regularization dial whose right setting is chosen on held-out data, and here the honest setting is no boosting at all. That zero-tree conclusion is conditional on this split and on keeping the pre-validation model frozen; a recursive or rolling design that refits at each forecast origin estimates a different deployment risk and could select a different stopping point. The validation block is the last 25% of quarters, so it contains the pandemic period, while the training data contain no comparable episode. In this split, additional flexibility fitted to pre-2020 fluctuations did not help and moved forecasts away from the training-sample mean, which happened to be the safer prediction under this break. The absolute MSE values are correspondingly large—the 2020 contraction alone generates residuals of tens of annualized percentage points that dominate the average. The gap between the training and validation levels is consistent with a pronounced distribution shift, although one realized split cannot by itself establish a formal regime change. The comparison with the synthetic example above is instructive precisely because it differs: there, the data-generating process is stable across the split and the validation curve shows the textbook U-shape.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.ensemble import GradientBoostingRegressor

raw = pd.read_csv("data/fred_qd_current.csv")
fred_qd = raw.iloc[2:].copy()
fred_qd["sasdate"] = pd.to_datetime(fred_qd["sasdate"])
for col in fred_qd.columns:
    if col != "sasdate":
        fred_qd[col] = pd.to_numeric(fred_qd[col], errors="coerce")
fred_qd = fred_qd.sort_values("sasdate").reset_index(drop=True)

gdp_growth  = 400 * np.log(fred_qd["GDPC1"]).diff()
inflation   = 400 * np.log(fred_qd["CPIAUCSL"]).diff()
term_spread = fred_qd["GS10"] - fred_qd["TB3MS"]
log_vix     = np.log(fred_qd["VIXCLSx"])

macro = pd.DataFrame({
    "y":           gdp_growth.shift(-1),
    "gdp_lag1":    gdp_growth,
    "infl_lag1":   inflation,
    "spread_lag1": term_spread,
    "vix_lag1":    log_vix,
    "unrate_lag1": fred_qd["UNRATE"],
}).dropna()

feature_cols = ["gdp_lag1", "infl_lag1", "spread_lag1", "vix_lag1", "unrate_lag1"]
split = int(len(macro) * 0.75)
X_train = macro.iloc[:split][feature_cols].values
y_train = macro.iloc[:split]["y"].values
X_validation = macro.iloc[split:][feature_cols].values
y_validation = macro.iloc[split:]["y"].values

configs = [(0.02, "C0"), (0.2, "C3")]
fig, axes = plt.subplots(1, 2, figsize=(12, 4.5), sharey=True)

for (lr, _), ax in zip(configs, axes):
    model = GradientBoostingRegressor(
        n_estimators=250, learning_rate=lr, max_depth=2,
        min_samples_leaf=6, random_state=42,
    )
    model.fit(X_train, y_train)
    initial_prediction = np.mean(y_train)
    train_mse = [np.mean((y_train - initial_prediction) ** 2)]
    val_mse = [np.mean((y_validation - initial_prediction) ** 2)]
    for p_tr, p_val in zip(model.staged_predict(X_train), model.staged_predict(X_validation)):
        train_mse.append(np.mean((y_train - p_tr) ** 2))
        val_mse.append(np.mean((y_validation - p_val) ** 2))
    tree_counts = np.arange(len(train_mse))
    ax.plot(tree_counts, train_mse, color="C0", linewidth=2, label="Training MSE")
    ax.plot(tree_counts, val_mse,   color="C3", linewidth=2, label="Validation MSE")
    ax.set_title(f"Learning rate = {lr}")
    ax.set_xlabel("Number of trees")
    ax.grid(True, alpha=0.3)
    ax.legend(frameon=False)

axes[0].set_ylabel("Mean squared error")
plt.tight_layout()
plt.show()

Figure 12.3: Training MSE (blue) and time-ordered validation MSE (red) for one-step-ahead GDP growth from FRED-QD. Both panels include the initial constant model at zero trees followed by 250 depth-2 trees; the left panel uses a learning rate of 0.02 and the right panel 0.2.

The zero-tree minima in Figure 12.3 show why an early-stopping rule must be allowed to retain the initial constant model rather than being forced to select at least one tree.

12.7 Early Stopping and Validation

Because boosting can continue to improve the training fit for many iterations, some external criterion is needed to decide when to stop. Early stopping uses a validation set or validation path to choose the number of trees, and it is one of the central regularization devices of boosting—not a minor implementation detail.

Econometric Warning

For forecasting applications, early stopping must be based on time-ordered validation data. Random K-fold validation contaminates the forecast-origin information set by allowing later observations into training folds. It can therefore select a materially different stopping point from a real-time design; when future observations reveal useful related structure, it often favors an overly flexible path, but the direction is data-dependent under drift or structural change.

The validation design should mirror the actual forecast exercise:

• rolling or expanding windows for time series
• group-aware validation for clustered or panel-style dependence
• leakage-safe preprocessing inside each training fold

12.8 Empirical Comparison: Credit-Default Probabilities

The three tree-based chapters introduce a single tree, a random forest, and gradient boosting separately. We can now compare them on one common task. The Taiwan credit-card default data contain 30,000 clients and a binary indicator for default in the month after six months of observed repayment histories (Yeh 2009; Yeh and Lien 2009). The predictors also include the granted credit limit, bill and payment amounts, age, and coded demographic characteristics.

The design exactly reuses the stratified 60% training, 20% validation, and 20% test split from the feed-forward networks chapter. This reuse makes the model scores directly comparable, but it does not create a second independent test of the neural network. The network architecture remains fixed at the eight-unit, one-hidden-layer candidate selected there; we do not alter it after seeing its earlier test result.

The single-tree, random-forest, and boosting settings are chosen by validation log loss from small candidate sets declared in the code. The comparison uses the same transformed predictors for every model: continuous variables are standardized and categorical variables are one-hot encoded. Standardization is needed by logistic regression and the network. It does not change a tree’s available split order because it is a strictly increasing transformation, but using one representation removes an incidental difference between the model pipelines. All transformations are estimated without the test observations.

We evaluate the final probability forecasts by log loss, by the Brier score N_{\mathrm{test}}^{-1}\sum_{i=1}^{N_{\mathrm{test}}}(y_i-\hat p_i)^2, and by accuracy after applying a 0.5 threshold. The first two columns evaluate the probabilities themselves; the last evaluates one particular binary decision rule.

import numpy as np
import pandas as pd
from book_tables import display_styled_table
from sklearn.base import clone
from sklearn.compose import ColumnTransformer
from sklearn.ensemble import GradientBoostingClassifier, RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import accuracy_score, brier_score_loss, log_loss
from sklearn.model_selection import train_test_split
from sklearn.neural_network import MLPClassifier
from sklearn.preprocessing import OneHotEncoder, StandardScaler
from sklearn.tree import DecisionTreeClassifier

credit_default = pd.read_csv("data/credit_card_default_taiwan.csv")
y = credit_default["default_next_month"]
X = credit_default.drop(columns=["id", "default_next_month"])

categorical_columns = [
    "sex_code",
    "education_code",
    "marriage_code",
    "repay_status_sep",
    "repay_status_aug",
    "repay_status_jul",
    "repay_status_jun",
    "repay_status_may",
    "repay_status_apr",
]
continuous_columns = [
    column for column in X.columns if column not in categorical_columns
]

X_development, X_test, y_development, y_test = train_test_split(
    X, y, test_size=0.20, stratify=y, random_state=42,
)
X_train, X_validation, y_train, y_validation = train_test_split(
    X_development,
    y_development,
    test_size=0.25,
    stratify=y_development,
    random_state=42,
)

def make_preprocessor():
    return ColumnTransformer(
        [
            ("continuous", StandardScaler(), continuous_columns),
            (
                "categorical",
                OneHotEncoder(
                    drop="first",
                    handle_unknown="ignore",
                    sparse_output=False,
                ),
                categorical_columns,
            ),
        ]
    )

# Use training-only transformations for model selection.
selection_preprocessor = make_preprocessor()
X_train_transformed = selection_preprocessor.fit_transform(X_train)
X_validation_transformed = selection_preprocessor.transform(X_validation)

def select_by_log_loss(candidates):
    fitted_candidates = []
    for setting, model in candidates:
        model.fit(X_train_transformed, y_train)
        probability = model.predict_proba(X_validation_transformed)[:, 1]
        fitted_candidates.append(
            (log_loss(y_validation, probability), setting, model)
        )
    return min(fitted_candidates, key=lambda result: result[0])

tree_candidates = [
    (
        f"depth {depth}, leaf {leaf}",
        DecisionTreeClassifier(
            max_depth=depth,
            min_samples_leaf=leaf,
            random_state=42,
        ),
    )
    for depth in [2, 4, 6, 8]
    for leaf in [20, 100]
]

forest_candidates = [
    (
        f"300 trees, max features {max_features}, leaf {leaf}",
        RandomForestClassifier(
            n_estimators=300,
            max_features=max_features,
            min_samples_leaf=leaf,
            n_jobs=-1,
            random_state=42,
        ),
    )
    for max_features in ["sqrt", 0.5]
    for leaf in [5, 20, 50]
]

boosting_candidates = [
    (
        f"rate {rate}, {trees} trees, depth {depth}",
        GradientBoostingClassifier(
            n_estimators=trees,
            learning_rate=rate,
            max_depth=depth,
            min_samples_leaf=20,
            random_state=42,
        ),
    )
    for rate, trees in [(0.1, 50), (0.05, 100), (0.025, 200)]
    for depth in [1, 2]
]

_, tree_setting, selected_tree = select_by_log_loss(tree_candidates)
_, forest_setting, selected_forest = select_by_log_loss(forest_candidates)
_, boosting_setting, selected_boosting = select_by_log_loss(
    boosting_candidates
)

# Re-estimate all preprocessing on the development sample before testing.
final_preprocessor = make_preprocessor()
X_development_transformed = final_preprocessor.fit_transform(X_development)
X_test_transformed = final_preprocessor.transform(X_test)

final_models = {
    "Logistic regression": (
        "No penalty",
        LogisticRegression(penalty=None, max_iter=1000),
    ),
    "FNN": (
        "Hidden layer 8",
        MLPClassifier(
            hidden_layer_sizes=(8,),
            activation="relu",
            solver="adam",
            alpha=0.0,
            batch_size=256,
            learning_rate_init=0.001,
            max_iter=500,
            tol=1e-4,
            n_iter_no_change=10,
            random_state=42,
        ),
    ),
    "Classification tree": (tree_setting, clone(selected_tree)),
    "Random forest": (forest_setting, clone(selected_forest)),
    "Gradient boosting": (boosting_setting, clone(selected_boosting)),
}

test_probabilities = {
    "Constant probability": np.full(len(y_test), y_development.mean())
}
selected_settings = {"Constant probability": "Development mean"}

for name, (setting, model) in final_models.items():
    model.fit(X_development_transformed, y_development)
    test_probabilities[name] = model.predict_proba(X_test_transformed)[:, 1]
    selected_settings[name] = setting

comparison = pd.DataFrame(
    [
        {
            "Model": name,
            "Selected setting": selected_settings[name],
            "Log loss": log_loss(y_test, probability),
            "Brier score": brier_score_loss(y_test, probability),
            "Accuracy at 0.5": accuracy_score(y_test, probability >= 0.5),
        }
        for name, probability in test_probabilities.items()
    ]
).set_index("Model")

best_cell_style = "background-color: #d9ead3; font-weight: bold;"

display_styled_table(
    comparison.style
    .format(
        {
            "Log loss": "{:.4f}",
            "Brier score": "{:.4f}",
            "Accuracy at 0.5": "{:.4f}",
        }
    )
    .highlight_min(
        subset=["Log loss", "Brier score"],
        props=best_cell_style,
    )
    .highlight_max(
        subset=["Accuracy at 0.5"],
        props=best_cell_style,
    )
)

Table 12.1: Probability-forecast performance on the common 6,000-client test sample from the Taiwan credit-card default data. Hyperparameters for the FNN and three tree-based methods are chosen by validation log loss; logistic regression has no tuned hyperparameter. All nonconstant models are then refitted on the combined training and validation samples. Lower log loss and Brier score are better; higher accuracy under the 0.5 threshold is better.

	Selected setting	Log loss	Brier score	Accuracy at 0.5
Model
Constant probability	Development mean	0.5284	0.1723	0.7788
Logistic regression	No penalty	0.4431	0.1390	0.8170
FNN	Hidden layer 8	0.4369	0.1376	0.8170
Classification tree	depth 8, leaf 100	0.4473	0.1407	0.8145
Random forest	300 trees, max features 0.5, leaf 20	0.4324	0.1356	0.8180
Gradient boosting	rate 0.1, 50 trees, depth 2	0.4386	0.1377	0.8152

The random forest has the lowest test log loss, 0.4324, and Brier score, 0.1356, in this comparison. The selected FNN is next on both probability scores, followed closely by gradient boosting. The single tree improves substantially on the constant-probability benchmark but has the weakest probability scores among the fitted models. Accuracy gives a less informative and slightly different ranking: the forest attains 81.80%, the FNN and logistic regression both attain 81.70%, and boosting attains 81.52%.

These results support a deliberately local conclusion. On this split and within these modest candidate sets, averaging randomized trees works better than relying on one partition, and it slightly outperforms the selected FNN. The table does not show that random forests universally dominate networks or boosting. The score gaps are small, the test observations are shared across models, and the data cover only one borrower cohort. A serious performance claim would require paired uncertainty assessment across fresh test observations or repeated deployment-relevant samples; a later-period cohort would also be needed to assess stability across credit regimes.

Question for Reflection

Why does the table justify saying that the random forest has the best realized log loss on this test sample, but not that it has the lowest population risk?

Suggested Answer

The first claim is a direct comparison of scores computed for the same 6,000 clients. The second concerns expected performance over new samples. One realized test set contains sampling variation, and the small score gaps could change for another cohort or credit regime. Establishing a population ranking requires uncertainty assessment or additional deployment-relevant test samples.

12.9 Strengths and Limits for Econometric Work

Boosting is a flexible off-the-shelf method for tabular prediction. It is useful when:

• nonlinearities are important
• interactions are present but unknown
• there are many candidate predictors
• forecast accuracy matters more than structural interpretability

Strengths

• stagewise fitting can combine many shallow trees into a flexible predictor
• automatic detection of nonlinearities and interactions
• flexible loss functions for means, probabilities, conditional quantiles via pinball loss, and robust location targets via Huber loss
• regularization knobs that can be tuned to the signal-to-noise environment

Limitations

• sensitive to hyperparameter tuning and validation design
• no natural extrapolation beyond the support of the training data
• slower and less transparent than a single tree
• variable importance and partial-effect summaries remain predictive, not causal

Like random forests, boosted trees should be understood as flexible predictive approximators, not as identification strategies.

12.10 Summary

Key Takeaways

1. Gradient boosting builds an additive ensemble by fitting successive trees to negative loss gradients.
2. The pseudo-residuals are ordinary residuals under squared loss and y-p under Bernoulli log-loss.
3. Learning rate, tree count, depth, subsampling and early stopping control the ensemble complexity.

Common Pitfalls

• Rapid improvement in training loss does not justify a large learning rate or additional boosting stages.
• Compare models on the same held-out observations and with the same evaluation criterion.
• Boosted-tree importance measures do not identify structural or causal effects.
• Adding more trees does not give the ensemble an ability to extrapolate smooth trends outside the training support.

12.11 Exercises

Exercise 12.1: First Step of Gradient Boosting with Shrinkage

Suppose these are four monthly forecast-origin/target pairs. At forecast origin t, the leading indicator x_t is observed and is used to predict next-month output growth y_{t+1}. You observe:

(x_t, y_{t+1}) = (1,2), (2,3), (3,2.5), (4,5).

Consider squared-error boosting with decision stumps as weak learners. Let the learning rate be \nu = 0.4.

1. Briefly explain the idea of gradient boosting under squared loss.
2. Compute the initial model F_0(x) and the residuals r_{t+1} = y_{t+1} - F_0(x_t).
3. Fit the optimal first stump h_1(x) by checking the candidate split points 1.5, 2.5, and 3.5.
4. Write the updated model F_1(x) = F_0(x) + 0.4\, h_1(x). Compute the fitted values at the four training points.
5. Compute the training sum of squared errors (SSE) of F_1. Compare it with the training SSE that would result from taking the full step \nu=1. Why can a smaller learning rate still be useful in practice?
6. Suppose the four pairs are part of a longer monthly sample ordered by forecast origin t, and the number of boosting iterations M is selected from a large grid. Explain why random K-fold validation does not estimate the intended real-time risk and can favor an overly large M. Describe a validation design that targets the intended one-step-ahead forecasting problem.
7. Explain why the direction of the random-fold distortion is not guaranteed under structural change. Use a break halfway through the sample to contrast a pre-break validation month with an early post-break validation month.

Exam level. Parts 1–5 combine boosting mechanics, shrinkage, and the regularization interpretation of the learning rate; Parts 6–7 diagnose how validation design and structural change interact with the iteration count.

Hint for Part 2

Under squared loss, the best constant predictor is the sample mean.

Hint for Part 3

For each candidate split, compute the mean residual in the left and right leaves and then the residual SSE.

Hint for Part 5

A smaller learning rate usually gives a worse fit after one step, but that is not the right comparison. Boosting is a multi-step procedure.

Hint for Parts 6–7

Ask whether every training fold would have been available before the observations in its validation fold at a real forecast origin. With a structural break halfway through the sample, compare a pre-break validation month whose training fold includes post-break observations with an early post-break validation month whose training fold includes later observations from the new regime.

Solution

Part 1: Core idea

Under squared loss, boosting starts from a simple predictor and then repeatedly fits a weak learner to the current residuals. Each new tree corrects part of the remaining error.

Part 2: Initial model and residuals

The initial model is the sample mean:

F_0(x)=\bar y=\frac{2+3+2.5+5}{4}=3.125.

Therefore the residuals are

r = (-1.125,\,-0.125,\,-0.625,\,1.875).

Part 3: Best first stump

Check the three candidate cutoffs.

For split 1.5:

• left mean residual: -1.125
• right mean residual: \frac{-0.125-0.625+1.875}{3}=0.375

Residual SSE:

0 + (-0.125-0.375)^2 + (-0.625-0.375)^2 + (1.875-0.375)^2 = 3.5.

For split 2.5:

• left mean residual: \frac{-1.125-0.125}{2}=-0.625
• right mean residual: \frac{-0.625+1.875}{2}=0.625

Residual SSE:

(-1.125+0.625)^2 + (-0.125+0.625)^2 + (-0.625-0.625)^2 + (1.875-0.625)^2 = 3.625.

For split 3.5:

• left mean residual: \frac{-1.125-0.125-0.625}{3}=-0.625
• right mean residual: 1.875

Residual SSE:

(-1.125+0.625)^2 + (-0.125+0.625)^2 + (-0.625+0.625)^2 + (1.875-1.875)^2 = 0.5.

So the best first stump is

h_1(x)= \begin{cases} -0.625, & x \leq 3.5, \\ 1.875, & x > 3.5. \end{cases}

Part 4: Shrinkage update

With \nu=0.4,

F_1(x)=F_0(x)+0.4\,h_1(x).

Therefore

F_1(x)= \begin{cases} 3.125 + 0.4(-0.625)=2.875, & x \leq 3.5, \\ 3.125 + 0.4(1.875)=3.875, & x > 3.5. \end{cases}

So the fitted values at the four training points are

(2.875,\ 2.875,\ 2.875,\ 3.875).

Part 5: Training SSE and interpretation

The training SSE is

(2-2.875)^2 + (3-2.875)^2 + (2.5-2.875)^2 + (5-3.875)^2 = 2.1875.

If instead \nu=1, then

F_1(x)= \begin{cases} 2.5, & x \leq 3.5, \\ 5, & x > 3.5, \end{cases}

and the resulting SSE is

(2-2.5)^2 + (3-2.5)^2 + (2.5-2.5)^2 + (5-5)^2 = 0.5.

So after one step, the full update fits the training sample better. But a smaller learning rate can still be useful because it regularizes the path of the algorithm. With many steps, cautious updates often generalize better out of sample than aggressive early corrections.

Part 6: Selecting the iteration count for a time series

Random K-fold validation mixes earlier and later months. A training fold can therefore contain observations dated after a validation observation. If the relationship drifts or the predictors and outcomes are serially dependent, this comparison does not reproduce the information set of a one-step-ahead forecaster and need not rank stopping points by their real-time risk. Because larger M gives the ensemble more flexibility to exploit patterns revealed by later observations, the contaminated design can make a long boosting path look too attractive and select an overly large M.

Use validation based on rolling or expanding windows instead. At every validation origin, fit the full preprocessing and boosting procedure only on observations available before that origin, evaluate the one-step-ahead loss, and choose M from the average loss across those ordered origins. This treats the iteration count as a hyperparameter and estimates it under the same timing protocol as deployment.

Part 7: Structural change and the direction of distortion

The direction is not guaranteed. Suppose a structural break occurs halfway through the sample. For a pre-break validation month, its random-fold training set mixes in post-break observations; a flexible large-M ensemble can fit the new-regime relationship and thereby predict the old-regime validation month poorly, favoring a smaller M. For an early post-break validation month, the training set can instead include later observations from the new regime that a real-time forecaster did not yet have; those observations can help a flexible ensemble and favor a larger M. In neither case does the random fold reproduce the relevant real-time training set.

Exercise 12.2: Bernoulli Log-Loss and the Initial Boosting Model

Consider binary outcomes y_i \in \{0,1\} and a boosting model with score F(x_i) and probability

p_i = \sigma(F(x_i)) = \frac{1}{1+e^{-F(x_i)}}.

The Bernoulli log-loss for observation i is

\ell(y_i,F(x_i)) = - y_i \log p_i - (1-y_i)\log(1-p_i).

1. Show that \frac{\partial \ell(y_i,F(x_i))}{\partial F(x_i)} = p_i - y_i. Conclude that the pseudo-residual is r_i = y_i - p_i.
2. Suppose the initial model is a constant, F_0(x)\equiv c, and assume the sample contains both classes, so that 0<\bar y<1. Show that the value minimizing the total log-loss satisfies \sigma(c)=\bar y, and hence c = \log\left(\frac{\bar y}{1-\bar y}\right).
3. Explain what happens to the minimization in Part 2 when \bar y = 0 or \bar y = 1. Why must software either reject a one-class training sample or impose a finite convention, for example by clipping probabilities?
4. For the sample y=(1,1,0,1), compute \bar y, the optimal initial constant F_0, the implied initial probability p_i, and the pseudo-residuals.
5. A Bernoulli leaf contains two observations with outcomes (1,0). Compute the one-step Newton leaf update when their current probabilities are (0.6,0.6) and again when they are (0.9,0.9). Use \widetilde\gamma =\frac{\sum_i(y_i-p_i)}{\sum_i p_i(1-p_i)}, where both sums run over the observations in the leaf. Explain why the second update has larger magnitude even though every pseudo-residual remains between -1 and 1.

Exam level: suitable as-is. Parts 1–2 derive the classification gradients, Part 3 checks the boundary, Part 4 applies the initialization, and Part 5 distinguishes bounded pseudo-residuals from a potentially large Newton leaf correction.

Hint for Part 1

Differentiate through the logistic link: \sigma'(z)=\sigma(z)(1-\sigma(z)).

Hint for Part 2

Differentiate the total loss with respect to the constant c and set the derivative equal to zero.

Solution

Part 1: Derivative of log-loss

Write p_i=\sigma(F_i) with F_i=F(x_i). Then

\ell_i = -y_i \log p_i - (1-y_i)\log(1-p_i).

Differentiate with respect to F_i:

\frac{\partial \ell_i}{\partial F_i} = -y_i \frac{1}{p_i}\frac{\partial p_i}{\partial F_i} - (1-y_i)\frac{1}{1-p_i}\left(-\frac{\partial p_i}{\partial F_i}\right).

Since \frac{\partial p_i}{\partial F_i}=p_i(1-p_i), this becomes

\frac{\partial \ell_i}{\partial F_i} = -y_i(1-p_i) + (1-y_i)p_i = p_i - y_i.

Therefore the negative gradient is

r_i = -\frac{\partial \ell_i}{\partial F_i} = y_i - p_i.

Part 2: Optimal constant initial model

If F_0(x)\equiv c, then every observation has the same probability p=\sigma(c). The total loss is

L(c)= -\sum_{i=1}^N \left[y_i \log p + (1-y_i)\log(1-p)\right].

From Part 1, the derivative with respect to c is

\frac{dL}{dc} = \sum_{i=1}^N (p - y_i) = Np - \sum_{i=1}^N y_i.

Setting this equal to zero gives

p = \frac{1}{N}\sum_{i=1}^N y_i = \bar y.

Since p=\sigma(c),

\sigma(c)=\bar y.

Applying the logit transformation, which is well defined because 0<\bar y<1 was assumed,

c = \log\left(\frac{\bar y}{1-\bar y}\right).

Part 3: The boundary cases

If \bar y = 1, every y_i equals one and the total loss is L(c) = -N\log\sigma(c), which is strictly decreasing in c and approaches its infimum of 0 only as c \to \infty. The first-order condition \sigma(c)=\bar y=1 has no solution, because the logistic function never attains 1: the loss has no finite minimizer, and \log(\bar y/(1-\bar y)) diverges. The case \bar y = 0 is symmetric, with c \to -\infty.

The statistical content is that an all-one or all-zero sample places the intercept-only Bernoulli likelihood at the boundary and cannot pin down a finite score. This non-existence is analogous to what happens in a logit model under complete separation, where covariates perfectly separate the two classes. Software must therefore choose how to handle a problem whose unpenalized objective has no finite minimizer. For example, implementations may clip the initial probability into [\varepsilon,1-\varepsilon].

Part 4: Numerical example

For y=(1,1,0,1),

\bar y = \frac{3}{4} = 0.75.

Hence

F_0 = \log\left(\frac{0.75}{0.25}\right)=\log 3 \approx 1.099.

The implied initial probability is

p_i = \sigma(F_0)=0.75 \qquad \text{for all } i.

So the pseudo-residuals are

r = (1-0.75,\ 1-0.75,\ 0-0.75,\ 1-0.75) = (0.25,\ 0.25,\ -0.75,\ 0.25).

Part 5: Newton leaf corrections

At probabilities (0.6,0.6), the numerator and denominator are

(1-0.6)+(0-0.6)=-0.2, \qquad 2(0.6)(0.4)=0.48,

so

\widetilde\gamma=-\frac{0.2}{0.48}\approx-0.417.

At probabilities (0.9,0.9),

(1-0.9)+(0-0.9)=-0.8, \qquad 2(0.9)(0.1)=0.18,

and therefore

\widetilde\gamma=-\frac{0.8}{0.18}\approx-4.44.

The pseudo-residuals are bounded, but the Newton step divides their sum by total Bernoulli curvature. When probabilities are near zero or one, that denominator can be small, so the log-odds correction can be large.

12.12 References

Bühlmann, Peter, and Bin Yu. 2003. “Boosting with the L_2 Loss: Regression and Classification.” Journal of the American Statistical Association 98 (462): 324–39. https://doi.org/10.1198/016214503000125.

Chen, Tianqi, and Carlos Guestrin. 2016. “XGBoost: A Scalable Tree Boosting System.” In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, 785–94. https://doi.org/10.1145/2939672.2939785.

Freund, Yoav, and Robert E. Schapire. 1997. “A Decision-Theoretic Generalization of On-Line Learning and an Application to Boosting.” Journal of Computer and System Sciences 55 (1): 119–39. https://doi.org/10.1006/jcss.1997.1504.

Friedman, Jerome H. 2001. “Greedy Function Approximation: A Gradient Boosting Machine.” Annals of Statistics 29 (5): 1189–1232. https://doi.org/10.1214/aos/1013203451.

———. 2002. “Stochastic Gradient Boosting.” Computational Statistics & Data Analysis 38 (4): 367–78. https://doi.org/10.1016/S0167-9473(01)00065-2.

Friedman, Jerome, Trevor Hastie, and Robert Tibshirani. 2000. “Additive Logistic Regression: A Statistical View of Boosting.” Annals of Statistics 28 (2): 337–407. https://doi.org/10.1214/aos/1016218223.

McCracken, Michael W., and Serena Ng. 2020. “FRED-QD: A Quarterly Database for Macroeconomic Research.” Working Paper 26872. National Bureau of Economic Research. https://doi.org/10.3386/w26872.

Yeh, I-Cheng. 2009. “Default of Credit Card Clients.” UCI Machine Learning Repository. https://doi.org/10.24432/C55S3H.

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.

Zhang, Tong, and Bin Yu. 2005. “Boosting with Early Stopping: Convergence and Consistency.” Annals of Statistics 33 (4): 1538–79. https://doi.org/10.1214/009053605000000255.

Footnotes

1. This is a pedagogical final-vintage FRED-QD example. A full real-time GDP forecasting evaluation would also need to account for release calendars and data revisions.