Workshop 4 - SOLUTION, Algorithms and Data Analysis

Introduction to Machine Learning: the logistic regression model

Abstract

This is an INDIVIDUAL workshop. We build our first Machine Learning model: a logistic regression that predicts whether a stock will beat the median return of its industry in the next year, using accounting information of the current year. We learn the logic of the logistic model (probabilities, odds and log-odds), how to interpret its coefficients, and the Machine Learning workflow: training and testing samples split by time, predictions, the confusion matrix, precision, sensitivity, specificity, the ROC curve, and how to evaluate whether a model is useful for an investor. Finally, we program the logistic regression ourselves with an optimization algorithm (gradient descent).

1 General directions

Work in a new Google Colab notebook named W4-Algorithms-YourFirstName-YourLastName. Same rules as before: re-type and run the code, write your notes, predict before you run, and verify your results.

You do not have to submit this workshop; however, by doing this workshop in detail you will be prepared for the final evidence and the final exam.

This workshop is the best preparation for the Evidence (in-class exam).

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.formula.api as smf

2 What is Machine Learning?

Machine Learning (ML) is a set of methods that learn patterns from data to make predictions about new data. There are two main families:

• Supervised learning: we have examples with the correct answer (the target or dependent variable Y), and the model learns the relationship between the features (explanatory variables X) and Y. If Y is a number (a return, a price), it is a regression problem; if Y is a category (default / no default, buy / do not buy), it is a classification problem.
• Unsupervised learning: there is no target; the model finds structure in the data (for example, clusters of similar firms).

Classification problems are everywhere in Finance: Will this client default on the loan? Is this transaction a fraud? Will this stock beat its industry next year? The logistic regression is the simplest and most widely used classification model, and it is the foundation of more complex models such as neural networks.

What makes ML different from the regression models you know is not the model itself, but the workflow: we care about how well the model predicts data that it has never seen. Then, we always estimate the model with one part of the data (the training sample) and evaluate it with another part (the testing sample).

3 From linear regression to logistic regression

Imagine a binary target Y that is 1 if the stock beats its industry next year and 0 otherwise. Could we use a linear regression, Y=b_0+b_1X? The prediction of a linear model can be any number, like −0.3 or 1.4, which cannot be a probability. We need a function that transforms any number into a value between 0 and 1. That function is the sigmoid (or logistic) function:

P(Y=1)=\frac{1}{1+e^{-z}}\qquad\text{where}\qquad z=b_0+b_1X_1+b_2X_2+\dots

z = np.linspace(-6, 6, 200)
sigmoid = 1 / (1 + np.exp(-z))
plt.figure(figsize=(7, 4))
plt.plot(z, sigmoid)
plt.axhline(0.5, color="grey", lw=0.7); plt.axvline(0, color="grey", lw=0.7)
plt.xlabel("z = b0 + b1 X1 + ..."); plt.ylabel("P(Y=1)")
plt.title("The sigmoid function transforms any number into a probability")
plt.grid(alpha=0.3); plt.show()

When z=0, the probability is 0.5. When z is large and positive, the probability approaches 1; when it is large and negative, it approaches 0.

3.1 Probability, odds and log-odds

The odds of an event are the probability that it happens divided by the probability that it does not happen:

odds=\frac{p}{1-p}

A probability of 0.8 means odds of 0.8/0.2=4 (“4 to 1”). A probability of 0.5 means odds of 1. If we solve the sigmoid equation for z, we get:

\ln\left(\frac{p}{1-p}\right)=b_0+b_1X_1+b_2X_2+\dots

Then, the logistic model is a linear model for the log of the odds (the log-odds or logit). This gives us the interpretation of the coefficients:

• The sign of b_1 tells whether X_1 increases (+) or decreases (−) the probability of Y=1.
• The p-value of b_1 tells whether the relationship is statistically significant, as in a linear regression.
• The magnitude: if X_1 increases by \Delta units, the odds are multiplied by e^{b_1\Delta}. This is called the odds ratio. Be careful: the odds are multiplied, not the probability.

# Example: how probabilities and odds change
for p in [0.2, 0.5, 0.6, 0.8]:
    print(f"p = {p:.1f} -> odds = {p / (1 - p):.2f} -> log-odds = {np.log(p / (1 - p)):+.2f}")

p = 0.2 -> odds = 0.25 -> log-odds = -1.39
p = 0.5 -> odds = 1.00 -> log-odds = +0.00
p = 0.6 -> odds = 1.50 -> log-odds = +0.41
p = 0.8 -> odds = 4.00 -> log-odds = +1.39

