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—halting the boosting run at the iteration count that performs best on held-out data—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 (FNN) using the same data split and the same probability-based 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 \hat F_m for the fitted ensemble prediction function after iteration m and F for a generic candidate function. This chapter-local notation is distinct from the use of F for a predictive cumulative distribution function in earlier chapters. In regression, \hat F_m(x) is a point prediction; in Bernoulli classification below, it is the log-odds of the event.

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—the boosting counterpart of the gradient-descent step size \eta in the optimization chapter; the symbol \nu follows the boosting literature and is unrelated to the degrees-of-freedom parameter of the distributional network chapter
  • 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 Monte Carlo error in the average Trees may remain highly correlated
Random forest Parallel bootstrap trees with feature subsampling max_features and leaf size; B controls Monte Carlo error in the average 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. The training sample is \{(x_i, y_i)\}_{i=1}^N, with predictor vector x_i \in \mathbb{R}^p and scalar outcome y_i—in the book’s case convention, the realized values of random pairs (X_i, Y_i), as in the decision trees chapter. Under squared loss, the training objective of a candidate prediction function F is

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 residuals are negative derivatives of the loss with respect to the prediction. Hold observation i fixed and regard its loss contribution \tfrac{1}{2}(y_i - z)^2 as an ordinary function of one real argument z, the prediction for that observation. Its derivative is z - y_i, so the negative derivative, evaluated at the current fitted value z = \hat F_{m-1}(x_i), is exactly the residual r_{im}. Boosting therefore fits the next tree \hat h_m(x) to the residuals by least squares. The ideal tree would solve

\min_{h}\ \sum_{i=1}^N \bigl(r_{im} - h(x_i)\bigr)^2,

where the minimization runs over all regression trees h of maximum depth d. Solving this problem exactly is computationally infeasible, so in practice \hat h_m is grown by recursive binary splitting with the residuals as the response. Each split is myopically optimal given the splits already made, so the fitted tree aims at a small value of the least-squares criterion without any guarantee of attaining its global minimum over depth-d trees. One part of the fit is exact, though: given the partition the greedy search selects, the prediction in each leaf is the mean of the residuals in that leaf, which is the least-squares-optimal constant there.

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.

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. The data are simulated: a single predictor, an outcome that combines a sine wave, a mild linear trend, and an upward level shift beyond x = 0.8, plus Gaussian noise.

Show the code
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 gains detail 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

Squared loss gives us a particularly convenient starting point. Under \ell(y,z)=\tfrac12(y-z)^2, the negative gradient y-z is an ordinary residual, measured in the same units as the outcome and the prediction. The residual therefore directly describes an additive correction on the target scale. Moreover, its average within a fixed tree leaf is already the loss-minimizing correction for that leaf. These properties make squared-error boosting especially simple.

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). As under squared loss, we start from a constant predictor that minimizes the training loss, assuming a finite minimizer exists.

At iteration m, hold y_i fixed, regard \ell(y_i,z) as a function of the scalar prediction z, and evaluate its negative derivative at the current fitted value:

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

A positive pseudo-residual indicates that a small increase in the prediction lowers the observation’s loss; a negative one indicates that a small decrease does so. Fitting a shallow regression tree to these pseudo-residuals gives a tree approximation to the negative-gradient direction. The update retains the same additive form,

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

For a general loss, the average pseudo-residual in a leaf provides a gradient-based correction but need not minimize the original loss within that leaf. We can therefore refine the leaf values after choosing the splits: choose each leaf’s correction using the original loss, then apply shrinkage (J. H. Friedman 2001).

The loss defines the population risk criterion: the ideal predictor minimizes \mathbb{E}[\ell(Y,F(X))]. Changing the loss can change the statistical target, although different losses can target the same quantity.

Classification with Bernoulli Log-Loss

For binary outcomes y_i\in\{0,1\}, squared-error boosting still targets the conditional mean, which is the event probability:

\mathbb{E}[Y\mid X=x]=\mathbb{P}(Y=1\mid X=x).

