⏳ Time Series Analysis

📚 What You'll Learn

By the end of this lesson, you will be able to:

Create and work with a DatetimeIndex using pd.to_datetime() and pd.date_range()
Change a series' frequency with resample() and asfreq() (downsampling and upsampling)
Smooth and summarize trends with rolling windows
Shift data in time with shift(), diff(), and pct_change() to build lags and growth rates
Decompose a series into trend, seasonality, and residual components
Forecast future values and score them with error metrics like MAE, RMSE, MAPE, and R²

⏱️ Estimated Time: 55–70 minutes

🎯 Project: Take four years of daily data, resample it to monthly, smooth it with rolling means, build several baseline forecasts, and rank them with proper time-series error metrics.

Work with data that has a clock attached — index by time, change frequency, smooth trends, and forecast what comes next.

🌟 Data That Moves Through Time

Sales per day, temperature per hour, a heartbeat per second, a stock price per tick — an enormous amount of the world's data is ordered by time. That ordering is not just metadata; it is the signal. Yesterday informs today, last December looks like this December, and the gap between two readings is itself information.

Pandas was built with time series in its bones. Once your index is a DatetimeIndex, a whole vocabulary opens up: slice by '2024-03', resample daily data into monthly totals, take a 7-day rolling average, or lag a column to compare it against its own past. In this lesson you will move from raw timestamps all the way to a scored forecast.

🗓️ Building a Datetime Index

Everything starts by telling pandas that a column is time. Once it is the index, time-aware slicing and resampling just work.

import pandas as pd
import numpy as np

# Parse strings into real timestamps
df['date'] = pd.to_datetime(df['date'])

# Make the timestamp the index
df = df.set_index('date').sort_index()

# Generate a regular range of dates
idx = pd.date_range('2024-01-01', periods=365, freq='D')
ts = pd.Series(np.random.randn(365).cumsum() + 100, index=idx)

# Time-aware slicing — no boolean masks needed
ts['2024-03']            # every day in March 2024
ts['2024-06':'2024-08']  # the summer quarter
ts.between_time('09:00', '17:00')  # intraday selection

📈 Interactive Time Series Explorer

Pick a date range and a sampling frequency to see how the same underlying signal reads at different resolutions.

From To

Tip: change the inputs and watch the same trend + seasonality + noise recombine at each frequency.

🔁 Resampling: Changing Frequency

resample() is a groupby for time. Downsampling (daily → monthly) aggregates; upsampling (monthly → daily) fills.

⬇️ Downsampling

Fewer, coarser buckets — you must aggregate.

# Daily -> monthly totals
ts.resample('M').sum()

# Daily -> weekly mean
ts.resample('W').mean()

# Open-high-low-close per month
ts.resample('M').ohlc()

⬆️ Upsampling

More, finer buckets — you must fill the gaps.

# Monthly -> daily, forward fill
monthly.resample('D').ffill()

# Interpolate between points
monthly.resample('D').interpolate()

# Change frequency without aggregating
ts.asfreq('H')

🧩 Seasonal Decomposition

Most real series are a sum of parts: a slow trend, a repeating seasonal pattern, and leftover residual noise. Separating them makes each easier to understand and model.

from statsmodels.tsa.seasonal import seasonal_decompose

result = seasonal_decompose(ts, model='additive', period=365)
result.trend      # long-run direction
result.seasonal   # repeating cycle
result.resid      # what's left over
result.plot()

Trend

Seasonal

Residual

🪟 Rolling Windows

A rolling window slides a fixed-size frame across the series, computing a statistic at each step. It is the classic tool for smoothing noise and revealing the trend underneath. Drag the slider to widen the window and watch the smoothed line calm down.

Window size 20 days

# 7-day rolling mean smooths daily noise
ts.rolling(window=7).mean()

# Rolling std as a volatility measure
ts.rolling(window=30).std()

# Center the window on each point
ts.rolling(window=7, center=True).mean()

# Require a minimum number of observations
ts.rolling(window=30, min_periods=10).mean()

↔️ Shifting, Lags & Growth

Moving a series forward or backward in time lets you compare each point against its own past — the foundation of growth rates, momentum, and lag features for models.

# Yesterday's value alongside today's
df['prev'] = ts.shift(1)

# Day-over-day change
df['change'] = ts.diff()

# Percentage growth
df['growth'] = ts.pct_change() * 100

# Year-over-year comparison (365-day lag)
df['yoy'] = ts.pct_change(periods=365) * 100

🚨 Spotting Anomalies

Points that stray far from a rolling mean — say, more than a few rolling standard deviations away — are candidate anomalies worth investigating.