4 Preparing the data

We will predict whether the stock return of a firm will beat the median return of its industry in the next year, using two accounting signals of the current year:

• oepspw: operating earnings per share deflated by price (winsorized), the “earnings yield”. Investors call stocks with a high earnings yield value stocks.
• bmrw: book-to-market ratio (winsorized), another classic value signal.

The question behind the model is a famous one in Finance: do cheap stocks (relative to their earnings and book value) earn higher returns than expensive stocks?

Why do we compare with the industry median and not with the market? Two reasons: (1) the median splits the firms of each industry-year in two halves, so our target will have approximately 50% of ones and 50% of zeros (a balanced target), and (2) the comparison removes the effect of the market and of the industry in each year: we are only trying to predict which firms will do better than their peers.

We start with the same data management of Workshop 3. The following code puts the steps together in one function. Read it carefully: you wrote every line in Workshop 3.

def import_csv(file, base_url="https://www.apradie.com/datos/"):
    """Downloads a file from base_url and import it into a data frame ."""
    import os
    import io
    import requests
    import pandas

    headers = {"User-Agent": "Mozilla/5.0 (Windows NT 10.0; Win64; x64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/91.0.4472.124 Safari/537.36"}
    try: 
        response = requests.get(base_url + file, headers= headers)
        response.raise_for_status() # Error if the response is 4xx or 5xx
        #Reading the content un a dataframe:
        df = pd.read_csv(io.StringIO(response.text))
        print(f"File {file} was downloaded and imported as a data frame")
    except requests.exceptions.HTTPError as e:
        print(f"HTTP Error occured: {e}")
        raise
    except requests.exceptions.RequestException as e:
        print(f"An error ocurred dufing the request: {e}")
        raise
    return df 

usdata = import_csv("usdata.csv")
usfirms = import_csv("usfirms.csv")

# We do the merge of both dataset as in W3:

df = pd.merge(usdata, usfirms[["empresa", "Nombre", "status", "naics1"]],
                  left_on="firm", right_on="empresa", how="left").drop(columns="empresa")
# We make sure that the length of the original usdata is the same as the df (the merged dataset)
assert len(df) == len(usdata)

File usdata.csv was downloaded and imported as a data frame
File usfirms.csv was downloaded and imported as a data frame

We make sure that the dataset is sorted by firm-year, so all shift calculations will work correctly:

df = df.sort_values(["firm", "year"])

Now we calculate the important financial variables that we will use to compute financial ratios for our model:

df["ebit"] = df["revenue"] - df["cogs"] - df["sgae"] - df["depreciation"]
df["marketvalue"] = df["originalprice"] * df["sharesoutstanding"]
df["bookvalue"] = df["totalassets"] - df["totalliabilities"]
df["stockreturn"] = np.log(df["adjprice"] / df.groupby("firm")["adjprice"].shift(1))

We generate columns (variables) for Operating Earnings per Share Deflated by Price (oepsp) and Book-to-market ratio (bmr) as measures to identify value stocks. Firms with high values of oepsp and bmr are considered as value stocks.

df["oepsp"] = (df["ebit"] / df["sharesoutstanding"]) / df["originalprice"]
df["bmr"] = df["bookvalue"] / df["marketvalue"]

When working with big datasets it is common to have extreme values when calculating ratios due to cases with zero values or close-to-zero values.

When Python calculates a ratio where the denominator is zero, the result is an Infinite value codified as np.inf.

In the cases of the oepsp and bmr ratio it is very common to have zero values for the denominators (sharesoutstanding, originalprice, marketvalue), so it is needed to replace the np.inf values for missing values (np.nan). We do this since the regression models can be estimated when we have np.nan values, but it is not possible to estimate a model with np.inf values!

Then we replace the np.inf values for np.nan values for any column of the dataset:

# infinite values appear when we divide by zero; we treat them as missing
df = df.replace([np.inf, -np.inf], np.nan)

In any multiple regression model and logistic regression model when any independent variable (X variable) has very extreme values it is likely that the regression model will estimate non-reliable coefficients and/or standard error of coefficients. Then, these extreme values can be influential for the final model.

The most practical measure to deal with these extreme values is to winsorize these variables. What is Winsorization?

Winsorization is a technique to treat extreme values (outliers): the values below a low percentile are replaced by the value of that percentile, and the values above a high percentile are replaced by the value of that percentile. Unlike deleting outliers, winsorization keeps all the observations.