The fitted values, however, need not stay in [0,1]. A classification forest preserves this range by averaging leaf probabilities. Boosting adds signed corrections, so a direct probability fit can leave the range even when its initial constant is a valid probability.

In the Bernoulli construction below, we let the ensemble predict unrestricted log-odds, just as the linear index x_i'\beta does in a logit model. The sigmoid converts the ensemble output into a probability:

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

The log-odds, also called the logit, are the logarithm of the odds p_i/(1-p_i):

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

The log-odds can take any real value, and the sigmoid maps every finite updated value into (0,1). We can therefore add tree corrections on the log-odds scale while preserving valid probabilities.

The sigmoid enforces the probability range; the loss is a separate choice. We use the negative log-likelihood of a Bernoulli model. For event probability p_i, this model assigns the observed outcome y_i probability

p_i^{y_i}(1-p_i)^{1-y_i}.

Taking its negative logarithm gives the Bernoulli log-loss:

\ell(y_i,F(x_i)) =-\log\!\left[p_i^{y_i}(1-p_i)^{1-y_i}\right] =-y_i\log p_i-(1-y_i)\log(1-p_i).

Boosting adds corrections to the model’s log-odds output F(x_i), so we need the derivative with respect to this output. Write z for a possible value of F(x_i) and p=\sigma(z) for the corresponding probability. Holding y_i fixed, the sigmoid derivative and the chain rule give

\begin{aligned} \frac{dp}{dz} &=\frac{e^{-z}}{(1+e^{-z})^2}\\ &=p(1-p),\\[0.5em] \frac{\partial\ell(y_i,z)}{\partial z} &=\left(-\frac{y_i}{p}+\frac{1-y_i}{1-p}\right)\frac{dp}{dz}\\ &=\left(-\frac{y_i}{p}+\frac{1-y_i}{1-p}\right)p(1-p)\\ &=-y_i(1-p)+(1-y_i)p\\ &=p-y_i\\ &=\sigma(z)-y_i. \end{aligned}

At iteration m, we evaluate this derivative at the current log-odds z=\hat F_{m-1}(x_i), whose corresponding probability is p_{i,m-1}=\sigma(\hat F_{m-1}(x_i)). The pseudo-residual is therefore

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

The next regression tree thus fits the discrepancy between the binary outcome and the current probability. When both classes occur in the training sample, the loss-minimizing initial probability is \bar y, giving \hat F_0(x)=\log[\bar y/(1-\bar y)].

The distinction between the residual and the update matters: y_i-p_{i,m-1} compares an outcome with a probability, but the tree’s correction is added to the log-odds \hat F_{m-1}. We then recompute the probability using the sigmoid. A common refinement uses a Newton step, which adjusts each leaf correction using the curvature of the Bernoulli loss. The Newton correction also lives on the log-odds scale. Newton’s method is an optional refinement of the gradient step; y_i-p_{i,m-1} is already the negative derivative with respect to log-odds. The optional calculation below gives the details.

Let R_{jm} denote leaf j of the tree at iteration m. Once the tree’s splits are fixed, we can choose a common correction \gamma for all observations in that leaf by minimizing its original loss (J. H. Friedman 2001):

Q_{jm}(\gamma) =\sum_{i:\,x_i\in R_{jm}} \ell\bigl(y_i,\hat F_{m-1}(x_i)+\gamma\bigr), \qquad \gamma_{jm}^{\star}\in\arg\min_\gamma Q_{jm}(\gamma).

For squared loss, the minimizing correction is the mean residual in the leaf. For Bernoulli loss, there is generally no closed-form solution. With finite current log-odds, a finite minimizer exists when both classes occur in the leaf. An all-one or all-zero leaf instead approaches its infimum as \gamma\to+\infty or \gamma\to-\infty, respectively.

A common approximation takes one Newton step from \gamma=0. The first and second derivatives there are

Q'_{jm}(0)=\sum_{i:\,x_i\in R_{jm}}(p_{i,m-1}-y_i), \qquad Q''_{jm}(0)=\sum_{i:\,x_i\in R_{jm}}p_{i,m-1}(1-p_{i,m-1}).