roll_mean = ts.rolling(30).mean()
roll_std  = ts.rolling(30).std()
upper = roll_mean + 3 * roll_std
lower = roll_mean - 3 * roll_std
anomalies = ts[(ts > upper) | (ts < lower)]

🔮 Forecasting & Evaluation

Once you understand a series you can project it forward. Always hold out the most recent slice as a test set (never shuffle time!), forecast it, and score the forecast with error metrics. Lower error and higher R² are better.

Naive / Seasonal Naive

Repeat the last value (or last season). The baseline every model must beat.

Exponential Smoothing

Weighted averages of the past; Holt-Winters adds trend and seasonality.

ARIMA / SARIMA

Model autocorrelation and differencing; SARIMA adds seasonal terms.

Best model on the test set

9.41

MAE

12.89

RMSE

2.8%

MAPE

0.92

R²

Time Series Forecasting Methods


import pandas as pd
import numpy as np
from statsmodels.tsa.holtwinters import ExponentialSmoothing
from statsmodels.tsa.arima.model import ARIMA
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')

# Create sample time series
np.random.seed(42)
dates = pd.date_range('2020-01-01', '2023-12-31', freq='D')
trend = np.linspace(100, 400, len(dates))
seasonal = 50 * np.sin(2 * np.pi * np.arange(len(dates)) / 365.25)
noise = np.random.randn(len(dates)) * 10
ts = pd.Series(trend + seasonal + noise, index=dates)

# Split into train and test
train_size = int(len(ts) * 0.8)
train, test = ts[:train_size], ts[train_size:]

print(f"Train set: {len(train)} observations")
print(f"Test set: {len(test)} observations")

# Forecasting Methods
#####################

# 1. Naive Methods
##################

# Last value (Naive)
naive_forecast = pd.Series([train.iloc[-1]] * len(test), index=test.index)