Then Winsorization “flattens” the extreme values of a variable. For the high values, we choose a percentile (for example, the 99th), and all the values above it are replaced by the value of that percentile; the same for the low values with a low percentile (for example, the 1st). The observations are not deleted: they keep their position (very high or very low), but they lose their extreme magnitude. It is important before a regression because the estimated coefficients (ordinary least squares, logistic regression) are very sensitive to extreme values: a few absurd observations, like an oepsp of 10,000 caused by a price of $0.001, can completely change the estimated beta coefficients and make them unreliable.

To quickly see whether a variable has very extreme values se can do a histogram of the variable.

Bellow are the histograms for the book-to-market and the operating epsp:

import matplotlib.pyplot as plt
import seaborn as sns

plt.figure(figsize=(10, 6))
sns.histplot(df['oepsp'].dropna(), bins=50, kde=True)
plt.title('Distribution of Operational EPS deflated by Stock Price (oepsp)')
plt.xlabel('oepsp')
plt.ylabel('Frequency')
plt.show()

We can see that there are very, very few extreme values for oepsp; the maximum value is more than 10 million! You have to see the 1e7 in the X axix indicating the scale of oepsp; 1e7 is scientific notation for 10 million (1 x 10^7)

import matplotlib.pyplot as plt
import seaborn as sns

plt.figure(figsize=(10, 6))
sns.histplot(df['bmr'].dropna(), bins=50, kde=True)
plt.title('Distribution of Book-to-market ratio (bmr)')
plt.xlabel('bmr')
plt.ylabel('Frequency')
plt.show()

We see that the same proble as in the oepsp ratio, there are very extreme values for bmr.

We must winsorize both ratios before including these ratios as X variables of any regression model.

We write a function to perform winsorization of any variable in a dataset:

def winsorize(x, low=0.01, high=0.99):
    """Replaces the values of x below its `low` quantile and above its `high` quantile
    by those quantiles."""
    # The default percentiles for winsorization are 1 and 99%
    lower = x.quantile(low)
    upper = x.quantile(high)
    return x.clip(lower=lower, upper=upper)
# the clip function replaces high or low values with specific percentiles

We do the winsorization of oepsp and bmr by applying the function winsorize using the transform method of Pandas.

We can apply the winsorization by year since we have extreme values for all years and it is better to do the process by groups for each year.

df["oepspw"] = df.groupby("year")["oepsp"].transform(winsorize)
df["bmrw"] = df.groupby("year")["bmr"].transform(winsorize)

We now plot the winsorized versions for both ratios to see whether extreme values where flattened:

import matplotlib.pyplot as plt
import seaborn as sns

plt.figure(figsize=(10, 6))
sns.histplot(df['oepspw'].dropna(), bins=50, kde=True)
plt.title('Distribution of WINSORIZED Operational EPS deflated by Stock Price (oepsp)')
plt.xlabel('oepspw')
plt.ylabel('Frequency')
plt.show()

import matplotlib.pyplot as plt
import seaborn as sns

plt.figure(figsize=(10, 6))
sns.histplot(df['bmrw'].dropna(), bins=50, kde=True)
plt.title('Distribution of WINSORIZED Book-to-market ratio (bmr)')
plt.xlabel('bmrW')
plt.ylabel('Frequency')
plt.show()

Now there are not so extreme values of the ratios; we could winsorize for a higher percentile to reduce some very negative values for oepsp, but for now we go with this level of winsorization.

NoteSurvivorship bias

We keep all firms, active and cancelled. If we kept only the firms that are active today, we would drop the firms that went bankrupt or were delisted, which are usually the ones with the worst returns. Our sample would look better than reality, and our model would learn from a biased sample. This is called survivorship bias, and it is one of the most common errors in investment research.

5 CHALLENGE 1: Building the target variable

1. Calculate the median stockreturn of each industry (naics1) in each year, and assign it to each firm-year (column ind_median_ret). Hint: transform("median"), as in Workshop 3.
2. Create beat_ind: 1 if the firm’s stockreturn is higher than ind_median_ret, 0 otherwise, and missing if the return is missing.
3. Create the target beat_ind_next: the value of beat_ind of the next year of the same firm. Hint: a negative shift, shift(-1), within each firm.
4. Create the modeling dataset data with the columns firm, year, naics1, oepspw, bmrw, beat_ind_next and stockreturn_next (the next-year return of the firm, which we will use in Challenge 5), dropping the rows with any missing value.
5. How many observations does data have? What percentage of them has beat_ind_next = 1?
6. Explain with your own words why we must use the next year’s return as the target, and not the return of the same year of the accounting data.

