Building models with scikit-learn

scikit-learn is where machine learning on a table of numbers happens in Python. Linear models, trees, forests, boosting, support-vector machines, clustering, dimensionality reduction, and around them the parts that matter more than the algorithm: preprocessing, cross-validation, hyperparameter search, metrics. It began as a Google Summer of Code project in 2007 and became what it is when a group at INRIA in France took it over in 2010 (about). The Python stack post placed it among its neighbours; this post is about using it well.

The example is the one the Machine Learning Series 1 used: five tickers, daily features computed from prices (recent returns, recent volatility, distance from the moving average), and the question will the next 20 trading days be a stress month, volatility above 30% a year. Train on 2020 to 2023, test on 2024 to 2025. Every number below comes from one script on that data.

The one idea: everything is an estimator

scikit-learn has one interface, and once you know it you know the library. Every object is an estimator: you construct it with its settings (hyperparameters), you call fit(X, y) to learn from data, and afterwards it has learned attributes whose names end in an underscore (coef_, feature_importances_, classes_). Beyond that, three kinds (Table 1):

The design was a deliberate choice, written up by the core developers in 2013 (Buitinck et al.): a shared interface for every unit, composition over inheritance, and sensible defaults so that a model works before it is tuned. Géron’s Hands-On Machine Learning summarises the same principles as consistency, inspection, non-proliferation of classes, composition and sensible defaults. The practical consequence is that swapping a model is a one-line change, and the code around it (splitting, scaling, scoring, searching) does not care which model it holds.

The code below runs as pasted, in order, in one Python session (pip install numpy pandas scikit-learn). This first block builds the data from the blog’s price file. Each row is one ticker on one trading day; X_train holds the seven features and the ticker, y_train the answer (1 if a stress month followed); X_test and y_test are the same for 2024 to 2025.

import numpy as np
import pandas as pd

URL = "https://raw.githubusercontent.com/ankurtrapasiya/examples/main/data/prices.csv"
prices = pd.read_csv(URL, index_col="date", parse_dates=True)

FEATURES = ["ret_1d", "ret_5d", "ret_20d", "vol_5d", "vol_20d", "sma_ratio", "bollinger"]
rows = []
for ticker in prices.columns:                    # one block of rows per ticker
    p = prices[ticker]
    r = p.pct_change()                           # daily returns
    sma = p.rolling(20).mean()                   # 20-day average price
    rows.append(pd.DataFrame({
        "ret_1d": r, "ret_5d": p / p.shift(5) - 1, "ret_20d": p / p.shift(20) - 1,
        "vol_5d": r.rolling(5).std(), "vol_20d": r.rolling(20).std(),
        "sma_ratio": p / sma - 1, "bollinger": (p - sma) / (2 * p.rolling(20).std()),
        "next_vol": r[::-1].rolling(20).std()[::-1].shift(-1),   # volatility of the NEXT 20 days
        "ticker": ticker,
    }).dropna())
df = pd.concat(rows).rename_axis("date").reset_index().sort_values(["date", "ticker"])
df["stress"] = (df["next_vol"] * np.sqrt(252) * 100 > 30).astype(int)   # 1 = stress month ahead

train, test = df[df["date"] < "2024-01-01"], df[df["date"] >= "2024-01-01"]   # time order
X_train, y_train = train[FEATURES + ["ticker"]], train["stress"]
X_test, y_test = test[FEATURES + ["ticker"]], test["stress"]
print(len(X_train), len(X_test), round(y_train.mean(), 3), round(y_test.mean(), 3))   # 4930 2205 0.368 0.241

It prints 4930 2205 0.368 0.241: the training and test rows, and the share of each with a stress month ahead. r[::-1].rolling(20).std()[::-1].shift(-1) turns the series round, takes the 20-day volatility, and turns it back, which gives the volatility of the next 20 days. Now the four models:

from sklearn.compose import ColumnTransformer
from sklearn.dummy import DummyClassifier
from sklearn.ensemble import HistGradientBoostingClassifier, RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import OneHotEncoder, StandardScaler

def prep():                       # a fresh preprocessor per pipeline (see below)
    return ColumnTransformer([("num", StandardScaler(), FEATURES),
                              ("ticker", OneHotEncoder(handle_unknown="ignore"), ["ticker"])])