# Seasonal Naive
seasonal_period = 365
seasonal_naive = train.iloc[-seasonal_period:].values
seasonal_naive_forecast = pd.Series(
    np.tile(seasonal_naive, len(test) // seasonal_period + 1)[:len(test)],
    index=test.index
)

# Average method
average_forecast = pd.Series([train.mean()] * len(test), index=test.index)

# 2. Moving Average Forecast
#############################

def moving_average_forecast(train, test, window=30):
    """Simple moving average forecast"""
    forecast = []
    history = list(train)

    for _ in range(len(test)):
        ma = np.mean(history[-window:])
        forecast.append(ma)
        # In practice, would add actual observation
        history.append(ma)

    return pd.Series(forecast, index=test.index)

ma_forecast = moving_average_forecast(train, test)

# 3. Exponential Smoothing
##########################

# Simple Exponential Smoothing
from statsmodels.tsa.holtwinters import SimpleExpSmoothing

ses_model = SimpleExpSmoothing(train)
ses_fit = ses_model.fit(smoothing_level=0.2, optimized=False)
ses_forecast = ses_fit.forecast(len(test))

# Double Exponential Smoothing (Holt's method)
from statsmodels.tsa.holtwinters import Holt

holt_model = Holt(train)
holt_fit = holt_model.fit()
holt_forecast = holt_fit.forecast(len(test))

# Triple Exponential Smoothing (Holt-Winters)
hw_model = ExponentialSmoothing(
    train,
    seasonal_periods=365,
    trend='add',
    seasonal='add',
    damped_trend=True
)
hw_fit = hw_model.fit()
hw_forecast = hw_fit.forecast(len(test))

# 4. ARIMA Models
#################

# Auto ARIMA (find best parameters)
from statsmodels.tsa.stattools import adfuller
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

# Check stationarity and difference if needed
adf_result = adfuller(train)
if adf_result[1] > 0.05:
    train_diff = train.diff().dropna()
    d = 1
else:
    train_diff = train
    d = 0

# Fit ARIMA model
arima_model = ARIMA(train, order=(2, d, 2))  # (p, d, q)
arima_fit = arima_model.fit()
arima_forecast = arima_fit.forecast(steps=len(test))

# SARIMA (Seasonal ARIMA)
from statsmodels.tsa.statespace.sarimax import SARIMAX

sarima_model = SARIMAX(
    train,
    order=(1, 1, 1),
    seasonal_order=(1, 1, 1, 365),
    enforce_stationarity=False,
    enforce_invertibility=False
)
sarima_fit = sarima_model.fit(disp=False)
sarima_forecast = sarima_fit.forecast(steps=len(test))

# 5. Prophet (Facebook's forecasting tool)
###########################################

from prophet import Prophet

# Prepare data for Prophet
prophet_train = pd.DataFrame({
    'ds': train.index,
    'y': train.values
})

# Fit Prophet model
prophet_model = Prophet(
    yearly_seasonality=True,
    weekly_seasonality=False,
    daily_seasonality=False,
    changepoint_prior_scale=0.05
)
prophet_model.fit(prophet_train)

# Make future dataframe
future = prophet_model.make_future_dataframe(periods=len(test), freq='D')
prophet_pred = prophet_model.predict(future)
prophet_forecast = prophet_pred.loc[prophet_pred['ds'].isin(test.index), 'yhat'].values

# Model Evaluation
##################

def evaluate_forecast(actual, forecast, model_name):
    """Calculate forecast metrics"""
    mae = mean_absolute_error(actual, forecast)
    rmse = np.sqrt(mean_squared_error(actual, forecast))
    mape = np.mean(np.abs((actual - forecast) / actual)) * 100
    r2 = r2_score(actual, forecast)

    return {
        'Model': model_name,
        'MAE': mae,
        'RMSE': rmse,
        'MAPE': mape,
        'R2': r2
    }

# Evaluate all models
results = []
models = {
    'Naive': naive_forecast,
    'Seasonal Naive': seasonal_naive_forecast[:len(test)],
    'Moving Average': ma_forecast,
    'Simple Exp Smoothing': ses_forecast,
    'Holt': holt_forecast,
    'Holt-Winters': hw_forecast,
    'ARIMA': arima_forecast
}

for name, forecast in models.items():
    if len(forecast) == len(test):
        results.append(evaluate_forecast(test.values, forecast.values, name))

results_df = pd.DataFrame(results)
results_df = results_df.sort_values('RMSE')
print("\nModel Performance Comparison:")
print(results_df.to_string(index=False))

# Visualization
###############

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Plot 1: All forecasts comparison
ax1 = axes[0, 0]
train.plot(ax=ax1, label='Training Data', color='blue', alpha=0.7)
test.plot(ax=ax1, label='Actual', color='black', linewidth=2)
ma_forecast.plot(ax=ax1, label='Moving Average', alpha=0.7)
holt_forecast.plot(ax=ax1, label='Holt', alpha=0.7)
hw_forecast.plot(ax=ax1, label='Holt-Winters', alpha=0.7)
ax1.set_title('Forecast Comparison')
ax1.legend(loc='upper left')
ax1.grid(True, alpha=0.3)

# Plot 2: Best model forecast with confidence intervals
ax2 = axes[0, 1]
best_model_name = results_df.iloc[0]['Model']
best_forecast = models[best_model_name]

train.plot(ax=ax2, label='Training Data', color='blue', alpha=0.7)
test.plot(ax=ax2, label='Actual', color='black', linewidth=2)
best_forecast.plot(ax=ax2, label=f'Best: {best_model_name}', color='red', linewidth=2)

# Add confidence intervals (example for ARIMA)
if best_model_name == 'ARIMA':
    forecast_result = arima_fit.get_forecast(steps=len(test))
    conf_int = forecast_result.conf_int()
    ax2.fill_between(test.index,
                     conf_int.iloc[:, 0],
                     conf_int.iloc[:, 1],
                     alpha=0.2, color='red')

ax2.set_title(f'Best Model: {best_model_name}')
ax2.legend()
ax2.grid(True, alpha=0.3)

# Plot 3: Residuals analysis
ax3 = axes[1, 0]
residuals = test - best_forecast
residuals.plot(ax=ax3, label='Residuals', color='green')
ax3.axhline(y=0, color='red', linestyle='--', alpha=0.5)
ax3.set_title('Forecast Residuals')
ax3.set_ylabel('Residual')
ax3.legend()
ax3.grid(True, alpha=0.3)

# Plot 4: Residuals distribution
ax4 = axes[1, 1]
residuals.hist(ax=ax4, bins=30, edgecolor='black', alpha=0.7)
ax4.axvline(x=0, color='red', linestyle='--', linewidth=2)
ax4.set_title('Residuals Distribution')
ax4.set_xlabel('Residual')
ax4.set_ylabel('Frequency')

# Add normal distribution overlay
from scipy import stats
mu, std = residuals.mean(), residuals.std()
xmin, xmax = ax4.get_xlim()
x = np.linspace(xmin, xmax, 100)
p = stats.norm.pdf(x, mu, std) * len(residuals) * (xmax - xmin) / 30
ax4.plot(x, p, 'k-', linewidth=2, label='Normal')
ax4.legend()

plt.tight_layout()
plt.show()

# Cross-Validation for Time Series
##################################

from sklearn.model_selection import TimeSeriesSplit

def time_series_cv(ts, model_func, n_splits=5):
    """
    Perform time series cross-validation
    """
    tscv = TimeSeriesSplit(n_splits=n_splits)
    scores = []

    for train_idx, test_idx in tscv.split(ts):
        train_cv = ts.iloc[train_idx]
        test_cv = ts.iloc[test_idx]

        # Apply model function
        forecast_cv = model_func(train_cv, len(test_cv))

        # Calculate RMSE
        rmse = np.sqrt(mean_squared_error(test_cv, forecast_cv))
        scores.append(rmse)

    return np.mean(scores), np.std(scores)

# Example: Cross-validate moving average
def ma_model(train, horizon):
    return pd.Series([train.rolling(30).mean().iloc[-1]] * horizon)

cv_mean, cv_std = time_series_cv(train, ma_model)
print(f"\nCross-validation RMSE: {cv_mean:.2f} ± {cv_std:.2f}")

✅ Time Series Best Practices

Always visualize first: Plot your data to identify patterns, trends, and anomalies
Check for stationarity: Many models require stationary data
Handle missing values: Use appropriate interpolation methods
Consider seasonality: Account for regular patterns in your data
Use proper train-test split: Respect temporal order, never shuffle
Validate on future data: Use walk-forward validation for realistic assessment
Monitor for regime changes: Time series properties can change over time
Combine multiple models: Ensemble methods often perform better
Document preprocessing: Keep track of all transformations applied
Consider external factors: Include relevant exogenous variables when available

📋 Quick Reference

Essential Functions:

pd.date_range() - Create date ranges
pd.to_datetime() - Convert to datetime
df.resample() - Change frequency
df.rolling() - Rolling window operations
df.shift() - Lag operations
df.diff() - Differencing
df.expanding() - Expanding window
df.ewm() - Exponential weighted functions
df.asfreq() - Convert to specified frequency
df.interpolate() - Fill missing values

Frequency Strings:

'D' - Daily
'B' - Business day
'W' - Weekly
'M' - Month end
'MS' - Month start
'Q' - Quarter end
'Y' - Year end
'H' - Hourly
'T' or 'min' - Minutely
'S' - Secondly

📓 Learning Journal

Keep a learning journal — digital or physical. After this lesson, take a few minutes to write down:

Key concepts you learned
Techniques that clicked for you
Questions or confusion points to revisit
Ideas you want to try
Your progress and feelings about learning this

✍️ This lesson's prompt: Think of something in your own life that is recorded over time — steps walked, money spent, hours studied. What frequency is it naturally recorded at, and what would you learn by resampling it to a coarser one or smoothing it with a rolling mean? Which pattern do you suspect is trend and which is seasonality?

📝 Lesson Summary

🎓 Key Takeaways

A DatetimeIndex unlocks time-aware slicing, resampling, and rolling operations — set it and sort it first.
resample() is a groupby over time: downsampling aggregates, upsampling fills.
Rolling windows smooth noise to reveal trend; shift(), diff(), and pct_change() compare each point to its own past.
Forecasts must be tested on a held-out future slice (never shuffled) and scored with metrics like MAE, RMSE, MAPE, and R².

🎉 What You've Accomplished

You can now take a raw stream of timestamped data, reshape its frequency, expose its trend and seasonality, engineer lag and growth features, and produce a forecast you can actually defend with error metrics. That is the full arc of applied time-series analysis.

❓ Common Questions at This Stage

What's the difference between resample() and rolling()?

resample() changes the frequency of the series — it re-buckets the timeline (e.g. daily into monthly) and returns fewer rows. rolling() keeps the same frequency and row count but computes a statistic over a sliding window at each point. Use resample to change resolution, rolling to smooth.

Which forecast error metric should I report?

Report more than one. MAE is in the data's own units and easy to explain; RMSE punishes large misses harder; MAPE gives a scale-free percentage (but blows up near zero); R² shows how much variance you captured. The right primary metric depends on whether big errors or percentage errors hurt your use case more.

Why must I never shuffle time series before splitting?

Because it leaks the future into the past. If training rows come from after your test rows, the model "sees" information it would never have at prediction time, and your scores become fantasy. Always split by time: train on the earlier portion, test on the later portion, and consider walk-forward validation.

🔭 Looking Ahead

Rolling means are one member of a whole family of window operations. Next you'll generalize them — expanding windows, exponentially weighted windows, and window-scoped aggregations — so you can compute exactly the moving statistic each problem calls for.

✅ Before the Next Lesson

Load a time-indexed dataset, resample it to two different frequencies, and compare the plots.
Add a 7-period and 30-period rolling mean and describe how the smoothing differs.
Build a naive forecast on a held-out tail and compute its MAE and RMSE.
Write your Learning Journal entry for this lesson.

🌟 Encouragement for the Journey

Time series is where data science meets prediction — the moment your analysis stops describing the past and starts anticipating the future. It can feel like a lot of moving parts, but every one you just met is reusable. Keep the timeline in view, and keep forecasting!