Minimizing the quadratic approximation Q_{jm}(0)+Q'_{jm}(0)\gamma+\tfrac12Q''_{jm}(0)\gamma^2 gives

\widetilde\gamma_{jm} =-\frac{Q'_{jm}(0)}{Q''_{jm}(0)} =\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})}.

The tilde distinguishes this approximate correction from an exact minimizer \gamma_{jm}^{\star}. Using \widetilde\gamma_{jm} as the leaf value in \hat h_m gives, for x_i\in R_{jm},

p_{im}=\sigma\!\left(\hat F_{m-1}(x_i)+\nu\widetilde\gamma_{jm}\right).

Even though each pseudo-residual is bounded, the Newton correction need not be: a non-negligible summed residual divided by very small total curvature produces a large log-odds correction. One Newton step can also overshoot and increase the loss. Shrinkage reduces the step size, but no fixed positive learning rate guarantees a loss decrease for every possible leaf and current fit.

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 and illustrates the same idea with a different loss (Freund and Schapire 1997; J. Friedman, Hastie, and Tibshirani 2000). For this comparison, use labels y_i\in\{-1,+1\} and let F be a real-valued classification score. AdaBoost uses exponential loss,

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

At iteration m, AdaBoost fits a classifier using observation weights

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

The product y_i\hat F_{m-1}(x_i) is the signed margin: it is positive for a correctly classified observation and negative for a misclassified one. A more negative margin gives a larger weight, so the next classifier places more emphasis on the current mistakes. After fitting that classifier, AdaBoost chooses its weight in the ensemble to minimize the exponential training loss. Only this one weight is optimized at that step; the selected classifier stays fixed.

Exponential loss can assign arbitrarily large weights to confidently misclassified observations. For Bernoulli boosting, with its 0/1 labels, the pseudo-residuals instead satisfy |y_i-p_{i,m-1}|<1. That bound limits each observation’s contribution to the gradient at a given iteration; it does not guarantee that the fitted ensemble is insensitive to individual observations.

With y_i\in\{-1,+1\} and exponential loss, the negative gradient is

r_{im} =-\left.\frac{\partial\ell(y_i,z)}{\partial z}\right|_{z=\hat F_{m-1}(x_i)} =y_i\,e^{-y_i\hat F_{m-1}(x_i)}=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}^Nw_{im}y_i h(x_i), \qquad C_m=\sum_{i=1}^N(w_{im}^2+1).

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.

Write \hat g_m(x)\in\{-1,+1\} for the selected classifier and \mathrm{err}_m=\sum_i w_{im}\mathbf{1}\{y_i\ne\hat g_m(x_i)\}/\sum_i w_{im} for its weighted misclassification rate. Holding \hat g_m fixed, the exponential loss after adding \alpha\hat g_m is proportional to

(1-\mathrm{err}_m)e^{-\alpha}+\mathrm{err}_m e^{\alpha}.

Differentiating with respect to \alpha and setting the derivative to zero gives e^{2\alpha}=(1-\mathrm{err}_m)/\mathrm{err}_m, hence

\alpha_m=\frac12\log\!\left(\frac{1-\mathrm{err}_m}{\mathrm{err}_m}\right), \qquad \hat F_m=\hat F_{m-1}+\alpha_m\hat g_m.

A finite positive coefficient requires 0<\mathrm{err}_m<1/2. Error 1/2 gives zero weight; a perfect classifier requires a stopping or finite-step convention. AdaBoost therefore builds the ensemble stage by stage by minimizing exponential loss over each new classifier’s coefficient (J. Friedman, Hastie, and Tibshirani 2000).

The score convention differs from Bernoulli boosting. For p(x)=\mathbb{P}(Y=+1\mid X=x)\in(0,1), minimizing the conditional risk p(x)e^{-F}+(1-p(x))e^F gives F=\tfrac12\log[p(x)/(1-p(x))]. The corresponding probability map is \sigma(2F), whereas Bernoulli boosting uses \sigma(F).

12.6 Regularization in Boosting