5.1 SOLUTION: Challenge 1

df["ind_median_ret"] = df.groupby(["naics1", "year"])["stockreturn"].transform("median")
df["beat_ind"] = (df["stockreturn"] > df["ind_median_ret"]).astype(float)
df.loc[df["stockreturn"].isna(), "beat_ind"] = np.nan

df["beat_ind_next"] = df.groupby("firm")["beat_ind"].shift(-1)
df["stockreturn_next"] = df.groupby("firm")["stockreturn"].shift(-1)

cols = ["firm", "year", "naics1", "oepspw", "bmrw", "beat_ind_next", "stockreturn_next"]
data = df[cols].dropna().copy()
print(f"Observations: {len(data):,}")
print(f"Percentage of ones: {data['beat_ind_next'].mean():.2%}")

Observations: 53,537
Percentage of ones: 49.90%

Explanation.

• transform("median") gives each firm-year the median return of its industry in that year, as in Workshop 3.
• astype(float) converts True/False into 1.0/0.0. We use float instead of int because an integer column cannot contain NaN, and we need NaN when the return is missing (otherwise, a missing return would be counted as “did not beat the industry”).
• shift(-1) moves the column one row up, within each firm: each row gets the value of the next year of the same firm. The last year of each firm gets NaN.
• The percentage of ones is close to (but a little below) 50%: by definition of the median, about half of the firms beat it, and firms whose return is exactly equal to the median are counted as zeros.

Why the next year? The financial statements of year t are published a few weeks or months after the end of year t. When an investor knows the earnings of year t, the stock return of year t has already happened. If we used the same-year return, the model would use information that was not available at the time of the decision (this is called look-ahead bias), and its performance would be an illusion. With the next-year return, the model answers a real investment question: with what I know today, which stocks should I buy for next year? (To be even more rigorous, we could measure returns from April to April, after the financial statements are published; we keep calendar years for simplicity.)

6 CHALLENGE 2: The logistic regression model

Using all the observations of data:

1. Run a logistic regression of beat_ind_next on oepspw and bmrw with statsmodels: smf.logit("beat_ind_next ~ oepspw + bmrw", data=data).fit(), and show the summary.
2. Interpret the sign and the significance of each coefficient with your own words.
3. Interpret the magnitude of the coefficient of oepspw with the odds ratio for an increase of +0.10 in oepspw (e^{b_1\times0.10}).
4. Using the model (method predict), calculate the probability of beating the industry for 2 hypothetical firms with bmrw = 0.5: firm A with oepspw = -0.20 (operating losses) and firm B with oepspw = +0.10.

6.1 SOLUTION: Challenge 2

logit_all = smf.logit("beat_ind_next ~ oepspw + bmrw", data=data).fit()
print(logit_all.summary())

Optimization terminated successfully.
         Current function value: 0.690904
         Iterations 5
                           Logit Regression Results                           
==============================================================================
Dep. Variable:          beat_ind_next   No. Observations:                53537
Model:                          Logit   Df Residuals:                    53534
Method:                           MLE   Df Model:                            2
Date:                Sat, 17 Oct 2026   Pseudo R-squ.:                0.003234
Time:                        11:34:23   Log-Likelihood:                -36989.
converged:                       True   LL-Null:                       -37109.
Covariance Type:            nonrobust   LLR p-value:                 7.639e-53
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept     -0.0675      0.012     -5.448      0.000      -0.092      -0.043
oepspw         0.3605      0.027     13.484      0.000       0.308       0.413
bmrw           0.1251      0.016      7.899      0.000       0.094       0.156
==============================================================================

How to read the output. The table has the same structure as the output of any regression model: coefficient, standard error, z-statistic (it plays the role of the t-statistic), p-value and confidence interval. The first lines show that the model was estimated by maximum likelihood with an iterative algorithm (“Optimization terminated successfully… Iterations”): the computer searched, step by step, for the coefficients that make the observed data most likely. Unlike linear regression, there is no closed formula; we will program a similar search algorithm in Challenge 6. The Pseudo R-squared is a measure of fit for logistic models; it is usually very small in return prediction models.

b = logit_all.params
p = logit_all.pvalues
print(b.round(4))
print(p.round(4))

