Workshop 4, 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. In this case, you do not have to submit your Notebook, but it is strongly recommended to do the complete workshop yourself in and out of class.

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, 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.

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.

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?

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?

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.

11 CHALLENGE 6 (Algorithms): Programming the logistic regression with gradient descent

How does the computer find the coefficients of the logistic model? There is no formula: it uses an optimization algorithm. The model looks for the coefficients that minimize the log-loss (the average error of the predicted probabilities):

LogLoss=-\frac{1}{N}\sum_{i=1}^{N}\left[y_i\ln(p_i)+(1-y_i)\ln(1-p_i)\right]

Gradient descent is the simplest optimization algorithm: start with any coefficients (for example, zeros), and repeat many times: calculate the slope (gradient) of the log-loss, and move the coefficients a small step in the opposite direction (downhill). For the logistic model, the gradient has a very simple formula:

gradient=\frac{1}{N}X^{T}(p-y)

where X is the matrix of features with a column of ones for the intercept, p the vector of predicted probabilities, and y the vector of real values.

SET b = vector of zeros (one per column of X)
REPEAT 5,000 times:
    p = sigmoid(X @ b)
    gradient = X.T @ (p − y) / N
    b = b − learning_rate × gradient
RETURN b

1. Using the training sample, standardize the 2 features (subtract their mean and divide by their standard deviation). Gradient descent works much better when all the features have a similar scale.
2. Build X as a NumPy array with a column of ones and the 2 standardized features, and y as an array of the target.
3. Program the algorithm with a learning rate of 0.5 and 5,000 iterations. Store the log-loss every 500 iterations to check that it goes down.
4. Run the same logistic model with statsmodels using the standardized features, and compare the coefficients with yours.
5. Explain with your own words: what is the connection between this algorithm and the bisection algorithm of Workshop 1?

12 DataCamp online course

It is RECOMMENDED TO TAKE Chapter 1: Classification and Chapter 3: Fine-Tunning your Model from the course Supervised Learning with scikit-learn.

13 W4 submission

In this case, you do not have to submit your Notebook. This workshop will help you to be prepared for the Final Exam.