Boosting can keep improving the training 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. Along the training-risk path, the interaction between \nu and M is roughly reciprocal: a smaller \nu requires correspondingly more iterations to reach the same training risk, so as a heuristic, halving \nu roughly doubles the number of trees needed. The validation-optimal M is a different matter—if small-\nu paths merely retraced the same fits more slowly, shrinkage could not improve generalization. What J. H. Friedman (2001) reports empirically is the qualitative pattern: smaller \nu tends to deliver better test error and requires a correspondingly larger M; this is a rule of thumb, not a theorem. Zhang and Yu (2005) supply theoretical support for treating the iteration count as a regularization device.

Number of Trees

The number of iterations M controls how long the algorithm keeps adding trees to 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. The mechanism is the one described in the random-forests chapter: candidate splits that would violate the requested ordering are rejected, and bounds on leaf values are propagated down the tree. In boosting, each successive tree is constrained in this way, so the accumulated sum of trees respects the ordering as well.

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, again on simulated data: one predictor, a smooth nonlinearity with a threshold shift, and independent training and validation samples of 280 observations each. The vertical axis is logarithmic, so proportional changes are equally visible everywhere and the behavior near the validation minimum is not compressed by the large initial errors. Because the held-out sample is inspected at every iteration, it is a tuning sample rather than a final test sample.

Show the code
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import NullFormatter, ScalarFormatter
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 (log scale)")
axes[0].set_yscale("log")
axes[0].set_yticks([0.1, 0.2, 0.5, 1, 2])
axes[0].yaxis.set_major_formatter(ScalarFormatter())
axes[0].yaxis.set_minor_formatter(NullFormatter())
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, shown on a logarithmic vertical axis. 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 gives a flatter, more forgiving validation-error profile—so flat here that its minimum has not clearly been reached within the 180 iterations plotted, so the grid’s upper bound rather than the validation curve would determine the stopping point. With the larger learning rate, validation error falls more quickly, reaches its minimum within the first few dozen trees, and then drifts mildly upward.

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?

A small learning rate with many trees is better supported when its validation curve is flatter near the minimum, making performance less sensitive to the chosen stopping point. A larger learning rate with fewer trees is better supported when it reaches similar validation error with fewer boosting steps and lower computational cost. For forecasting, compare validation loss over the relevant forecast origins; for prediction with independent and identically distributed (i.i.d.) observations, use an appropriate independent validation sample. Training fit alone does not settle the comparison.

Real macroeconomic data can be less forgiving. Figure 12.3 shows training and validation MSE curves—again on a logarithmic vertical axis—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 of the 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—forecast origins 1962Q3 through 2009Q2—and evaluated, without re-estimation, on the remaining 25%, origins 2009Q3 through 2025Q1.1 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. On the logarithmic axis this deterioration looks mild—about 6% of the zero-tree level at \nu=0.02 and 16% at \nu=0.2—while the training MSE falls to roughly half and a tenth of its starting value, respectively; the widening gap between the two curves is the clearest visual sign that the added flexibility buys nothing out of sample.

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 model fitted on the training quarters fixed rather than re-estimating it; a recursive or rolling design that refits at each forecast origin estimates a different deployment risk and could select a different stopping point. The result is not driven solely by the four 2020 outcomes: dropping those validation targets leaves the minimum of both validation curves at zero trees. Those quarters strongly affect the level of the validation MSE—the 2020 contraction and rebound generate residuals of tens of annualized percentage points, and removing those quarters cuts the zero-tree validation MSE by an order of magnitude. Under this fixed-estimation design, these predictors and boosting settings yield no validation gain over the constant. Weak usable signal, instability in the predictive relationship, and estimation error could all contribute; the exclusion check does not distinguish them. In the synthetic example above, the data-generating process is known to be stable and to contain nonlinear signal, and validation error improves substantially on the constant model before flattening out.

