Building your first statistical model

A note on reading time: The reading time shown above is too high. It also counts the data and code behind the interactive charts. The article itself takes about seven minutes to read.

Picture a buyer at a retail chain in September, deciding how much stock to order for December. In the U.S. retail data used in this article, December alone brought in 9.7 percent of all sales in 2015, nearly a tenth of the year in one month. Order too little and the shelves run empty. With too much, the rest goes out in January at a discount.

The buyer needs a number for a month that has not happened yet, and the only material at hand is the past. Here is that past: 293 months of U.S. retail sales, from January 1992 to May 2016.

The past repeats itself, though never exactly

Most of what a forecaster needs is visible in the chart. Sales climb over time, from about 146,000 in January 1992 to about 471,000 in May 2016, in the file's own units. Each year also has the same shape: in 23 of the 24 full years, December was the strongest month, and in 23 of them the weakest month was January or February. Then, in 2008 and 2009, the climb stops, and annual sales in the file fall 1.1 and then 7.4 percent.

A series like this, one value recorded at regular intervals and kept in time order, is called a time series. Order matters because this month's sales resemble last month's. The same structure sits behind electricity demand, website traffic and hospital admissions, so the method in this article carries over to them.

The climb and the yearly shape are regular, which means a method that learns them from the past can extend them forward. A recession is a different kind of thing. Nothing in the years before 2008 announced it.

Forecasting methods range from rules of thumb to neural networks

The simplest forecast needs no model at all. Take last December's sales and call it the forecast for this December. Forecasters call this the seasonal naive rule, and it is the baseline any serious method has to beat.

One step up, smoothing methods average recent months and give the newest ones the most weight, so the forecast adapts as the series moves. Statistical models such as ARIMA go further and describe how each month depends on the months before it. Machine-learning methods, from decision trees to neural networks, can take in many more inputs, such as prices, holidays or weather, and find patterns nobody wrote down in advance.

Each step buys flexibility and costs transparency. A neural network can fit the past closely and still give no reason for its forecast. For a first model, the opposite trade is the more useful one.

A linear regression turns the calendar into a recipe

This article uses linear regression, one of the oldest tools in statistics. It draws the straight line that fits a cloud of points best, and the result reads like a recipe. Here the recipe has two ingredients: a fixed amount of growth per month, and a fixed adjustment for each calendar month, large and positive for December and negative for February. Every number the model learns can be printed and checked.

The recipe has an obvious limit. It can only extend patterns it has already seen, so a straight line cannot see a recession coming.

Under twenty lines of Python build and test the model

For readers who code, the full script follows. Pandas and scikit-learn, two free Python libraries for tables and machine learning, are the only extras. If you do not code, skip the box, because the text describes each step.

import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression

# 1. Load a real series: monthly U.S. retail sales, January 1992 to May 2016
url = "https://raw.githubusercontent.com/facebook/prophet/main/examples/example_retail_sales.csv"
df = pd.read_csv(url, parse_dates=["ds"], index_col="ds")
y = df["y"]

# 2. Turn the calendar into columns: a month counter and the month of the year
df["t"] = np.arange(len(df))
df["month"] = df.index.month
X = pd.get_dummies(df[["t", "month"]], columns=["month"], dtype=float)

# 3. Split by date, never at random: hold out the last 60 months
cut = len(df) - 60
X_train, X_test = X.iloc[:cut], X.iloc[cut:]
y_train, y_test = y.iloc[:cut], y.iloc[cut:]

# 4. Fit the line and predict the unseen months
model = LinearRegression().fit(X_train, y_train)
pred = model.predict(X_test)

# 5. Compare with the naive guess: "same month one year earlier"
naive = y.shift(12).iloc[cut:]
print(f"model MAE: {np.abs(y_test - pred).mean():.0f}")
print(f"naive MAE: {np.abs(y_test - naive).mean():.0f}")

Step 1 loads the data, one number per month. A model cannot read a calendar, so step 2 puts the data into a form it can read. It adds a counter that numbers the months from 0 to 292 and carries the climb.