Intercept   -0.0675
oepspw       0.3605
bmrw         0.1251
dtype: float64
Intercept    0.0
oepspw       0.0
bmrw         0.0
dtype: float64

Sign and significance. The coefficient of oepspw is 0.3605 with a p-value of 0.0000, and the coefficient of bmrw is 0.1251 with a p-value of 0.0000. A positive and significant coefficient (p-value < 0.05) means that firms with a higher value of that variable have a higher probability of beating their industry in the next year. A negative and significant coefficient means the opposite. If the p-value is higher than 0.05, we cannot conclude that the variable is related to the probability.

odds_ratio = np.exp(b["oepspw"] * 0.10)
print(f"Odds ratio for +0.10 in oepspw: {odds_ratio:.4f}")

Odds ratio for +0.10 in oepspw: 1.0367

Magnitude. A change of +1.00 in oepspw is unrealistic (most values are between −0.3 and +0.3), so we interpret a change of +0.10. If a firm has an oepspw 0.10 higher than another firm (with the same bmrw), its odds of beating its industry next year are multiplied by 1.0367, that is, they change by +3.67%. Remember: the odds are multiplied, not the probability. To talk about probabilities, we calculate them directly:

firms = pd.DataFrame({"oepspw": [-0.20, 0.10], "bmrw": [0.5, 0.5]}, index=["Firm A", "Firm B"])
firms["prob_beat"] = logit_all.predict(firms)
firms

	oepspw	bmrw	prob_beat
Firm A	-0.2	0.5	0.480741
Firm B	0.1	0.5	0.507765

predict applies the equation p=\frac{1}{1+e^{-(b_0+b_1X_1+b_2X_2)}} to each row. We can check it by hand for firm B:

z_b = b["Intercept"] + b["oepspw"] * 0.10 + b["bmrw"] * 0.5
print(f"Probability of firm B by hand: {1 / (1 + np.exp(-z_b)):.4f}")

Probability of firm B by hand: 0.5078

The difference between the two probabilities is the effect, in probability terms, of moving from operating losses to a healthy earnings yield. Notice that the probabilities are close to 50%: accounting signals give only a small edge in predicting returns. This is normal: if it were easy to predict stock returns, everybody would do it, and the opportunity would disappear.

7 The Machine Learning workflow

The model of Challenge 2 was estimated with all the data, so it tells us about relationships in the past, but not how well it would predict new data. The ML workflow is:

1. Split the data into a training sample and a testing sample.
2. Train (estimate) the model using only the training sample.
3. Predict the probabilities for the testing sample.
4. Classify: if the probability is higher than a threshold (usually 0.5), predict 1; otherwise, predict 0.
5. Evaluate the predictions against the real values of the testing sample.

7.1 Why we split by time, not at random

Many ML tutorials split the data at random (for example, 80% and 20%). With financial panel data this is a mistake: a random split puts observations of 2020 in the training sample and observations of 2015 in the testing sample, so the model “learns from the future” to predict the past. Also, the returns of all firms in the same year are affected by the same events (a crisis, a rally), so a random split leaks information about each year into the training sample.

The correct approach mimics real life: train with the past, test with the future. We will use the years up to 2016 for training, and the years after 2016 for testing.

7.2 Measures of performance: the confusion matrix

The confusion matrix compares the predictions with the real values in the testing sample:

	Predicted 0	Predicted 1
Actual 0	True Negatives (TN)	False Positives (FP)
Actual 1	False Negatives (FN)	True Positives (TP)

From it we calculate:

• Accuracy =\frac{TP+TN}{\text{all}}: the proportion of correct predictions.
• Precision =\frac{TP}{TP+FP}: of the stocks that the model predicted as winners, what proportion were really winners? For an investor who buys the stocks predicted as 1, this is the hit rate of the buy signal.
• Sensitivity (recall, true positive rate) =\frac{TP}{TP+FN}: of the real winners, what proportion did the model identify?
• Specificity (true negative rate) =\frac{TN}{TN+FP}: of the real losers, what proportion did the model identify correctly?

Always compare the accuracy with a naive benchmark: a “model” that always predicts the most frequent class. If 52% of the test observations are zeros, predicting always 0 gives an accuracy of 52% without any model at all.

8 CHALLENGE 3: Training, testing and the confusion matrix

1. Split data into train (years up to 2016) and test (years after 2016). How many observations does each sample have?
2. Train the logistic model with the training sample only, using scikit-learn, the most popular ML library in Python:

from sklearn.linear_model import LogisticRegression
features = ["oepspw", "bmrw"]
model = LogisticRegression(C=1e10, max_iter=1000)
model.fit(train[features], train["beat_ind_next"])

The parameter C controls regularization, a technique that shrinks the coefficients to prevent overfitting. scikit-learn applies regularization by default; with a very large C we practically turn it off, so that the coefficients are the same as those of statsmodels. Check that they are.

3. Predict the probabilities of the testing sample (model.predict_proba(test[features])[:, 1]) and classify them with a threshold of 0.5.
4. Show the confusion matrix and calculate accuracy, precision, sensitivity and specificity. Interpret each one in terms of an investor.
5. Compare the accuracy with the naive benchmark. Is the model better?

8.1 SOLUTION: Challenge 3

SPLIT_YEAR = 2016
train = data[data["year"] <= SPLIT_YEAR]
test = data[data["year"] > SPLIT_YEAR]
print(f"Training: {len(train):,} observations ({train['year'].min()}-{train['year'].max()})")
print(f"Testing:  {len(test):,} observations ({test['year'].min()}-{test['year'].max()})")

Training: 33,837 observations (2000-2016)
Testing:  19,700 observations (2017-2022)

Remember that the features are from year t and the target is from t+1. The last training year is 2016, so its target uses returns of 2017. The first testing year is 2017 (target in 2018). The training and testing targets do not overlap.

from sklearn.linear_model import LogisticRegression

features = ["oepspw", "bmrw"]
model = LogisticRegression(C=1e10, max_iter=1000)
model.fit(train[features], train["beat_ind_next"])

# Same model with statsmodels, to compare the coefficients
logit_train = smf.logit("beat_ind_next ~ oepspw + bmrw", data=train).fit(disp=0)
print("scikit-learn:", model.intercept_.round(4), model.coef_.round(4))
print("statsmodels: ", logit_train.params.round(4).values)

scikit-learn: [-0.1145] [[0.1938 0.1927]]
statsmodels:  [-0.1144  0.1939  0.1927]

Both libraries give the same coefficients (up to small numerical differences). statsmodels is designed for statistical inference (standard errors, p-values), while scikit-learn is designed for prediction (it has all the tools of the ML workflow, with the same syntax for dozens of models). disp=0 hides the optimization messages.

The syntax of scikit-learn is the same for every model: create the model, .fit(X, y), and .predict_proba(X). The features must be a data frame (or a 2D array), which is why we use train[features] with a list of names.

prob_test = model.predict_proba(test[features])[:, 1]
pred_test = (prob_test > 0.5).astype(int)
y_test = test["beat_ind_next"].astype(int)
print(prob_test[:5].round(4), pred_test[:5])

[0.484  0.4837 0.4818 0.479  0.4781] [0 0 0 0 0]

predict_proba returns a matrix with 2 columns: the probability of 0 and the probability of 1. [:, 1] selects all the rows of the second column (the probability of 1).

from sklearn.metrics import confusion_matrix, ConfusionMatrixDisplay

cm = confusion_matrix(y_test, pred_test)
ConfusionMatrixDisplay(cm, display_labels=["0: does not beat", "1: beats"]).plot(cmap="Blues")
plt.title("Confusion matrix (testing sample)")
plt.show()

tn, fp, fn, tp = cm.ravel()
accuracy = (tp + tn) / cm.sum()
precision = tp / (tp + fp)
sensitivity = tp / (tp + fn)
specificity = tn / (tn + fp)
naive = max(y_test.mean(), 1 - y_test.mean())

print(f"TN={tn}, FP={fp}, FN={fn}, TP={tp}")
print(f"Accuracy:    {accuracy:.4f}")
print(f"Precision:   {precision:.4f}")
print(f"Sensitivity: {sensitivity:.4f}")
print(f"Specificity: {specificity:.4f}")
print(f"Naive benchmark accuracy: {naive:.4f}")

TN=6183, FP=3671, FN=5991, TP=3855
Accuracy:    0.5095
Precision:   0.5122
Sensitivity: 0.3915
Specificity: 0.6275
Naive benchmark accuracy: 0.5002

Explanation. cm.ravel() “flattens” the 2×2 matrix into 4 numbers, reading from top-left to bottom-right: TN, FP, FN, TP. We calculate each measure with its definition (scikit-learn also has functions like precision_score, but calculating them by hand ensures that you understand them).

Interpretation (in terms of an investor who buys the stocks predicted as 1):