Show the code
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.ticker import NullFormatter, ScalarFormatter
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 (log scale)")
axes[0].set_yscale("log")
axes[0].set_yticks([1, 2, 5, 10, 20, 40])
axes[0].yaxis.set_major_formatter(ScalarFormatter())
axes[0].yaxis.set_minor_formatter(NullFormatter())
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, shown on a logarithmic vertical axis. 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—folds built so that all observations from the same group (borrower, firm, region) stay together
  • leakage-safe preprocessing inside each training fold

Once validation selects the learning rate, depth, and stopping point, assess the chosen procedure on an untouched test sample or a subsequent forecast period. The smallest loss on the tuning path is not an independent estimate of the selected model’s performance.

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 records for 30,000 clients, including 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 and random-forest settings are chosen by validation log loss from small candidate sets declared in the code. For boosting, the learning rate and depth come from a small grid, and the number of trees is chosen by early stopping: the validation log loss is tracked from the initial constant prediction through each grid point’s staged predictions, and the iteration with the smallest loss is kept. Exact ties along a path favor fewer trees. If zero trees wins, we retain a constant-probability model and re-estimate its probability on the combined training and validation samples before testing. The iteration count is thus tuned jointly with the learning rate. The comparison uses the same transformed predictors for every model: continuous variables are standardized and categorical variables are one-hot encoded. Standardization helps numerical optimization for 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. One consequence of the encoding is not incidental, though: the repayment-status codes are ordered categories (months of payment delay), and one-hot encoding discards the ordering, so no single split can implement a threshold rule such as “delay of two months or more”; the tree-based models must assemble such rules from several level-versus-rest splits. We accept this to keep one representation for all models; an ordinal encoding could well favor the tree-based methods. 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. Log loss and the Brier score evaluate the probabilities themselves; accuracy evaluates one particular binary decision rule.

Show the code
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.dummy import DummyClassifier
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]
]

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

# Boosting: for each (rate, depth) grid point, the number of trees is
# chosen by early stopping, including the initial constant at zero trees.
initial_validation_loss = log_loss(
    y_validation, np.full(len(y_validation), y_train.mean())
)
boosting_results = []
for rate in [0.1, 0.05]:
    for depth in [1, 2, 3, 4]:
        path_model = GradientBoostingClassifier(
            n_estimators=400,
            learning_rate=rate,
            max_depth=depth,
            min_samples_leaf=20,
            random_state=42,
        )
        path_model.fit(X_train_transformed, y_train)
        staged_losses = [initial_validation_loss]
        staged_losses.extend(
            log_loss(y_validation, staged[:, 1])
            for staged in path_model.staged_predict_proba(
                X_validation_transformed
            )
        )
        # Index zero is the constant; argmin favors fewer trees on exact ties.
        best_m = int(np.argmin(staged_losses))
        boosting_results.append((min(staged_losses), rate, depth, best_m))

_, boosting_rate, boosting_depth, boosting_m = min(boosting_results)
if boosting_m == 0:
    boosting_setting = "0 trees (constant probability)"
    # fit() below re-estimates the class probability on the development sample.
    selected_boosting = DummyClassifier(strategy="prior")
else:
    boosting_setting = f"rate {boosting_rate}, depth {boosting_depth}, {boosting_m} trees"
    selected_boosting = GradientBoostingClassifier(
        n_estimators=boosting_m,
        learning_rate=boosting_rate,
        max_depth=boosting_depth,
        min_samples_leaf=20,
        random_state=42,
    )

# 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 the three tree-based methods are chosen by validation log loss; for gradient boosting, the number of trees is additionally chosen by early stopping along the validation path. 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, depth 4, 91 trees 0.4311 0.1354 0.8178

Gradient boosting attains the lowest test log loss, 0.4311, and Brier score, 0.1354, in this comparison, with the random forest close behind on both probability scores and the selected FNN third. The single tree improves substantially on the constant-probability benchmark but has the weakest probability scores among the fitted models. The margin between the two tree ensembles is numerically small; no paired uncertainty estimate is reported here. Accuracy gives a less informative and slightly different ranking: the forest attains 81.80%, boosting 81.78%, and the FNN and logistic regression both attain 81.70%.