models = {
    "always 'calm'": DummyClassifier(strategy="prior"),
    "logistic regression": Pipeline([
        ("prep", prep()), ("model", LogisticRegression(max_iter=1000))]),
    "random forest": Pipeline([
        ("prep", prep()), ("model", RandomForestClassifier(
            n_estimators=300, max_features=3, min_samples_leaf=20, random_state=0))]),
    "gradient boosting": Pipeline([
        ("prep", prep()), ("model", HistGradientBoostingClassifier(
            max_depth=3, learning_rate=0.05, max_iter=200, random_state=0))]),
}
for name, m in models.items():
    m.fit(X_train, y_train)                       # the same two calls for every model
    p = m.predict_proba(X_test)[:, 1]             # chance of a stress month, per test row
    print(f"{name:20} accuracy {m.score(X_test, y_test):.3f}  AUC {roc_auc_score(y_test, p):.3f}")

Four models, one loop. Two things in Table 2 before going further. First, the dummy classifier that always answers “calm” scores 76% accuracy, because 76% of test rows have a calm month ahead; accuracy is the wrong score for this problem, and the next section says what to use. Second, the plain logistic regression sits within 0.003 of AUC of the random forest, leads all three on average precision, and loses to the tree models on accuracy by 0.003 and 0.004. A linear model with scaled inputs is the baseline to beat. Do not skip it.

Preprocessing belongs inside the pipeline

The ColumnTransformer above scales the seven numeric columns and one-hot encodes the ticker, and the Pipeline glues that to the model. The glue is not there for tidiness. A transformer fitted on the whole dataset has seen the test set. StandardScaler learns a mean and a standard deviation; fit it on all the rows and the test rows’ statistics leak into training. Inside a pipeline, fit runs the scaler on the training rows only and predict reuses those numbers on the test rows, and the same holds inside every cross-validation fold. The same reasoning covers SimpleImputer (the fill value for missing data is learned) and OneHotEncoder (the set of categories is learned; handle_unknown="ignore" is what keeps a ticker that was not in training from crashing prediction).

Two things the library does not do for you:

• Estimators are stateful, and a pipeline fits its steps in place. The first version of the script built one ColumnTransformer and passed the same object to four pipelines. Each fit refitted it, which was harmless while every pipeline saw the same columns, and then a fifth pipeline fed it rescaled columns and silently broke the other four: their test AUC dropped a dozen points with no error. Build a fresh object per pipeline, or use sklearn.base.clone. Meta-estimators such as GridSearchCV and cross_val_score clone for you; Pipeline does not.
• Custom steps are a small class. Anything with fit and transform can go in a pipeline; inheriting from BaseEstimator and TransformerMixin gives it fit_transform, get_params and a sensible repr for free. Here is one that turns daily volatility into the yearly percentage the rest of the site quotes, so that the numbers mean something when the model is inspected:

from sklearn.base import BaseEstimator, TransformerMixin
from sklearn.pipeline import make_pipeline

class Annualise(BaseEstimator, TransformerMixin):
    def __init__(self, columns=("vol_5d", "vol_20d")):
        self.columns = columns        # stored as given, never changed

    def fit(self, X, y=None):
        return self                   # nothing to learn

    def transform(self, X):
        X = X.copy()                  # never change the caller's table
        X[list(self.columns)] = X[list(self.columns)] * (252 ** 0.5) * 100   # daily -> % a year
        return X

pipe = make_pipeline(Annualise(), prep(), LogisticRegression(max_iter=1000))
pipe.fit(X_train, y_train)
print(round(roc_auc_score(y_test, pipe.predict_proba(X_test)[:, 1]), 3))   # 0.890, as before

It prints 0.890. The constructor stores its arguments untouched; that is the one rule the library’s developer guide insists on, because get_params and set_params (which grid search uses) read and write those attributes by name. The model with this step scores the same $\key{0.890}$ as before, since the scaler undoes a constant factor; the step is there so that vol_20d reads as “31% a year” rather than “0.0195” when you look inside.

Scores: accuracy is the wrong one here

The confusion matrix of the logistic regression on the test rows, at the default cut-off of 0.5 (Table 3):

TP, stress cases caughtFN, stress cases missedFP, false alarms