The month of the year needs different handling. Written as a plain number from 1 to 12, it would force the model to raise its adjustment in equal steps from January to December, as if December were simply more month than November. So each calendar month gets its own column instead, set to 1 when a row belongs to that month and 0 otherwise. The columns work like light switches. In a December row only the December switch is on, so the model adds December's adjustment and no other. The table below shows what the data looks like once step 2 is done.

Step 4 fits the line, and the learned growth comes out at 11,624 per year.

The test comes from the future, and the baseline gets a head start

Step 3 trains on the 233 months from January 1992 to May 2011 and hides the 60 months from June 2011 to May 2016. A random split would show the model the months on both sides of every test month, which turns forecasting into filling gaps, so the split follows the calendar. The score is the mean absolute error, the average size of the miss.

A score means little on its own, so the model needs an opponent. The obvious one is the seasonal naive rule from earlier, which step 5 computes: it forecasts each month with the real value from the same month a year before. Take April 2013, when sales were 408,381. The model predicted 411,600 and missed by 3,219. For the naive rule, the forecast was April 2012’s 393,035, a miss of 15,346.

Repeat that for all 60 test months and average the misses. The model's average miss is 8,411, about 2.0 percent of a typical month's sales, and the naive rule's is 16,466, or 3.9 percent. In 49 of the 60 months the model came closer, and it cut the average error almost in half.

The contest is also tilted against the model. For a month in 2015, the naive rule gets to use real sales from 2014. The model has to forecast all five years without ever seeing 2012, 2013 or 2014. Beating an opponent with fresher information suggests the line has learned something real.

The chart lays the model over the real data. Left of the dotted line, the dashed orange line is the fit to the training years. To the right, the solid orange line is the forecast, drawn without the model ever seeing those months. Switch to the test window to compare the two lines month by month.

The fine points sit in four places this introduction only touches

This article is an introduction, and the model is the simplest one that does the job. A 2 percent miss sounds precise, yet sales grew 4.6 percent a year on average from 2010 to 2015, so the typical miss was about two-fifths of a year's growth. Where that error comes from is where forecasting gets interesting, and each of the following topics deserves a closer look than this piece can give it.

• How the computer finds the β values. The fit solves a least-squares problem, either in one step with linear algebra or in many small steps with gradient descent, the same idea that trains neural networks. Search for ordinary least squares and gradient descent.
• Errors that come in streaks. On the training data, one month's miss correlates with the next month's at 0.86, a sign that the line leaves a pattern behind. Autocorrelation is the term, and ARIMA models are built to use it.
• Breaks such as a recession. Refit every June from 2006 to 2015, and the 12-month forecasts made in 2008 and 2009 missed by 34,650 and 41,516, against 6,796 to 11,035 in calm years. Look up structural breaks and rolling-origin evaluation.
• A seasonal swing that grows. The model learned a December bump 76,709 above February, while the real gap in the test years averaged 111,858. Multiplicative seasonality and feature engineering are the places to start.

The line is a baseline, and a baseline is worth keeping

For a first model, the line does its job. It runs in seconds, came closer than last year's value in 49 of 60 test months, and shows in one chart where the economy moved without it. Whether a 2 percent miss is good enough depends on the decision it feeds, and only the person making that decision can say.

Sources and further reading

• Data: example_retail_sales.csv in the examples folder of the Prophet repository on GitHub. The file states no original source or unit.
• scikit-learn: LinearRegression documentation. The class fits the model by ordinary least squares.
• pandas: get_dummies documentation, the function that builds the twelve month switches.
• How the algorithm works: Rob J Hyndman and George Athanasopoulos, Forecasting: Principles and Practice, 3rd edition, chapter 7 on time series regression. Section 7.2 explains least squares estimation, and section 7.4 covers trend and seasonal dummy variables, the exact setup used here. The book is readable online.
• Visual intuition: 3Blue1Brown, Gradient descent, how neural networks learn (https://www.3blue1brown.com/lessons/gradient-descent). The video shows how a model searches for the numbers that make its squared errors as small as possible, the same idea this model uses, applied to a neural network.