Feature Engineering · Notebook 5 companion articleFeature Engineering for Date and Time
A timestamp like 2012-06-14 17:00:00 holds a lot of information, but a model can't use most of it as written. To a model, it's one number counting up from 1970, and all it can learn from that is "later means more" or "later means less." The useful parts, like it being a Thursday, rush hour or summer, are hidden inside it. You have to pull them out on purpose.
Start hereIntroduction#
Every notebook so far worked with columns that describe something fixed: a contract type, a fare, a cabin. A timestamp is different. It's not one feature. It's many packed together: hour of day, day of week, season, time since some key event. A model given the raw number can't unpack any of that by itself.
Two lessons from notebook 4 come straight back. The right encoding depends on the model, and section 3 shows it again with hours. And notebook 4's warning about values that only appear after deployment has a time version: the future always brings dates the model has never seen.
This notebook shows how to pull those features out, and how the process can quietly go wrong if you don't check when each piece of information was actually known.
Most of the work uses UCI Bike Sharing: 17,379 hourly records of bike rentals in Washington, DC in 2011 and 2012. The target is the number of bikes rented in each hour. It's a good fit because time really drives the result. Demand depends on the hour, day of the week, season, weather and whether people are commuting, so every feature we build can be scored against real numbers.
Two things to know about the file. The weather columns come already scaled to 0 to 1 (for example, temp is the temperature in Celsius divided by 41), which is why they look so small. And two full years would be 17,544 hours, so 165 hours have no row at all. Keep that in mind whenever you assume one row per hour.
A second, smaller e-commerce dataset joins later. Bike rentals have no customer to follow from visit to visit, and section 4 needs one.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.ensemble import RandomForestRegressor
from sklearn.linear_model import Ridge
from sklearn.metrics import r2_score
from sklearn.model_selection import KFold, TimeSeriesSplit, cross_val_score
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
plt.style.use('seaborn-v0_8-whitegrid')
pd.set_option('display.max_columns', 30)
pd.set_option('display.width', 170)
RANDOM_STATE = 42
bikes = pd.read_csv('datasets/bike_sharing.csv', parse_dates=['dteday'])
print(f'Rows : {len(bikes):,}')
print(f'Date range : {bikes["dteday"].min().date()} to {bikes["dteday"].max().date()}')
print(f'Target (cnt): {bikes["cnt"].min()} to {bikes["cnt"].max()} rentals per hour, mean {bikes["cnt"].mean():.0f}')
print()
bikes[['dteday', 'hr', 'workingday', 'holiday', 'temp', 'hum', 'windspeed', 'cnt']].head()
Rows : 17,379 Date range : 2011-01-01 to 2012-12-31 Target (cnt): 1 to 977 rentals per hour, mean 189
dteday hr workingday holiday temp hum windspeed cnt 0 2011-01-01 0 0 0 0.24 0.81 0.0 16 1 2011-01-01 1 0 0 0.22 0.80 0.0 40 2 2011-01-01 2 0 0 0.22 0.80 0.0 32 3 2011-01-01 3 0 0 0.24 0.75 0.0 13 4 2011-01-01 4 0 0 0.24 0.75 0.0 1
Section 1Parsing, and the traps that fail silently#
Parsing is one function call. The real risk is the three ways it goes wrong without any error
Dates arrive as strings from CSVs, APIs and database exports. Until you parse them, pandas treats them as plain text, and you can't pull anything out of text.
unparsed = pd.read_csv('datasets/datetime_data.csv')
print('As loaded from CSV:')
print(unparsed.dtypes.to_string())
print()
print('These look like dates to us, but they are strings to pandas:')
print(unparsed['order_date_str'].head(3).to_list())
print()
try:
unparsed['order_date_str'].dt.month
except AttributeError as err:
print(f'Asking for the month raises: AttributeError: {str(err)[:60]}...')
As loaded from CSV: order_id int64 order_date_str str signup_date_str str delivery_date_str str order_amount float64 These look like dates to us, but they are strings to pandas: ['2024-05-10 08:23:00', '2024-11-10 22:04:00', '2022-05-02 13:57:00'] Asking for the month raises: AttributeError: Can only use .dt accessor with datetimelike values...
orders = unparsed.copy()
for col in ['order_date', 'signup_date', 'delivery_date']:
orders[col] = pd.to_datetime(orders[f'{col}_str'])
orders = orders.drop(columns=[c for c in orders.columns if c.endswith('_str')])
print('After parsing:')
print(orders.dtypes.to_string())
print()
print(f'Now the components are available: month of the first order = {orders["order_date"].dt.month.iloc[0]}')
After parsing: order_id int64 order_amount float64 order_date datetime64[us] signup_date datetime64[us] delivery_date datetime64[us] Now the components are available: month of the first order = 5
Parsing itself is easy. What matters are three ways a date can turn into a believable but wrong value, with no warning. That makes them more dangerous than a date that simply fails to parse.
Trap one: day and month mixed up#
03/04/2023 means 3 April in most of the world and 4 March in the United States. Both are valid dates, so nothing errors. You just get the wrong answer on every ambiguous row. A string with a day after the 12th, like 25/04/2023, can only be read one way, so those rows look fine. That hides the problem on the others if you only check a few rows.
ambiguous = pd.Series(['03/04/2023', '05/06/2023', '11/12/2023'])
us_style = pd.to_datetime(ambiguous, format='%m/%d/%Y')
uk_style = pd.to_datetime(ambiguous, format='%d/%m/%Y')
print(pd.DataFrame({
'raw_string': ambiguous,
'read as MM/DD/YYYY': us_style.dt.strftime('%d %B %Y'),
'read as DD/MM/YYYY': uk_style.dt.strftime('%d %B %Y'),
}).to_string(index=False))
print()
print('Same three strings, six different dates, and no error either way.')
raw_string read as MM/DD/YYYY read as DD/MM/YYYY 03/04/2023 04 March 2023 03 April 2023 05/06/2023 06 May 2023 05 June 2023 11/12/2023 12 November 2023 11 December 2023 Same three strings, six different dates, and no error either way.
Three strings gave six different dates, with no error. The fix is to pass format= whenever the source isn't clear ISO format. It's one extra argument, and it turns silent errors into loud ones: a string that doesn't match the format fails right away instead of guessing.
When you control the format yourself, use ISO 8601 (YYYY-MM-DD). It sorts correctly as text, it can't be misread, and pandas parses it without help.
Trap two: values that aren't dates#
Real columns have blanks, placeholders like N/A or 9999-12-31, and dates that don't exist, like 30 February.
messy = pd.Series(['2023-01-15', 'not a date', '2023-02-30', '2023-03-01'])
print('Input:', messy.to_list())
print()
try:
pd.to_datetime(messy)
except (ValueError, TypeError) as err:
print(f'Default behaviour raises: {type(err).__name__}: {str(err)[:70]}...')
print()
coerced = pd.to_datetime(messy, errors='coerce')
print('With errors="coerce":')
print(coerced.to_string())
print()
print(f'Unparseable values became NaT: {coerced.isna().sum()} of {len(messy)}')
Input: ['2023-01-15', 'not a date', '2023-02-30', '2023-03-01'] Default behaviour raises: ValueError: time data "not a date" doesn't match format "%Y-%m-%d". You might want... With errors="coerce": 0 2023-01-15 1 NaT 2 NaT 3 2023-03-01 Unparseable values became NaT: 2 of 4
Look at 2023-02-30. It's a well-formed string, but February has no 30th, so the calendar catches it, not a format check. errors='coerce' turns every failure into NaT, the date version of NaN. That's usually right for a pipeline, but only if you then count the NaTs. Coercing without checking turns an error into silent data loss. Notebook 1 did this right: after coercing TotalCharges, it counted the 11 blanks it had created and then found out why they were there.
Trap three: time zones#
A timestamp without a time zone is unclear about the one thing people most want from it: the hour.
moment = pd.Timestamp('2023-11-05 06:30:00', tz='UTC')
print('One single instant in time, viewed from three places:')
for zone in ['UTC', 'America/New_York', 'Asia/Kolkata']:
local = moment.tz_convert(zone)
print(f' {zone:20s} {local} hour = {local.hour}')
print()
print('If your feature is "hour of day", the answer depends entirely on the zone.')
print('A model trained on UTC hours learns nothing about human routines.')
One single instant in time, viewed from three places: UTC 2023-11-05 06:30:00+00:00 hour = 6 America/New_York 2023-11-05 01:30:00-05:00 hour = 1 Asia/Kolkata 2023-11-05 12:00:00+05:30 hour = 12 If your feature is "hour of day", the answer depends entirely on the zone. A model trained on UTC hours learns nothing about human routines.
The rule: store times in UTC, and convert to local time only when you extract something people care about, like the hour. Server logs are usually in UTC, and rush hour happens in local time, so convert before pulling out the hour, never after.
The date above was chosen on purpose. On 5 November 2023, US clocks went back an hour, so 01:30 local time happened twice that night. Because of daylight saving, local times aren't always unique, and some don't exist at all (in the US, clocks jump from 1:59 to 3:00 every spring). That's another reason to keep UTC as the source of truth.
Bike Sharing avoids all this, because the hour already comes as a local number in its own column. Whoever prepared the dataset did the time zone work so users wouldn't have to.
Section 2Extracting components, and reading what they show#
The .dt accessor gives you every part of a timestamp. Which parts to keep depends on what you're predicting
Once a column is a real datetime, pandas' .dt accessor gives you every part of it: year, month, day, day of week, quarter and more. Each part answers a different question, and there's no fixed set that's always right to keep.
bikes['date'] = bikes['dteday']
bikes['year'] = bikes['dteday'].dt.year
bikes['month'] = bikes['dteday'].dt.month
bikes['day'] = bikes['dteday'].dt.day
bikes['dayofweek'] = bikes['dteday'].dt.dayofweek # Monday = 0
bikes['dayofyear'] = bikes['dteday'].dt.dayofyear
bikes['quarter'] = bikes['dteday'].dt.quarter
bikes['is_weekend'] = (bikes['dayofweek'] >= 5).astype(int)
bikes['days_since_launch'] = (bikes['dteday'] - bikes['dteday'].min()).dt.days
# 'hr' is already supplied; derive a coarser part-of-day from it
bikes['part_of_day'] = pd.cut(bikes['hr'], bins=[-1, 5, 11, 17, 23],
labels=['night', 'morning', 'afternoon', 'evening'])
print('Derived from one date column plus the hour:')
print(bikes[['dteday', 'hr', 'year', 'month', 'dayofweek', 'is_weekend',
'part_of_day', 'days_since_launch']].head(3).to_string(index=False))
print()
print(f'One timestamp column became 9 usable features, plus the hour already supplied.')
Derived from one date column plus the hour:
dteday hr year month dayofweek is_weekend part_of_day days_since_launch
2011-01-01 0 2011 1 5 1 night 0
2011-01-01 1 2011 1 5 1 night 0
2011-01-01 2 2011 1 5 1 night 0
One timestamp column became 9 usable features, plus the hour already supplied.Extracting the parts is easy. Reading what they show is what shapes the rest of the notebook.
fig, axes = plt.subplots(2, 2, figsize=(14, 8))
by_hour = bikes.groupby('hr')['cnt'].mean()
axes[0, 0].bar(by_hour.index, by_hour.values, color='#3b76af')
axes[0, 0].set(xlabel='hour of day', ylabel='mean rentals', title='Two peaks, not one')
by_month = bikes.groupby('month')['cnt'].mean()
axes[0, 1].bar(by_month.index, by_month.values, color='#c4453c')
axes[0, 1].set(xlabel='month', ylabel='mean rentals', title='Seasonal, and it wraps around')
by_dow = bikes.groupby('dayofweek')['cnt'].mean()
axes[1, 0].bar(range(7), by_dow.values, color='#5c9e5c')
axes[1, 0].set(xticks=range(7), xticklabels=['Mon', 'Tue', 'Wed', 'Thu', 'Fri', 'Sat', 'Sun'],
ylabel='mean rentals', title='Day of week barely matters on its own')
for working, label, colour in [(1, 'working day', '#3b76af'), (0, 'non-working day', '#c4453c')]:
subset = bikes[bikes['workingday'] == working].groupby('hr')['cnt'].mean()
axes[1, 1].plot(subset.index, subset.values, marker='o', ms=3, label=label, color=colour)
axes[1, 1].set(xlabel='hour of day', ylabel='mean rentals',
title='The same hour means different things')
axes[1, 1].legend()
plt.tight_layout()
plt.show()
Three things in these charts shape the rest of this notebook.
There are two daily peaks, not one. Demand rises to 359 rentals at 8am, drops during the middle of the day, then rises to 461 at 5pm. That double peak comes from commuting, and it's why section 3 goes further than most guides.
Day of week looks unimportant on its own, ranging only from about 177 to 196 rentals across the week. Looking at that alone, you might drop it. The bottom-right chart shows why that would be a mistake.
The same hour means almost opposite things on different days. At 8am, a working day averages 477 rentals and a non-working day 106. By noon it flips: 201 against 366. Commuters and leisure riders are two different groups sharing one column, and neither hr nor workingday can show that alone. Section 5 covers what can.
Section 3Cyclical encoding, and where the usual advice stops#
The standard fix is right, but it doesn't go as far as most guides suggest
If you give a model hr as a number from 0 to 23, you're telling it hour 23 is 23 steps away from hour 0. In reality, they're one hour apart. The same problem applies to months, weekdays and anything else that loops back around.
The standard fix puts each value on a circle and uses its coordinates:
xsin = sin(2πt / T), xcos = cos(2πt / T)
Here T is the length of the cycle: 24 for hours, 12 for months. You need both sine and cosine, because sine alone would give two different hours the same value.
def cyclical(values, period, harmonic=1):
angle = 2 * np.pi * harmonic * values / period
return np.sin(angle), np.cos(angle)
pairs = [(23, 0), (23, 22), (11, 12), (0, 12)]
print(f'{"hours":>12s} {"raw distance":>14s} {"cyclical distance":>19s}')
print('-' * 48)
for a, b in pairs:
sa, ca = cyclical(a, 24) # scalars in, scalars out
sb, cb = cyclical(b, 24)
circular = np.hypot(sa - sb, ca - cb)
print(f'{f"{a:02d} vs {b:02d}":>12s} {abs(a - b):>14d} {circular:>19.3f}')
hours raw distance cyclical distance
------------------------------------------------
23 vs 00 23 0.261
23 vs 22 1 0.261
11 vs 12 1 0.261
00 vs 12 12 2.000The encoding works as promised. Hours 23 and 0 are 0.261 apart, the same as 23 and 22, and the same as 11 and 12. Every pair of neighboring hours is equally close. Hours 0 and 12 are 2.000 apart, the maximum, because they're opposite points on a clock.
fig, axes = plt.subplots(1, 2, figsize=(13, 5))
hours = np.arange(24)
axes[0].scatter(hours, np.zeros(24), c=hours, cmap='twilight', s=90)
for h in [0, 6, 12, 18, 23]:
axes[0].annotate(f'{h}', (h, 0), textcoords='offset points', xytext=(0, 12), ha='center')
axes[0].annotate('', xy=(23, 0.04), xytext=(0, 0.04),
arrowprops=dict(arrowstyle='<->', color='red', lw=2))
axes[0].text(11.5, 0.06, 'model thinks these are 23 apart', color='red', ha='center', fontsize=10)
axes[0].set(title='Raw integer hour', xlabel='hr', yticks=[], ylim=(-0.1, 0.12))
sin_h, cos_h = cyclical(hours, 24)
axes[1].scatter(cos_h, sin_h, c=hours, cmap='twilight', s=90)
for h in [0, 6, 12, 18, 23]:
axes[1].annotate(f'{h}', (cos_h[h], sin_h[h]), textcoords='offset points',
xytext=(9, 5), fontsize=9)
axes[1].annotate('', xy=(cos_h[23], sin_h[23]), xytext=(cos_h[0], sin_h[0]),
arrowprops=dict(arrowstyle='<->', color='red', lw=2))
axes[1].text(0.55, 0.35, 'now adjacent', color='red', fontsize=10)
axes[1].set(title='Cyclical encoding', xlabel='cos', ylabel='sin')
axes[1].set_aspect('equal')
plt.tight_layout()
plt.show()
Midnight and 11pm are neighbors on the circle, as they should be. So far this matches the standard advice. What the advice usually skips: does it help a real model on real data, and how does it compare with other options, not just with the raw number?
target = bikes['cnt']
weather = bikes[['temp', 'hum', 'windspeed', 'workingday']]
cv_random = KFold(n_splits=5, shuffle=True, random_state=RANDOM_STATE)
def ridge_r2(X):
return cross_val_score(make_pipeline(StandardScaler(), Ridge()), X, target,
cv=cv_random, scoring='r2').mean()
s1, c1 = cyclical(bikes['hr'], 24, harmonic=1)
s2, c2 = cyclical(bikes['hr'], 24, harmonic=2)
s3, c3 = cyclical(bikes['hr'], 24, harmonic=3)
encodings = {
'no hour feature at all': weather,
'raw integer hour': weather.assign(hr=bikes['hr']),
'cyclical, 1 harmonic': weather.assign(sin1=s1, cos1=c1),
'cyclical, 2 harmonics': weather.assign(sin1=s1, cos1=c1, sin2=s2, cos2=c2),
'cyclical, 3 harmonics': weather.assign(sin1=s1, cos1=c1, sin2=s2, cos2=c2, sin3=s3, cos3=c3),
'one-hot, 24 hours': pd.concat([weather, pd.get_dummies(bikes['hr'], prefix='h').astype(int)], axis=1),
}
print(f'{"encoding":26s} {"columns":>8s} {"Ridge R2":>10s}')
print('-' * 46)
for name, X in encodings.items():
print(f'{name:26s} {X.shape[1]:>8d} {ridge_r2(X):>10.4f}')
encoding columns Ridge R2 ---------------------------------------------- no hour feature at all 4 0.2503 raw integer hour 5 0.3243 cyclical, 1 harmonic 6 0.4236 cyclical, 2 harmonics 8 0.5116 cyclical, 3 harmonics 10 0.5452 one-hot, 24 hours 28 0.6006
This is the key result, and the usual advice would have stopped too early. Cyclical encoding beats the raw number, 0.4236 against 0.3243, as promised. But plain one-hot encoding reaches 0.6006 and beats both. Adding harmonics closes most of the gap: two harmonics reach 0.5116, three reach 0.5452.
The reason shows in the first chart from section 2. One sine and cosine pair can only draw one smooth peak per cycle. That's all it can do. Bike demand has two peaks a day, the morning and evening commutes, and no single wave can make that shape. A second harmonic adds a wave that repeats twice a day, which is what a double peak needs. The next chart measures this.
profile = bikes.groupby('hr')['cnt'].mean()
hours = np.arange(24)
fig, ax = plt.subplots(figsize=(11, 4.5))
ax.plot(hours, profile.values, 'o-', color='black', lw=2, label='actual mean demand')
for n_harmonics, colour in [(1, '#c4453c'), (2, '#3b76af'), (3, '#5c9e5c')]:
design = [np.ones(24)]
for k in range(1, n_harmonics + 1):
s, c = cyclical(hours, 24, harmonic=k)
design += [s, c]
design = np.column_stack(design)
coefficients, *_ = np.linalg.lstsq(design, profile.values, rcond=None)
fitted = design @ coefficients
explained = 1 - ((profile.values - fitted) ** 2).sum() / ((profile.values - profile.mean()) ** 2).sum()
ax.plot(hours, fitted, '--', color=colour, lw=1.8,
label=f'{n_harmonics} harmonic{"s" if n_harmonics > 1 else ""} (explains {explained:.1%})')
ax.set(xlabel='hour of day', ylabel='mean rentals',
title='What each number of harmonics can actually represent')
ax.legend(fontsize=9)
plt.tight_layout()
plt.show()
The red curve is all one sine and cosine pair can do: a single smooth hump that explains 67.4% of the daily pattern. It can't make two peaks, whatever the coefficients. It fits the bigger evening rush and mostly misses the 8am spike, predicting about 149 rentals where the real number is 359. Two harmonics reach 82.2%, and a second bump appears. Three reach 88.9%, and the morning peak lands close to where it should.
Notice that the two- and three-harmonic curves dip a little below zero in the early morning. A sum of waves doesn't know that rentals can't be negative. They're smooth approximations, not models that follow real-world limits.
This applies beyond bikes: use one harmonic when a cycle has one peak, and add more only when it doesn't. Temperature over a year has one peak, so one harmonic is usually enough. Human activity over a day often has two. Other cycles may need more.
For this dataset, one-hot encoding beat every cyclical version. It has 24 free parameters and can fit any daily shape, but it uses 24 columns and never knows that hour 23 and hour 0 are neighbors. Cyclical encoding is compact and smooth. One-hot is flexible. Neither is always the right choice.
One more comparison changes the answer completely.
forest = RandomForestRegressor(n_estimators=80, random_state=RANDOM_STATE, n_jobs=1)
print(f'{"encoding":26s} {"Ridge R2":>10s} {"RandomForest R2":>18s}')
print('-' * 56)
for name in ['raw integer hour', 'cyclical, 1 harmonic', 'one-hot, 24 hours']:
X = encodings[name]
tree_score = cross_val_score(forest, X, target, cv=cv_random, scoring='r2').mean()
print(f'{name:26s} {ridge_r2(X):>10.4f} {tree_score:>18.4f}')
encoding Ridge R2 RandomForest R2 --------------------------------------------------------
raw integer hour 0.3243 0.8058
cyclical, 1 harmonic 0.4236 0.8030
one-hot, 24 hours 0.6006 0.7874
The order flips. For the random forest, the raw hour number is best at 0.8058, cyclical is close behind at 0.8030, and one-hot is worst at 0.7874. A tree never measures distance between values, so it doesn't care that 23 and 0 are far apart as numbers. It just keeps splitting and can cut the hour range into whatever groups help. Spreading that across 24 one-hot columns, or squeezing it into two smooth ones, only makes its job harder.
Notice also how much higher the forest scores overall, about 0.80 against Ridge's 0.60 at best. Part of the reason is that trees find interactions like "hour depends on whether it's a working day" on their own. Section 5 has to build that by hand for Ridge.
This is the same lesson notebook 4 found with label encoding: the right encoding depends on the model, not just the data. Cyclical encoding is a fix for models that measure distance. Calling it "the right way to handle time cycles" without that condition is how people end up making their tree models worse.
Section 4Time deltas, durations and recency#
Often the useful number isn't a date, but the gap between two dates. Pandas makes the safe and the dangerous version equally easy to write
How long has this customer been with us? How long did delivery take? How long since their last purchase? Subtracting one datetime column from another gives a Timedelta, which you convert to a plain number before modeling.
orders['delivery_days'] = (orders['delivery_date'] - orders['order_date']).dt.days
orders['account_age_days'] = (orders['order_date'] - orders['signup_date']).dt.days
print(orders[['order_date', 'signup_date', 'delivery_date',
'account_age_days', 'delivery_days']].head(3).to_string(index=False))
print()
print(orders[['account_age_days', 'delivery_days']].describe().round(1).to_string())
order_date signup_date delivery_date account_age_days delivery_days
2024-05-10 08:23:00 2023-05-19 2024-05-11 357 0
2024-11-10 22:04:00 2023-09-30 2024-11-17 407 6
2022-05-02 13:57:00 2020-07-04 2022-05-05 667 2
account_age_days delivery_days
count 500.0 500.0
mean 351.6 3.9
std 217.1 2.7
min 0.0 0.0
25% 154.5 1.0
50% 345.5 4.0
75% 545.2 6.0
max 729.0 8.0account_age_days ranges from 0 days to 729, just under two years. It's a different feature from the raw signup_date, not just a new format. A model given the signup date learns "people who joined in 2021 behave this way," which stops being useful as 2021 gets further away. A model given account age learns "people two years in behave this way," which keeps working no matter what year it is.
That matters: dates go out of date, but durations don't. A feature based on year has to guess as soon as the model sees a year it wasn't trained on. A duration measured from each customer's own start, like account age, stays in a range the model has seen, however far the calendar moves.
But not every duration is safe. days_since_launch counts from one fixed day, so it keeps growing forever. After the training period, every new row has a larger value than anything the model saw, which is the same problem as year. Section 6 comes back to when that's still useful.
Recency features need the second dataset, because a bike rental has no customer to track.
ecom = pd.read_csv('datasets/ecommerce_transactions.csv', parse_dates=['order_date'])
ecom = ecom.sort_values(['customer_id', 'order_date']).reset_index(drop=True)
# Gap to that customer's PREVIOUS order. shift() inside a groupby keeps customers separate.
ecom['days_since_prev_order'] = ecom.groupby('customer_id')['order_date'].diff().dt.days
ecom['order_number'] = ecom.groupby('customer_id').cumcount() + 1
print(f'Customers: {ecom["customer_id"].nunique()} Orders: {len(ecom):,}')
print(f'Orders per customer: {ecom["customer_id"].value_counts().min()} to '
f'{ecom["customer_id"].value_counts().max()}')
print()
print('One customer, in order:')
print(ecom[ecom['customer_id'] == 1][
['customer_id', 'order_date', 'order_number', 'days_since_prev_order']].head(5).to_string(index=False))
print()
print(f'First orders have no previous order, so they are NaN: {int(ecom["days_since_prev_order"].isna().sum())}')
print(f'Median gap between orders: {ecom["days_since_prev_order"].median():.0f} days')
Customers: 200 Orders: 2,000
Orders per customer: 3 to 19
One customer, in order:
customer_id order_date order_number days_since_prev_order
1 2023-01-06 1 NaN
1 2023-01-22 2 16.0
1 2023-01-23 3 1.0
1 2023-02-02 4 10.0
1 2023-03-05 5 31.0
First orders have no previous order, so they are NaN: 200
Median gap between orders: 21 daysTwo details in that cell matter. groupby has to come before diff. Without it, a customer's first order would be compared with a different customer's last order, which is meaningless. Whenever you calculate differences across rows, check what happens where one group ends and the next begins.
The resulting NaN is correct. Don't fill it. A first order really has no previous order. Filling it with the median gap would say something false about a customer for whom the question doesn't apply yet. Better options: leave it as NaN for a model that handles it, fill it with a placeholder plus an is_first_order flag, or use a very large value meaning "never before."
The leak recency features invite#
There's a trap here serious enough to deserve its own warning. It also finishes the keep='first' vs keep='last' discussion that notebook 2 started. diff() looks backward, which is safe. As soon as anything looks forward, future information leaks into the past.
customer = ecom[ecom['customer_id'] == 1].head(5).copy()
customer['SAFE_days_since_prev'] = customer['order_date'].diff().dt.days # backward
customer['LEAK_days_until_next'] = -customer['order_date'].diff(-1).dt.days # forward
print(customer[['order_date', 'SAFE_days_since_prev', 'LEAK_days_until_next']].to_string(index=False))
print()
print('The second column answers "when will they order next?".')
print('At prediction time that has not happened yet, so the value does not exist.')
order_date SAFE_days_since_prev LEAK_days_until_next 2023-01-06 NaN 16.0 2023-01-22 16.0 1.0 2023-01-23 1.0 10.0 2023-02-02 10.0 31.0 2023-03-05 31.0 NaN The second column answers "when will they order next?". At prediction time that has not happened yet, so the value does not exist.
days_until_next_order answers "when will this customer order again?" At the moment you need a prediction, that hasn't happened yet, so the value doesn't exist. It would be very good at predicting churn, because it's built from the answer. It would give a great validation score for a useless model.
This is notebook 1's rule in its clearest form: a feature is only allowed if you'd know its value at the moment you make the prediction. With time data, this rule is easy to state and easy to break, because pandas makes looking forward just as easy as looking backward.
The rule reaches further than you might expect. This notebook uses each hour's actual weather. That's fine for learning, but a real forecast for tomorrow can only use tomorrow's weather forecast, because the actual weather isn't known yet.
It also answers notebook 2's duplicate question. If you're rebuilding what was known at a past decision point, keep='last' takes the most recent record, which might not have existed yet at that time. The safe approach is to sort by time and keep the latest record dated on or before the decision date.
Section 5Calendar interactions#
Section 2 showed two columns disagreeing about what an hour means. That disagreement gives one of the biggest gains in this notebook
Section 2's fourth chart showed workingday and hr disagreeing about what an hour means for demand. That's an interaction. Notebook 3 explained why a plain linear model can't find one by itself: it can only add effects together in straight lines, unless you give it the product of the two columns.
hour_dummies = pd.get_dummies(bikes['hr'], prefix='h').astype(int)
with_hour = pd.concat([weather, hour_dummies], axis=1)
interaction = hour_dummies.mul(bikes['workingday'], axis=0).add_prefix('workday_x_')
with_interaction = pd.concat([with_hour, interaction], axis=1)
with_holiday = with_interaction.assign(holiday=bikes['holiday'])
print(f'{"features":40s} {"columns":>8s} {"Ridge R2":>10s}')
print('-' * 60)
print(f'{"hour one-hot + workingday":40s} {with_hour.shape[1]:>8d} {ridge_r2(with_hour):>10.4f}')
print(f'{"+ workingday x hour interaction":40s} {with_interaction.shape[1]:>8d} {ridge_r2(with_interaction):>10.4f}')
print(f'{"+ holiday flag":40s} {with_holiday.shape[1]:>8d} {ridge_r2(with_holiday):>10.4f}')
features columns Ridge R2 ------------------------------------------------------------ hour one-hot + workingday 28 0.6006 + workingday x hour interaction 52 0.7331 + holiday flag 53 0.7336
The interaction raises R-squared from 0.6006 to 0.7331, a gain of 0.13. Only the switch from a raw hour number to one-hot hours added more (0.3243 to 0.6006). It works because the model can finally say: the effect of 8am depends on whether people are commuting that day.
Before, the model had one coefficient for 8am, and it had to cover both a 477-rental workday morning and a 106-rental weekend morning. It split the difference and got both wrong.
print(f'Holiday hours: {int(bikes["holiday"].sum())} of {len(bikes):,} ({bikes["holiday"].mean():.1%})')
print(f'Mean rentals on holidays : {bikes[bikes["holiday"] == 1]["cnt"].mean():.0f}')
print(f'Mean rentals on non-holidays : {bikes[bikes["holiday"] == 0]["cnt"].mean():.0f}')
print()
print('For real projects, the `holidays` package generates these by country and year:')
print(' pip install holidays')
print(' import holidays; us = holidays.US(years=[2011, 2012])')
print(' df["is_holiday"] = df["date"].isin(us).astype(int)')
print()
print('Hardcoding a list is fine for a fixed historical extract and breaks silently')
print('the first time your model runs in a year you did not enumerate.')
Holiday hours: 500 of 17,379 (2.9%)
Mean rentals on holidays : 157
Mean rentals on non-holidays : 190
For real projects, the `holidays` package generates these by country and year:
pip install holidays
import holidays; us = holidays.US(years=[2011, 2012])
df["is_holiday"] = df["date"].isin(us).astype(int)
Hardcoding a list is fine for a fixed historical extract and breaks silently
the first time your model runs in a year you did not enumerate.The holiday flag adds 0.0005, which is basically nothing, and it's better to say so than to hide the result. The numbers explain why: only 500 of the 17,379 hours are holidays, and workingday already contains it. In this dataset, a working day is defined as a day that's neither a weekend nor a holiday, so every holiday is already marked as a non-working day.
A feature that repeats what the model already has can only add what's left over. That's notebook 2's missing-indicator lesson again.
For real projects, the holidays Python package creates these flags automatically by country and year, so you don't need a hand-typed list:
# first, in a terminal: pip install holidays
import holidays
us = holidays.US(years=[2011, 2012])
# .dt.date matters: the package stores plain dates, so comparing
# full timestamps matches nothing and silently gives all zeros
df["is_holiday"] = df["date"].dt.date.isin(us).astype(int)Note the .dt.date. The package stores plain dates, so comparing full timestamps against it matches nothing and quietly gives all zeros, with no error. (The printed hint inside the notebook cell above has that bug; this version is the corrected one.)
A hand-typed list of dates is fine for a fixed historical dataset like this one. But it breaks silently the first time the model runs in a year nobody added to the list.
Section 6Time-based validation#
Every cross-validation in this series so far used a random split, and for four notebooks that was right. Here it stops being right
KFold(shuffle=True) mixes rows randomly into folds, ignoring time order. With time data, that means a model might train on Tuesday and Thursday to predict Wednesday. It's filling in a gap inside a period it has basically already seen, which is much easier than the real task. In production, a model always has to predict a future it hasn't seen. The problem gets worse when the data has a trend.
quarterly = bikes.groupby(bikes['dteday'].dt.to_period('Q'))['cnt'].mean().round(0)
print('Mean rentals per hour, by quarter:')
print(quarterly.to_string())
print()
print(f'2011 average: {bikes[bikes["year"] == 2011]["cnt"].mean():.0f}')
print(f'2012 average: {bikes[bikes["year"] == 2012]["cnt"].mean():.0f}')
print(f'Growth: {bikes[bikes["year"] == 2012]["cnt"].mean() / bikes[bikes["year"] == 2011]["cnt"].mean() - 1:.0%}')
Mean rentals per hour, by quarter: dteday 2011Q1 73.0 2011Q2 171.0 2011Q3 185.0 2011Q4 142.0 2012Q1 168.0 2012Q2 263.0 2012Q3 288.0 2012Q4 219.0 Freq: Q-DEC 2011 average: 144 2012 average: 235 Growth: 63%
Ridership grew 63% from 2011 to 2012. A model validated by mixing both years in every fold never has to deal with that growth, because examples from the busier future are already in its training data.
# Sort by date AND hour. Sorting on 'dteday' alone is not enough: pandas'
# default sort is not stable, so rows within a day come back in arbitrary order.
bikes_sorted = bikes.sort_values(['dteday', 'hr']).reset_index(drop=True)
X_sorted = pd.concat([
bikes_sorted[['temp', 'hum', 'windspeed', 'workingday']],
pd.get_dummies(bikes_sorted['hr'], prefix='h').astype(int),
], axis=1)
y_sorted = bikes_sorted['cnt']
random_scores = cross_val_score(make_pipeline(StandardScaler(), Ridge()), X_sorted, y_sorted,
cv=KFold(5, shuffle=True, random_state=RANDOM_STATE), scoring='r2')
forward_scores = cross_val_score(make_pipeline(StandardScaler(), Ridge()), X_sorted, y_sorted,
cv=TimeSeriesSplit(n_splits=5), scoring='r2')
is_2011 = bikes_sorted['year'] == 2011
holdout = make_pipeline(StandardScaler(), Ridge()).fit(X_sorted[is_2011], y_sorted[is_2011])
holdout_r2 = r2_score(y_sorted[~is_2011], holdout.predict(X_sorted[~is_2011]))
print(f'KFold(shuffle=True) R2 = {random_scores.mean():.4f} folds {np.round(random_scores, 3)}')
print(f'TimeSeriesSplit R2 = {forward_scores.mean():.4f} folds {np.round(forward_scores, 3)}')
print(f'Train 2011, test 2012 R2 = {holdout_r2:.4f}')
KFold(shuffle=True) R2 = 0.6006 folds [0.6 0.596 0.585 0.616 0.607] TimeSeriesSplit R2 = 0.5072 folds [0.524 0.563 0.472 0.506 0.471] Train 2011, test 2012 R2 = 0.4096
Three numbers for the same model and features. Only the way the rows were split changes.
Random KFold reports 0.6006. The model has seen hours from every week of both years, so predicting a held-out hour is mostly filling in between neighbors it already knows.
TimeSeriesSplit reports 0.5072. This trains on an early block, tests on the next block, then extends the training window and repeats. Each test fold is truly in the future compared with its training data, which is what happens in production.
Training on 2011 and testing on 2012 gives 0.4096. This is the hardest and most realistic of the three, because ridership grew 63% and the model never saw 2012-level demand in training.
The gap between 0.6006 and 0.4096 isn't noise. It's the difference between a number that flatters a model and a number you could defend to someone making a decision. Reporting the first number for a forecasting problem is a real and common mistake.
What this means in practice: use TimeSeriesSplit whenever rows have a meaningful order and you're predicting forward. Never shuffle time data before splitting. Watch for trends. When the overall level shifts like it does here, a model trained on the past will under-predict the future.
A steadily growing feature like days_since_launch can help a linear model here, because it lets the model continue a trend in a straight line. It doesn't help a tree model, since trees can't predict beyond the range of values they were trained on.
A related rule applies to grouped data, like several orders from each customer. Split by customer, using GroupKFold, so the same person never appears in both training and test.
A fair objection#
Sections 3 and 5 compared encodings using KFold(shuffle=True), and this section just argued that's wrong for time data. That's a fair objection to this notebook's own earlier results, and it deserves a direct answer.
Were those earlier numbers inflated? Yes, all of them. But the real question is whether the ranking of the encodings changes, because every earlier conclusion was a comparison between options, not a claim about one exact score.
X_sorted_hour = bikes_sorted[['temp', 'hum', 'windspeed', 'workingday']]
hr_sorted = bikes_sorted['hr']
sin_a, cos_a = cyclical(hr_sorted, 24, 1)
sin_b, cos_b = cyclical(hr_sorted, 24, 2)
oh_sorted = pd.get_dummies(hr_sorted, prefix='h').astype(int)
inter_sorted = oh_sorted.mul(bikes_sorted['workingday'], axis=0).add_prefix('workday_x_')
to_check = {
'raw integer hour': X_sorted_hour.assign(hr=hr_sorted),
'cyclical, 1 harmonic': X_sorted_hour.assign(s1=sin_a, c1=cos_a),
'cyclical, 2 harmonics': X_sorted_hour.assign(s1=sin_a, c1=cos_a, s2=sin_b, c2=cos_b),
'one-hot, 24 hours': pd.concat([X_sorted_hour, oh_sorted], axis=1),
'one-hot + interaction': pd.concat([X_sorted_hour, oh_sorted, inter_sorted], axis=1),
}
print(f'{"encoding":24s} {"random KFold":>14s} {"TimeSeriesSplit":>17s}')
print('-' * 58)
for name, X in to_check.items():
random_r2 = cross_val_score(make_pipeline(StandardScaler(), Ridge()), X, y_sorted,
cv=KFold(5, shuffle=True, random_state=RANDOM_STATE), scoring='r2').mean()
forward_r2 = cross_val_score(make_pipeline(StandardScaler(), Ridge()), X, y_sorted,
cv=TimeSeriesSplit(n_splits=5), scoring='r2').mean()
print(f'{name:24s} {random_r2:>14.4f} {forward_r2:>17.4f}')
encoding random KFold TimeSeriesSplit ---------------------------------------------------------- raw integer hour 0.3243 0.2055
cyclical, 1 harmonic 0.4236 0.3165
cyclical, 2 harmonics 0.5116 0.4100
one-hot, 24 hours 0.6006 0.5072
one-hot + interaction 0.7331 0.6517
Every score drops by about 0.08 to 0.12 when the split respects time. But the order of all five encodings is the same in both columns: raw number worst, harmonics improving on a single wave, one-hot beating them all, and the interaction beating everything. So the conclusions from sections 3 and 5 hold. Claims about the exact scores wouldn't have.
That's the key point: random cross-validation on time data can still rank options correctly while overstating how well each one will actually do. It isn't guaranteed, though. Features built from nearby rows, like the previous hour's rentals, get a big unfair boost from random splits, because the neighboring hours are sitting in the training data. With features like that, even the ranking can flip.
That's why this mistake survives in so many projects. It answers "which option is better?" well enough that nobody checks "how accurate will this be?" until someone asks what to promise. Using the random split earlier was a deliberate simplification to keep those comparisons easy to read, and this section settles that debt. In your own work, once you know the data has a time order, use TimeSeriesSplit from the start instead of adding it later.
Section 7Putting it in a function#
Building these features by hand each time leads to differences between training and serving, which is its own kind of production bug
One function that defines all these features in one place removes a whole type of bug: training calculated dayofweek one way, and a serving script written weeks later by someone else calculated it slightly differently.
def datetime_features(df, column, prefix=None, cycles=('month', 'dayofweek')):
# Derive standard calendar features from one datetime column.
# Returns a new DataFrame; does not modify the input.
prefix = prefix or column
ts = pd.to_datetime(df[column])
out = pd.DataFrame(index=df.index)
out[f'{prefix}_year'] = ts.dt.year
out[f'{prefix}_month'] = ts.dt.month
out[f'{prefix}_day'] = ts.dt.day
out[f'{prefix}_dayofweek'] = ts.dt.dayofweek
out[f'{prefix}_dayofyear'] = ts.dt.dayofyear
out[f'{prefix}_quarter'] = ts.dt.quarter
out[f'{prefix}_is_weekend'] = (ts.dt.dayofweek >= 5).astype(int)
out[f'{prefix}_is_month_start'] = ts.dt.is_month_start.astype(int)
out[f'{prefix}_is_month_end'] = ts.dt.is_month_end.astype(int)
# dayofyear uses 365 and is therefore very slightly off in leap years,
# which matters for almost nothing but is worth knowing.
periods = {'month': 12, 'dayofweek': 7, 'hour': 24, 'dayofyear': 365}
for name in cycles:
values = getattr(ts.dt, name)
angle = 2 * np.pi * values / periods[name]
out[f'{prefix}_{name}_sin'] = np.sin(angle)
out[f'{prefix}_{name}_cos'] = np.cos(angle)
return out
generated = datetime_features(orders, 'order_date')
print(f'One datetime column produced {generated.shape[1]} features:')
print(list(generated.columns))
print()
print(generated.head(3).round(3).to_string())
One datetime column produced 13 features: ['order_date_year', 'order_date_month', 'order_date_day', 'order_date_dayofweek', 'order_date_dayofyear', 'order_date_quarter', 'order_date_is_weekend', 'order_date_is_month_start', 'order_date_is_month_end', 'order_date_month_sin', 'order_date_month_cos', 'order_date_dayofweek_sin', 'order_date_dayofweek_cos'] order_date_year order_date_month order_date_day order_date_dayofweek order_date_dayofyear order_date_quarter order_date_is_weekend order_date_is_month_start order_date_is_month_end order_date_month_sin order_date_month_cos order_date_dayofweek_sin order_date_dayofweek_cos 0 2024 5 10 4 131 2 0 0 0 0.5 -0.866 -0.434 -0.901 1 2024 11 10 6 315 4 1 0 0 -0.5 0.866 -0.782 0.623 2 2022 5 2 0 122 2 0 0 0 0.5 -0.866 0.000 1.000
Two choices in that function are worth copying.
It returns a new DataFrame instead of changing the input. A function that changes its input makes the notebook depend on run order. Running the same cell twice could give different results, depending on the input's current state.
The cyclical columns are optional, using the cycles argument. Section 3 showed cyclical encoding helps some models and hurts others, so adding it every time would treat one modeling choice as a universal default.
The function also doesn't create every possible feature. Notebook 3 showed 83 polynomial columns doing worse than a smaller set of 27, and the same restraint applies here. Build what you have a reason to think will help, then measure whether it did. Notebook 8 covers removing features that don't earn their place.
Wrap-upKey things to remember from this notebook#
Parse dates explicitly, and pass format= whenever the source isn't ISO. 03/04/2023 is a valid date two different ways and neither raises an error, so an unclear format corrupts data silently. errors='coerce' turns bad values into NaT, which is only safe once you've counted them.
Store timestamps in UTC, and convert to local time before extracting anything people care about. One moment is hour 1 in New York and hour 12 in Kolkata, so "hour of day" means nothing without a time zone, and rush hour happens in local time.
Cyclical encoding solves a real problem, but it isn't always the best choice. It made hours 23 and 0 as close as 23 and 22, both 0.261 apart, and beat the raw number for Ridge regression, 0.4236 against 0.3243. One-hot then reached 0.6006, because one sine and cosine pair can only show one peak per cycle and bike demand has two. A second harmonic reached 0.5116 and a third 0.5452, matching how much of the daily pattern each can capture: 67.4%, 82.2%, 88.9%.
The best encoding depends on the model, not just the data. The random forest ranked the same three encodings in reverse: raw number best at 0.8058, cyclical 0.8030 and one-hot 0.7874, because trees split on thresholds instead of measuring distance. That's notebook 4's label-encoding lesson again.
Durations outlast dates. A model trained on signup_year has to guess once the calendar passes the training years, while account_age_days stays in a familiar range. A count from one fixed date, like days_since_launch, keeps growing and has the same problem as a year. Use groupby before diff so one customer's gap never mixes with another's, and treat the NaN on a first order as correct, not as missing data.
Only look backward. days_since_previous_order is fine. days_until_next_order is the answer in disguise, and pandas makes both equally easy to write. A feature is only allowed if you'd know its value when the prediction is made, and that includes things like the weather.
Calendar interactions were the biggest gain after encoding the hour itself. Multiplying the hour columns by workingday raised R-squared from 0.6006 to 0.7331, because 8am means 477 rentals on a working day and only 106 otherwise, and one shared coefficient can't cover both. The holiday flag added only 0.0005, since workingday already marks every holiday as a non-working day.
Never shuffle time data before splitting. Random KFold reported 0.6006, TimeSeriesSplit reported 0.5072, and training on 2011 to predict 2012 gave 0.4096. Ridership grew 63% between the two years, and only that last number reflects what a deployed model would face.
Notebook 6 asks the same kind of question about a different input: free text instead of a timestamp. The column isn't a number at all, and the first job is deciding what in a sentence is worth counting.
This article goes with Notebook 5 of a feature engineering series built on the UCI Bike Sharing dataset, with a small e-commerce table for recency features. Every code cell above runs exactly as shown in the notebook, with its real output.