$$\text{precision} = \frac{\good{TP}}{\good{TP} + \risk{FP}} = \frac{387}{387 + 141} = 0.73$$ $$\text{recall} = \frac{\good{TP}}{\good{TP} + \bad{FN}} = \frac{387}{387 + 145} = 0.73$$

When the model says “stress” it is right $\key{73\%}$ of the time, and it catches $\key{73\%}$ of the stress months. The 0.5 cut-off is a default, not a decision. A risk desk that wants to catch most stress months accepts more false alarms and moves the cut-off down. class_weight="balanced" does that at training time by weighting the rare class up: recall rises to $\key{76\%}$ (404 of 532 caught) and precision falls to $\key{69\%}$, with the AUC unchanged (0.890 both ways), because AUC measures the ranking of the probabilities, not the cut-off. The move is small on this data because stress months are 37% of the training rows, so the classes are not far from balanced before any weighting. That is why the tables here report AUC and average precision, not accuracy alone: both score how well the model separates the classes, whatever cut-off is chosen later. classification_report prints precision, recall and F1 per class in one call; use it before quoting an accuracy.

Cross-validation: do not shuffle time

Cross-validation is how a score is estimated without touching the test set: split the training rows into folds, fit on all but one, score on the held-out one, repeat. The default splitter, KFold with shuffling, assigns rows to folds at random. On a time series that is a leak. Rows from the same week share almost everything (the same market, nearly the same 20-day volatility), so a model validated on a random Tuesday has trained on the Monday and Wednesday around it. The honest splitter is TimeSeriesSplit: every fold trains on the past and validates on the block that follows.

from sklearn.model_selection import KFold, TimeSeriesSplit, cross_val_score

gb = models["gradient boosting"]
shuffled = KFold(5, shuffle=True, random_state=0)                  # random folds: leaks
for name, cv in [("shuffled", shuffled), ("time order", TimeSeriesSplit(5))]:
    scores = cross_val_score(gb, X_train, y_train, cv=cv, scoring="roc_auc")
    print(f"{name:10} folds {scores.round(3)}  mean {scores.mean():.3f}")   # 0.929 / 0.899

It prints each splitter’s five fold scores and their mean, the rows of Table 4:

Panel 1 of Figure 1 shows the two. The shuffled folds agree with each other and promise $\key{0.93}$; the model then delivers $\key{0.89}$ on real future data. The time-ordered folds promise less, $\key{0.90}$, and spread more (0.86 to 0.96), and that is useful information: the lowest fold, the third, trained on 2020 and 2021 and was asked about January to September 2022, the year the market fell without a crash; the highest, the fifth, was asked about the calm May to December of 2023. No fold is a disaster on this window, because the March 2020 crash sits inside the first training block, so every fold has seen one. The gap between the two means is still the shuffled splitter’s optimism, and it is the shuffled number that is wrong.

Note. The rule is older than scikit-learn: validate on data that comes after the data you trained on, the way the real test will. The library will happily shuffle time series for you, because KFold has no idea the rows are dated. Choosing the splitter is your job, and it is the single setting that most changes what a cross-validated score means.

Grid search, on the right splitter

GridSearchCV takes an estimator, a dictionary of settings to try, a splitter and a scoring rule, fits every combination on every fold, and keeps the best. The model__learning_rate syntax (double underscore) reaches into a pipeline step by name, which is why pipelines and grid search fit together. The search below tries 12 combinations for the boosting model, and a second search tries 11 values of the regularisation strength C for the logistic regression.

from sklearn.model_selection import GridSearchCV

settings = {"model__learning_rate": [0.01, 0.03, 0.1, 0.3],       # 4 x 3 = 12 combinations
            "model__max_depth": [2, 3, 5]}
grid = GridSearchCV(gb, settings, cv=TimeSeriesSplit(5), scoring="roc_auc", n_jobs=-1)
grid.fit(X_train, y_train)
print(grid.best_params_, round(grid.best_score_, 3))               # chosen on validation folds only
p = grid.best_estimator_.predict_proba(X_test)[:, 1]
print(round(roc_auc_score(y_test, p), 3))                          # looked at once, after the choice

It prints the chosen settings with their validation AUC, then the test AUC, the first row of Table 5 (the logistic regression search is the same pattern with model__C):