These results support a deliberately local conclusion. On this split, both tree ensembles improve on a single tree and outperform the FNN on the probability scores, while the table alone does not establish which ensemble has lower expected loss. The table does not show that boosting universally dominates forests or networks. 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.

One data caveat applies to every model in the table alike. The dataset documentation does not pin down precisely how the most recent repayment-status code is timed relative to the outcome month; if that code partly reflects payment behavior in the outcome month itself, every model receives an informative signal that a real-time scoring system would not yet have. Such leakage could make the reported scores optimistic and also alter the relative ranking, because the models can exploit the same unavailable signal to different degrees. The table therefore compares methods on a fixed information set; it does not establish deployable performance.

Question for Reflection

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

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.

12.9 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.
  • Bernoulli pseudo-residuals are y-p, but leaf corrections are added to the log-odds F, which the sigmoid then converts into a probability.
  • Adding more trees does not give the ensemble an ability to extrapolate smooth trends outside the training support.

12.10 Exercises

Exercise 12.1: First Step of Gradient Boosting with Shrinkage

Consider four monthly forecast-origin/target pairs, indexed by the forecast origins t=1,\dots,4. 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).

In this small sample the indicator values happen to coincide with the time index; the stump below splits on the indicator value x, not on time. Consider squared-error boosting with decision stumps as weak learners. Let the learning rate be \nu = 0.4.

  1. Compute the initial model \hat F_0(x) and the residuals r_{t+1} = y_{t+1} - \hat F_0(x_t), t=1,\dots,4, and explain why these residuals are exactly the negative gradients of the squared-error training loss \frac{1}{2}\sum_{t=1}^{4} \bigl(y_{t+1} - F(x_t)\bigr)^2 with respect to the fitted values.
  2. Fit the optimal first stump \hat h_1(x) by checking the candidate split points 1.5, 2.5, and 3.5.
  3. Write the updated model \hat F_1(x) = \hat F_0(x) + 0.4\, \hat h_1(x). Compute the fitted values at the four training points.
  4. The stump from Part 2 is fitted by least squares, so its leaf values are the leaf means of the residuals. Show that the training sum of squared errors (SSE) after one shrunken update satisfies \mathrm{SSE}(\nu) = \sum_{t=1}^{4} r_{t+1}^2 \;-\; \nu(2-\nu)\sum_{t=1}^{4} \hat h_1(x_t)^2, and verify the formula at \nu=0.4 and \nu=1 by comparing with the directly computed SSE. Conclude that on the training sample, the full step \nu=1 always fits at least as well as any shrunken step—and explain why this does not make \nu=1 the right choice out of sample.
  5. For this part only, replace the scalar predictor by two binary economic indicators x_1,x_2\in\{0,1\}, with all four combinations possible. Suppose the conditional mean of next-month output growth is f(x_1,x_2)=x_1x_2. Consider a stump ensemble \hat F_M=\hat F_0+\nu\sum_{m=1}^{M}\hat h_m with a constant initial model \hat F_0, a learning rate \nu>0, and decision stumps \hat h_1,\dots,\hat h_M. Show that it can be written as \hat F_M(x_1,x_2)=a+g_1(x_1)+g_2(x_2), where a is a constant, g_1 is a function of x_1 alone, and g_2 is a function of x_2 alone. Compute the contrast F(1,1)-F(1,0)-F(0,1)+F(0,0) both for F=\hat F_M and for F=f. Can increasing M or changing \nu make a stump ensemble match all four conditional means? Construct a depth-2 tree that represents f exactly.

Exam level: suitable as-is. Parts 1–4 combine boosting mechanics with a projection argument for shrinkage, and Part 5 isolates an interaction restriction that additional stumps cannot remove.

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

Expand \sum_{t=1}^{4} \bigl(r_{t+1} - \nu \hat h_1(x_t)\bigr)^2. For the cross term, use that within each leaf \hat h_1 is constant and equal to the mean of the residuals it fits, so \sum_{t=1}^{4} r_{t+1}\hat h_1(x_t) = \sum_{t=1}^{4} \hat h_1(x_t)^2.

Each stump splits on only one of the two indicators. Group the stumps by that indicator, then substitute the additive form into the contrast.