• Accuracy of 50.95%: the model classifies correctly this proportion of the firm-years of the testing sample.
• Precision of 51.22%: of the stocks that the model recommended buying, this proportion actually beat their industry the next year. It is the hit rate of our buy signal.
• Sensitivity of 39.15%: of all the stocks that actually beat their industry, the model identified this proportion.
• Specificity of 62.75%: of all the stocks that did not beat their industry, the model correctly avoided this proportion.

Is the model better than the naive benchmark? Always predicting the most frequent class gives an accuracy of 50.02%. Our model’s accuracy is 50.95%, a difference of +0.93 percentage points. In stock return prediction, even a small edge can be economically important when applied to thousands of stocks, but we should not be impressed by accuracies close to 50%. In Challenge 5 we will evaluate the model with the metric that matters to an investor: returns.

9 CHALLENGE 4: The ROC curve and the choice of the threshold

The threshold of 0.5 is arbitrary. With a higher threshold, the model predicts “1” only when it is very confident: fewer buy signals, probably more precise. With a lower threshold, it predicts “1” more often. The ROC curve (Receiver Operating Characteristic) shows the sensitivity (true positive rate) against 1 − specificity (false positive rate) for all possible thresholds. The area under the ROC curve (AUC) summarizes the model’s ability to rank the observations: an AUC of 0.5 is a random guess, and 1.0 is a perfect model.

1. Plot the ROC curve of the testing sample and calculate the AUC (from sklearn.metrics import roc_curve, roc_auc_score).
2. Using the training sample, calculate precision, sensitivity and specificity for thresholds from 0.40 to 0.60 (step 0.01), and plot them. Explain the trade-off between sensitivity and specificity.
3. Why must we choose the threshold with the training sample, and not with the testing sample?

9.1 SOLUTION: Challenge 4

from sklearn.metrics import roc_curve, roc_auc_score

fpr, tpr, thresholds = roc_curve(y_test, prob_test)
auc = roc_auc_score(y_test, prob_test)

plt.figure(figsize=(6, 6))
plt.plot(fpr, tpr, label=f"Logistic model (AUC = {auc:.4f})")
plt.plot([0, 1], [0, 1], linestyle="--", color="grey", label="Random guess (AUC = 0.5)")
plt.xlabel("False positive rate (1 - specificity)")
plt.ylabel("True positive rate (sensitivity)")
plt.title("ROC curve (testing sample)")
plt.legend(); plt.grid(alpha=0.3); plt.show()
print(f"AUC = {auc:.4f}")

AUC = 0.5317

Interpretation. The AUC can be interpreted as the probability that the model assigns a higher probability to a randomly chosen winner than to a randomly chosen loser. An AUC of 0.5317 means that the model ranks the stocks somewhat better than a random guess (0.5), but far from perfectly. The ROC curve is only slightly above the diagonal. This is typical of return prediction models.

prob_train = model.predict_proba(train[features])[:, 1]
y_train = train["beat_ind_next"].astype(int)

rows = []
for th in np.arange(0.40, 0.605, 0.01):
    pred = (prob_train > th).astype(int)
    tn_, fp_, fn_, tp_ = confusion_matrix(y_train, pred, labels=[0, 1]).ravel()
    rows.append({"threshold": round(th, 2),
                 "precision": tp_ / (tp_ + fp_) if (tp_ + fp_) > 0 else np.nan,
                 "sensitivity": tp_ / (tp_ + fn_),
                 "specificity": tn_ / (tn_ + fp_)})
th_table = pd.DataFrame(rows)

th_table.plot(x="threshold", y=["precision", "sensitivity", "specificity"], marker="o", figsize=(9, 5))
plt.title("Performance measures vs threshold (training sample)")
plt.grid(alpha=0.3); plt.show()
th_table.round(4).head()

	threshold	precision	sensitivity	specificity
0	0.40	0.4988	0.9963	0.0049
1	0.41	0.4989	0.9955	0.0064
2	0.42	0.4989	0.9938	0.0080
3	0.43	0.4990	0.9922	0.0098
4	0.44	0.4992	0.9897	0.0131

Explanation. The loop evaluates each threshold, and stores the results in a list of dictionaries (our usual pattern). labels=[0, 1] forces a 2×2 confusion matrix even if, for an extreme threshold, the model predicts only one class. The condition if (tp_ + fp_) > 0 else np.nan avoids a division by zero when the model predicts no ones at all.