Panel 2 of Figure 1 shows the boosting grid: the slowest learning rate won, with shallow trees (depth 3, and depth 2 a hair behind), which is the usual answer on small noisy data, and the whole surface spans 0.886 to 0.910 with fold standard deviations around 0.03, so the differences between neighbours are inside the noise. The tuned boosting model reaches $\key{0.897}$ on the test set, 0.006 above untuned; the tuned logistic regression reaches $\key{0.891}$, 0.001 above untuned. Tuning bought almost nothing here, and the grid says why: every cell is within one standard deviation of every other. The test numbers were looked at once, after the choice was made. That is what makes them test numbers. cv_results_ holds every combination’s mean and spread, and reading it is more useful than reading best_params_.

Two habits from the same section of the library. RandomizedSearchCV samples settings instead of trying all of them, and usually finds as good an answer in a fraction of the fits once the grid has more than two or three dimensions. And n_jobs=-1 runs the folds on every core; scikit-learn has no GPU support, and parallel folds is where its speed comes from.

Inspection: what did it learn?

The learned attributes with the trailing underscore are the first place to look (coef_ for a linear model, feature_importances_ for trees). The better tool is permutation importance: shuffle one column of the test set, score again, and see how much was lost. For the tuned boosting model, shuffling the 20-day volatility costs $\risk{0.152}$ of AUC and shuffling the ticker costs $\risk{0.054}$; no other column costs more than $\risk{0.002}$. The model is close to a function of two things: this month’s volatility predicts next month’s, and the name does too, because NVDA is above the 30% line in 83% of its test months and SPY in 5%, so the one-hot ticker column is a volatility level in disguise. The random forests post uses the same tool. Panel 3 of Figure 1 shows the three models’ ROC curves on the test rows, close enough that the choice between them comes down to simplicity.

Saving the model

pickle and joblib write any fitted estimator to disk, and the library’s persistence guide gives the two warnings that matter. There is no supported way to load a model into a different version of scikit-learn, so record the version next to the file. And a pickle executes code when loaded, so never load one from a source you do not trust; skops is the format for that case, and ONNX for serving predictions without Python. Save the training script, the data snapshot and the cross-validation score with the model. Without them the file is a number with no history.

Where it stops

scikit-learn works on a table that fits in memory, on CPU, with classical methods. That covers more real problems than the attention paid to the alternatives suggests. The Python stack post says where to go when it is not enough: XGBoost and LightGBM for the last few points of accuracy on a table (both plug into pipelines and grid search unchanged), PyTorch when the input is images, text or audio, Dask or Spark when the table does not fit one machine. Two more limits are not in the manual: there is no time-series forecasting in the library (you build the lagged features yourself, as this example did, and then it is a table), and the defaults are tuned for independent rows, so every splitter, every shuffle and every importance measure needs a second look when the rows are dated.

Why banks use it

• The interface is auditable: a pipeline is a list of named steps with stored settings, and get_params() prints the whole thing. A model-risk reviewer can read it.
• Determinism: random_state everywhere, so a result reproduces.
• Linear models and trees, the methods a regulator will accept an explanation for, are first-class, with the inspection tools beside them.
• It runs on the machine the data is already on. No cluster, no GPU, no new infrastructure to approve.

The script behind all the numbers

Every table and the chart come from examples/scikit-learn/stress_month.py, which also holds the confusion matrix, the balanced class weights, the logistic regression search and the permutation importance. To run it yourself, clone the examples repo: its README has a three-command Docker setup, or run uv run scikit-learn/stress_month.py.

Summary

Table 6 is the list of habits.

References

• Aurélien Géron, Hands-On Machine Learning with Scikit-Learn, Keras, and TensorFlow (2nd ed., O’Reilly, 2019), ch. 2 (the estimator design, transformation pipelines, custom transformers, grid search) and ch. 3 (confusion matrix, precision and recall, ROC curves). The author’s companion notebooks are at github.com/ageron/handson-mlp.
• James, Witten, Hastie, Tibshirani & Taylor, An Introduction to Statistical Learning, with Applications in Python (Springer, 2023), 5.1 on cross-validation.
• Lars Buitinck et al., “API design for machine learning software: experiences from the scikit-learn project” (2013), arXiv:1309.0238.
• scikit-learn user guide, the four pages this post leans on: pipelines and composite estimators, cross-validation, developing estimators and model persistence.