Part 1: Initial model, residuals, and the gradient view

The initial model is the sample mean:

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

Therefore the residuals are

(r_2,\,r_3,\,r_4,\,r_5) = (-1.125,\,-0.125,\,-0.625,\,1.875).

These are negative gradients because differentiating the training loss \frac{1}{2}\sum_{t=1}^{4} (y_{t+1}-F(x_t))^2 with respect to the fitted value F(x_t) gives -(y_{t+1}-F(x_t)); evaluating at the current fit F=\hat F_0 and flipping the sign yields exactly the residual r_{t+1}. Fitting the next stump to the residuals is therefore fitting it to the steepest-descent direction of the training loss.

Part 2: 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

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

Part 3: Shrinkage update

With \nu=0.4,

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

Therefore

\hat 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 4: The effect of shrinkage on the training fit

After the update, the new residuals are r_{t+1}-\nu\hat h_1(x_t), so

\mathrm{SSE}(\nu) =\sum_{t=1}^{4}\bigl(r_{t+1}-\nu\hat h_1(x_t)\bigr)^2 =\sum_{t=1}^{4} r_{t+1}^2-2\nu\sum_{t=1}^{4} r_{t+1}\hat h_1(x_t)+\nu^2\sum_{t=1}^{4}\hat h_1(x_t)^2.

Denote the two leaves of the stump by R_{11} and R_{21}. On leaf R_{j1}, \hat h_1 is constant and equal to the mean residual \gamma_{j1}^{\star} in that leaf, so \sum_{t:\,x_t\in R_{j1}} r_{t+1}=\sum_{t:\,x_t\in R_{j1}}\gamma_{j1}^{\star}. Summing leaf by leaf therefore gives

\sum_{t=1}^{4} r_{t+1}\hat h_1(x_t) =\sum_{j=1}^{2}\gamma_{j1}^{\star}\sum_{t:\,x_t\in R_{j1}}r_{t+1} =\sum_{j=1}^{2}\sum_{t:\,x_t\in R_{j1}}\bigl(\gamma_{j1}^{\star}\bigr)^2 =\sum_{t=1}^{4}\hat h_1(x_t)^2,

and substituting this identity into the expansion yields

\mathrm{SSE}(\nu)=\sum_{t=1}^{4} r_{t+1}^2-\nu(2-\nu)\sum_{t=1}^{4}\hat h_1(x_t)^2.

Numerically, \sum_{t=1}^{4} r_{t+1}^2=1.125^2+0.125^2+0.625^2+1.875^2=5.1875 and \sum_{t=1}^{4}\hat h_1(x_t)^2=3\times 0.625^2+1.875^2=4.6875. At \nu=0.4 the formula gives 5.1875-0.4\times 1.6\times 4.6875=2.1875, which matches the direct computation using the fitted values from Part 3:

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

At \nu=1 the formula gives 5.1875-4.6875=0.5, again matching the direct computation with fitted values (2.5,\,2.5,\,2.5,\,5):

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

Because \nu(2-\nu) is strictly increasing on (0,1], the training SSE after one update is strictly decreasing in \nu: the full step always fits the training sample at least as well as any shrunken step. Out of sample this is beside the point. Boosting is a multi-step procedure, the learning rate regularizes the whole path of updates, and the combination of \nu and the stopping iteration M is chosen on validation data; cautious updates often generalize better than aggressive early corrections even though they lose the one-step training comparison by construction.

Part 5: An interaction that stumps cannot represent

Each stump has a single split, on either x_1 or x_2, and is therefore a function of that indicator alone. Grouping the terms in \hat F_M=\hat F_0+\nu\sum_{m=1}^M\hat h_m by their split variable gives

\hat F_M(x_1,x_2)=a+g_1(x_1)+g_2(x_2),

with a=\hat F_0, g_1(x_1) equal to \nu times the sum of the stumps that split on x_1, and g_2(x_2) equal to \nu times the sum of the stumps that split on x_2; if no stump splits on an indicator, the corresponding function is zero. Therefore