The trade-off. When the threshold goes up, the model predicts “1” less often: sensitivity goes down (we miss more winners) and specificity goes up (we avoid more losers). They move in opposite directions: we cannot improve both by changing the threshold; only a better model (better features) can improve both. Precision usually increases with the threshold, because the model only says “buy” when it is more confident, but it gives fewer buy signals. The right threshold depends on the cost of each type of error: a bank that grants loans cares much more about false positives (lending to a client who defaults) than about false negatives.

Why the training sample? The testing sample must simulate new data that the model has never seen. If we choose the threshold that gives the best results in the testing sample, we are using the test data to make modeling decisions, and the test results will be too optimistic (this is a form of data snooping or overfitting). All decisions (features, threshold, parameters) are made with the training data; the testing sample is used only once, at the end, to evaluate the final model.

10 CHALLENGE 5: Is the model useful for an investor?

The best way to evaluate a return-prediction model is to invest with it. In the testing sample:

1. For each year, classify the firms into 5 groups (quintiles) according to their predicted probability: group 1 with the lowest probabilities and group 5 with the highest. Hint: test.groupby("year")["prob"].transform(lambda x: pd.qcut(x, 5, labels=False)) + 1.
2. Calculate the average next-year return (stockreturn_next) of each quintile, averaging first within each year and then across years.
3. Calculate the difference between quintile 5 and quintile 1 (a long-short strategy: buy the stocks with the highest probability, and short-sell the ones with the lowest). Is the difference positive? Plot the average return by quintile.

10.1 SOLUTION: Challenge 5

test = test.copy()
test["prob"] = prob_test
test["quintile"] = test.groupby("year")["prob"].transform(lambda x: pd.qcut(x, 5, labels=False)) + 1

by_year_q = test.groupby(["year", "quintile"])["stockreturn_next"].mean().reset_index()
q_returns = by_year_q.groupby("quintile")["stockreturn_next"].mean()
print(q_returns.round(4))
long_short = q_returns.loc[5] - q_returns.loc[1]
print(f"Average annual return of the long-short strategy (Q5 - Q1): {long_short:.4f}")

q_returns.plot(kind="bar", figsize=(7, 4), title="Average next-year cc return by quintile of predicted probability")
plt.xlabel("Quintile (1 = lowest probability, 5 = highest)"); plt.ylabel("Average return")
plt.grid(alpha=0.3, axis="y"); plt.show()

quintile
1   -0.2832
2   -0.0813
3   -0.0450
4   -0.0223
5    0.0016
Name: stockreturn_next, dtype: float64
Average annual return of the long-short strategy (Q5 - Q1): 0.2848

Explanation.

• test.copy(): test was created by filtering data, so we make an independent copy before adding columns.
• A lambda is a small anonymous function written in one line: lambda x: pd.qcut(x, 5, labels=False) is equivalent to def f(x): return pd.qcut(x, 5, labels=False). pd.qcut divides the values into 5 groups with the same number of observations, labeled 0 to 4 (we add 1 to get 1 to 5). With transform, this is done within each year, so we compare firms only with the firms of the same year.
• We average first within each year-quintile, and then across years, so that each year has the same weight (years with more firms do not dominate the result).

Interpretation. If the model is useful, the average returns should increase from quintile 1 to quintile 5, and the long-short return should be positive. The difference between Q5 and Q1 is 0.2848 per year in continuously compounded terms. The average returns increase monotonically from quintile 1 to quintile 5: the model ranks stocks in a useful way, even though its accuracy is only slightly above 50%. The ranking matters more than the accuracy.

You may be surprised that the average cc returns of most quintiles are negative. This is not an error. Our sample includes thousands of small and very volatile stocks, and cc returns penalize volatility (remember the volatility drag of Workshop 2): a stock that goes from $1 to $0.10 has a cc return of −230%. Quintile 1 is full of firms with operating losses and high distress, which are the ones that suffer these huge losses.

Sorting stocks into portfolios by a signal is how quantitative hedge funds evaluate their signals: not with accuracy, but with the returns of portfolios sorted by the signal. You will study these strategies (long-short, market-neutral, factor investing) in the Hedge Funds course. Keep in mind that we have not considered transaction costs, the difficulty of short-selling small stocks, or the risk of the strategy.

11 DataCamp online course

You MUST TAKE Chapter 1: Classification from the course Supervised Learning with scikit-learn.

12 W4 submission

1. Submit the link of your Google Colab notebook in Canvas, and make sure you shared it with cdorante@tec.mx (Edit privileges).
2. Answer the Canvas Quiz W4 before the deadline, with your notebook open.