\begin{aligned} &\hat F_M(1,1)-\hat F_M(1,0)-\hat F_M(0,1)+\hat F_M(0,0)\\ &=[a+g_1(1)+g_2(1)]-[a+g_1(1)+g_2(0)]\\ &\quad-[a+g_1(0)+g_2(1)]+[a+g_1(0)+g_2(0)] =0. \end{aligned}

For f(x_1,x_2)=x_1x_2, the same contrast is 1-0-0+0=1. Thus no stump ensemble can match all four conditional means, regardless of the number of iterations or the learning rate. More iterations cannot remove this restriction on interactions.

A depth-2 tree can represent f: split first at x_1=0.5 and assign the left leaf the value 0. In the right branch, split at x_2=0.5, assigning values 0 and 1 to its left and right leaves. The tree predicts 1 exactly when both indicators equal one.

Exercise 12.2: Bernoulli Log-Loss and Initialization

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

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

The per-observation Bernoulli log-loss is

\ell(y_i,z)=-y_i\log\sigma(z)-(1-y_i)\log(1-\sigma(z)).

  1. Gradient and pseudo-residuals. Derive the derivative of \ell(y_i,z) with respect to z and evaluate it at z=F(x_i). Hence obtain the pseudo-residual r_i. State whether the tree’s correction is added to the probability or to the log-odds.
  2. Initialization and numerical example. Suppose both classes occur in the training sample and the initial model is a constant, \hat F_0(x)\equiv c. Derive the value of c minimizing total log-loss and show that it is the unique minimizer. For the sample y=(1,1,0,1), compute \bar y, the optimal \hat F_0, the implied initial probability for each observation, and the pseudo-residuals. Exact logarithms are acceptable.

Exam level: suitable as-is. The two parts derive and apply the classification gradient and initialization.

A first-order condition identifies a candidate constant. To establish uniqueness, examine the second derivative of the total loss with respect to c.

Part 1: Gradient and pseudo-residuals

Hold y_i fixed and write p=\sigma(z). Since dp/dz=p(1-p), the chain rule gives

\begin{aligned} \frac{\partial\ell(y_i,z)}{\partial z} &=\left(-\frac{y_i}{p}+\frac{1-y_i}{1-p}\right)p(1-p)\\ &=-y_i(1-p)+(1-y_i)p\\ &=p-y_i. \end{aligned}

At the current log-odds z=F(x_i), the corresponding probability is p_i. The pseudo-residual is the negative derivative,

\begin{aligned} r_i &=-\left.\frac{\partial\ell(y_i,z)}{\partial z}\right|_{z=F(x_i)}\\ &=y_i-p_i. \end{aligned}

The tree’s correction is added to the log-odds F(x_i). The sigmoid then converts the updated log-odds into a new probability.

Part 2: Optimal constant and numerical example

For a constant model, all observations have the same probability p=\sigma(c). Write L(c)=\sum_{i=1}^N\ell(y_i,c) for its total loss. Using Part 1,

\begin{aligned} L'(c) &=\sum_{i=1}^N(p-y_i)\\ &=N(p-\bar y). \end{aligned}

The first-order condition therefore requires p=\bar y. Because both classes occur, 0<\bar y<1, so the corresponding log-odds are finite:

\hat F_0(x)=\log\left(\frac{\bar y}{1-\bar y}\right).

For every finite c, the sigmoid probability lies in (0,1), and

L''(c)=Np(1-p)>0.

Thus L is strictly convex, and the stationary point is its unique global minimizer.

For y=(1,1,0,1), the sample mean is \bar y=3/4, so

\begin{aligned} \hat F_0(x) &=\log\left(\frac{3/4}{1/4}\right)\\ &=\log 3. \end{aligned}

Every initial probability is p_i=3/4. The pseudo-residuals are

\begin{aligned} (r_1,r_2,r_3,r_4) &=(1-3/4,\;1-3/4,\;0-3/4,\;1-3/4)\\ &=(1/4,\;1/4,\;-3/4,\;1/4). \end{aligned}

12.11 References

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.↩︎