This tutorial provides a comprehensive introduction to time series analysis using Python, covering the same topics as the original R tutorial but using modern Python libraries like pandas, numpy, matplotlib, statsmodels, and scikit-learn.
Required Libraries
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.holtwinters import ExponentialSmoothing
from statsmodels.tsa.stattools import acf, pacf
from statsmodels.stats.diagnostic import acorr_ljungbox
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import warnings
warnings.filterwarnings('ignore')
# Set style for better plots
plt.style.use('seaborn-v0_8')
sns.set_palette("husl")
1. Reading and Creating Time Series Data
Reading Time Series from Files
# Example: Reading rainfall data (equivalent to R's read.table)
def read_time_series_data(url_or_path, delimiter=None):
"""
Read time series data from URL or file
"""
try:
if delimiter:
data = pd.read_csv(url_or_path, delimiter=delimiter, header=None)
else:
data = pd.read_csv(url_or_path, header=None)
return data.iloc[:, 0].values # Return as numpy array
except:
print(f"Error reading data from {url_or_path}")
return None
# Example data: Kings of England death ages
kings_ages = [60, 43, 67, 50, 56, 42, 50, 65, 68, 43, 65, 34, 47, 34, 49, 41, 13, 35, 53, 56, 16, 43, 69, 59, 48, 59, 86, 55, 68, 51, 33, 49, 67, 77, 81, 67, 71, 81, 68, 70, 77, 56]
# Convert to pandas time series
def create_time_series(data, start_year=1066, freq='A'):
"""
Create a pandas time series object
Parameters:
data: list or array of values
start_year: starting year
freq: frequency ('A' for annual, 'M' for monthly, 'Q' for quarterly)
"""
dates = pd.date_range(start=str(start_year), periods=len(data), freq=freq)
return pd.Series(data, index=dates)
# Create time series for kings data
kings_ts = create_time_series(kings_ages, start_year=1066, freq='A')
print("Kings time series:")
print(kings_ts.head(10))
Creating Monthly and Quarterly Time Series
# Example: Monthly data
monthly_data = np.random.randn(60) + 10 # 5 years of monthly data
monthly_ts = create_time_series(monthly_data, start_year=2019, freq='M')
# Example: Quarterly data
quarterly_data = [4.8, 4.1, 6.0, 6.5, 5.8, 5.2, 6.8, 7.4, 6.0, 5.7, 7.7, 8.0]
quarterly_ts = create_time_series(quarterly_data, start_year=2020, freq='Q')
print("\nMonthly time series (first 12 values):")
print(monthly_ts.head(12))
print("\nQuarterly time series:")
print(quarterly_ts)
2. Plotting Time Series
Basic Time Series Plots
def plot_time_series(ts, title="Time Series Plot", ylabel="Value", figsize=(12, 6)):
"""
Plot a time series with proper formatting
"""
fig, ax = plt.subplots(figsize=figsize)
ts.plot(ax=ax, linewidth=2)
ax.set_title(title, fontsize=16, fontweight='bold')
ax.set_ylabel(ylabel, fontsize=12)
ax.set_xlabel("Time", fontsize=12)
ax.grid(True, alpha=0.3)
plt.tight_layout()
return fig, ax
# Plot the kings data
fig, ax = plot_time_series(kings_ts,
title="Age at Death of English Kings",
ylabel="Age at Death (years)")
plt.show()
# Plot with trend line
fig, ax = plot_time_series(kings_ts,
title="Age at Death of English Kings (with trend)",
ylabel="Age at Death (years)")
# Add trend line
x = np.arange(len(kings_ts))
z = np.polyfit(x, kings_ts.values, 1)
p = np.poly1d(z)
ax.plot(kings_ts.index, p(x), "r--", alpha=0.8, linewidth=2, label='Trend')
ax.legend()
plt.show()
3. Smoothing Time Series
Moving Averages
def simple_moving_average(ts, window):
"""
Calculate simple moving average
"""
return ts.rolling(window=window, center=True).mean()
# Apply different moving averages to kings data
fig, axes = plt.subplots(2, 2, figsize=(15, 10))
axes = axes.flatten()
windows = [3, 5, 8, 12]
for i, window in enumerate(windows):
ax = axes[i]
# Plot original data
kings_ts.plot(ax=ax, alpha=0.7, label='Original', color='lightblue')
# Plot smoothed data
smoothed = simple_moving_average(kings_ts, window)
smoothed.plot(ax=ax, label=f'MA({window})', linewidth=2)
ax.set_title(f'Simple Moving Average (window={window})')
ax.set_ylabel('Age at Death')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
Exponential Smoothing
def exponential_smoothing(ts, alpha=0.3):
"""
Simple exponential smoothing
"""
result = []
result.append(ts.iloc[0]) # First value
for i in range(1, len(ts)):
result.append(alpha * ts.iloc[i] + (1 - alpha) * result[i-1])
return pd.Series(result, index=ts.index)
# Compare different alpha values
fig, ax = plt.subplots(figsize=(12, 8))
kings_ts.plot(ax=ax, alpha=0.7, label='Original', color='lightblue')
alphas = [0.1, 0.3, 0.7, 0.9]
colors = ['red', 'green', 'orange', 'purple']
for alpha, color in zip(alphas, colors):
smoothed = exponential_smoothing(kings_ts, alpha)
smoothed.plot(ax=ax, label=f'α={alpha}', color=color, linewidth=2)
ax.set_title('Exponential Smoothing with Different Alpha Values')
ax.set_ylabel('Age at Death')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
4. Decomposing Time Series
Seasonal Decomposition
# Create sample seasonal data
np.random.seed(42)
dates = pd.date_range('2015-01-01', periods=96, freq='M')
trend = np.linspace(100, 150, 96)
seasonal = 10 * np.sin(2 * np.pi * np.arange(96) / 12)
noise = np.random.normal(0, 5, 96)
seasonal_ts = pd.Series(trend + seasonal + noise, index=dates)
def decompose_time_series(ts, model='additive', period=None):
"""
Decompose time series into trend, seasonal, and residual components
"""
if period is None:
# Try to infer period from frequency
freq = pd.infer_freq(ts.index)
if freq and 'M' in freq:
period = 12
elif freq and 'Q' in freq:
period = 4
else:
period = 12 # Default
decomposition = seasonal_decompose(ts, model=model, period=period)
return decomposition
# Perform decomposition
decomp = decompose_time_series(seasonal_ts, model='additive')
# Plot decomposition
fig, axes = plt.subplots(4, 1, figsize=(15, 12))
decomp.observed.plot(ax=axes[0], title='Original Time Series')
axes[0].set_ylabel('Value')
decomp.trend.plot(ax=axes[1], title='Trend Component', color='red')
axes[1].set_ylabel('Trend')
decomp.seasonal.plot(ax=axes[2], title='Seasonal Component', color='green')
axes[2].set_ylabel('Seasonal')
decomp.resid.plot(ax=axes[3], title='Residual Component', color='orange')
axes[3].set_ylabel('Residual')
for ax in axes:
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# Print seasonal factors
print("Seasonal factors by month:")
seasonal_factors = decomp.seasonal.groupby(decomp.seasonal.index.month).mean()
for month, factor in seasonal_factors.items():
print(f"Month {month:2d}: {factor:6.2f}")
Seasonally Adjusted Time Series
def seasonal_adjustment(ts, method='additive'):
"""
Remove seasonal component from time series
"""
decomp = decompose_time_series(ts, model=method)
if method == 'additive':
adjusted = ts - decomp.seasonal
else: # multiplicative
adjusted = ts / decomp.seasonal
return adjusted, decomp
# Apply seasonal adjustment
adjusted_ts, decomp = seasonal_adjustment(seasonal_ts)
# Plot original vs seasonally adjusted
fig, axes = plt.subplots(2, 1, figsize=(15, 10))
seasonal_ts.plot(ax=axes[0], title='Original Time Series', alpha=0.7)
decomp.trend.plot(ax=axes[0], color='red', linewidth=2, label='Trend')
axes[0].legend()
axes[0].set_ylabel('Value')
adjusted_ts.plot(ax=axes[1], title='Seasonally Adjusted Time Series', color='green')
decomp.trend.plot(ax=axes[1], color='red', linewidth=2, alpha=0.7, label='Trend')
axes[1].legend()
axes[1].set_ylabel('Adjusted Value')
for ax in axes:
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
5. Exponential Smoothing Methods
Simple Exponential Smoothing
# Create rainfall data (similar to London rainfall example)
np.random.seed(42)
rainfall_data = np.random.normal(25, 5, 100) # London average ~25 inches
rainfall_dates = pd.date_range('1813-01-01', periods=100, freq='A')
rainfall_ts = pd.Series(rainfall_data, index=rainfall_dates)
def fit_simple_exponential_smoothing(ts, alpha=None, forecast_steps=10):
"""
Fit simple exponential smoothing model
"""
from statsmodels.tsa.holtwinters import SimpleExpSmoothing
model = SimpleExpSmoothing(ts)
fitted_model = model.fit(smoothing_level=alpha, optimized=True if alpha is None else False)
# Generate forecasts
forecast = fitted_model.forecast(steps=forecast_steps)
# Create confidence intervals (approximate)
residuals = ts - fitted_model.fittedvalues
mse = np.mean(residuals**2)
std_error = np.sqrt(mse)
return fitted_model, forecast, std_error
# Fit the model
model, forecast, std_error = fit_simple_exponential_smoothing(rainfall_ts)
print(f"Optimal alpha: {model.params['smoothing_level']:.4f}")
print(f"MSE: {np.mean((rainfall_ts - model.fittedvalues)**2):.4f}")
# Plot results
fig, ax = plt.subplots(figsize=(15, 8))
# Plot original data
rainfall_ts.plot(ax=ax, label='Observed', alpha=0.7)
# Plot fitted values
model.fittedvalues.plot(ax=ax, label='Fitted', color='red', linewidth=2)
# Plot forecasts
forecast_dates = pd.date_range(start=rainfall_ts.index[-1] + pd.DateOffset(years=1),
periods=len(forecast), freq='A')
forecast_series = pd.Series(forecast, index=forecast_dates)
forecast_series.plot(ax=ax, label='Forecast', color='green', linewidth=2, linestyle='--')
# Add confidence intervals for forecasts
upper_ci = forecast + 1.96 * std_error * np.sqrt(np.arange(1, len(forecast)+1))
lower_ci = forecast - 1.96 * std_error * np.sqrt(np.arange(1, len(forecast)+1))
ax.fill_between(forecast_dates, lower_ci, upper_ci, alpha=0.3, color='green', label='95% CI')
ax.set_title('Simple Exponential Smoothing - London Rainfall')
ax.set_ylabel('Rainfall (inches)')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
Holt’s Linear Exponential Smoothing
# Create data with trend (skirts example)
skirt_data = [608, 617, 615, 617, 619, 619, 620, 628, 634, 635, 639, 640, 641]
skirt_dates = pd.date_range('1866-01-01', periods=len(skirt_data), freq='A')
skirt_ts = pd.Series(skirt_data, index=skirt_dates)
def fit_holt_smoothing(ts, alpha=None, beta=None, forecast_steps=10):
"""
Fit Holt's linear exponential smoothing
"""
from statsmodels.tsa.holtwinters import Holt
model = Holt(ts)
fitted_model = model.fit(smoothing_level=alpha, smoothing_trend=beta,
optimized=True if alpha is None else False)
forecast = fitted_model.forecast(steps=forecast_steps)
# Calculate prediction intervals
residuals = ts - fitted_model.fittedvalues
mse = np.mean(residuals**2)
std_error = np.sqrt(mse)
return fitted_model, forecast, std_error
# Fit Holt's model
holt_model, holt_forecast, holt_std = fit_holt_smoothing(skirt_ts)
print(f"Optimal alpha (level): {holt_model.params['smoothing_level']:.4f}")
print(f"Optimal beta (trend): {holt_model.params['smoothing_trend']:.4f}")
# Plot results
fig, ax = plt.subplots(figsize=(15, 8))
skirt_ts.plot(ax=ax, label='Observed', marker='o', linewidth=2)
holt_model.fittedvalues.plot(ax=ax, label='Fitted', color='red', linewidth=2)
# Forecasts
forecast_dates = pd.date_range(start=skirt_ts.index[-1] + pd.DateOffset(years=1),
periods=len(holt_forecast), freq='A')
forecast_series = pd.Series(holt_forecast, index=forecast_dates)
forecast_series.plot(ax=ax, label='Forecast', color='green', linewidth=2, linestyle='--', marker='s')
ax.set_title("Holt's Linear Exponential Smoothing - Skirt Hem Diameters")
ax.set_ylabel('Diameter (mm)')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
Holt-Winters Seasonal Exponential Smoothing
# Create seasonal data (souvenirs example)
np.random.seed(42)
months = pd.date_range('1987-01-01', '1993-12-01', freq='M')
base_trend = np.linspace(1000, 8000, len(months))
seasonal_pattern = 1000 * np.sin(2 * np.pi * np.arange(len(months)) / 12) + 500
noise = np.random.normal(0, 200, len(months))
souvenirs_ts = pd.Series(base_trend + seasonal_pattern + noise, index=months)
def fit_holt_winters(ts, seasonal='additive', seasonal_periods=12, forecast_steps=24):
"""
Fit Holt-Winters exponential smoothing
"""
model = ExponentialSmoothing(ts, trend='add', seasonal=seasonal,
seasonal_periods=seasonal_periods)
fitted_model = model.fit()
forecast = fitted_model.forecast(steps=forecast_steps)
return fitted_model, forecast
# Fit Holt-Winters model
hw_model, hw_forecast = fit_holt_winters(souvenirs_ts)
print("Holt-Winters Parameters:")
print(f"Alpha (level): {hw_model.params['smoothing_level']:.4f}")
print(f"Beta (trend): {hw_model.params['smoothing_trend']:.4f}")
print(f"Gamma (seasonal): {hw_model.params['smoothing_seasonal']:.4f}")
# Plot results
fig, ax = plt.subplots(figsize=(15, 8))
souvenirs_ts.plot(ax=ax, label='Observed', alpha=0.7)
hw_model.fittedvalues.plot(ax=ax, label='Fitted', color='red', linewidth=2)
# Forecasts
forecast_dates = pd.date_range(start=souvenirs_ts.index[-1] + pd.DateOffset(months=1),
periods=len(hw_forecast), freq='M')
forecast_series = pd.Series(hw_forecast, index=forecast_dates)
forecast_series.plot(ax=ax, label='Forecast', color='green', linewidth=2, linestyle='--')
ax.set_title('Holt-Winters Seasonal Exponential Smoothing')
ax.set_ylabel('Sales')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
6. Model Validation and Diagnostics
Residual Analysis
def analyze_residuals(model, ts, title="Residual Analysis"):
"""
Comprehensive residual analysis
"""
residuals = ts - model.fittedvalues
fig, axes = plt.subplots(2, 2, figsize=(15, 10))
# Time plot of residuals
residuals.plot(ax=axes[0,0], title='Residuals vs Time')
axes[0,0].axhline(y=0, color='red', linestyle='--')
axes[0,0].set_ylabel('Residuals')
axes[0,0].grid(True, alpha=0.3)
# Histogram of residuals
residuals.hist(ax=axes[0,1], bins=20, density=True, alpha=0.7)
axes[0,1].set_title('Distribution of Residuals')
axes[0,1].set_xlabel('Residuals')
axes[0,1].set_ylabel('Density')
# Add normal distribution overlay
mu, sigma = residuals.mean(), residuals.std()
x = np.linspace(residuals.min(), residuals.max(), 100)
axes[0,1].plot(x, stats.norm.pdf(x, mu, sigma), 'r-', linewidth=2, label='Normal')
axes[0,1].legend()
axes[0,1].grid(True, alpha=0.3)
# Q-Q plot
from scipy import stats
stats.probplot(residuals, dist="norm", plot=axes[1,0])
axes[1,0].set_title('Q-Q Plot')
axes[1,0].grid(True, alpha=0.3)
# ACF of residuals
plot_acf(residuals.dropna(), ax=axes[1,1], title='ACF of Residuals')
plt.suptitle(title, fontsize=16, fontweight='bold')
plt.tight_layout()
plt.show()
# Statistical tests
print("=== Residual Analysis ===")
print(f"Mean of residuals: {residuals.mean():.6f}")
print(f"Std of residuals: {residuals.std():.6f}")
# Ljung-Box test for autocorrelation
ljung_box = acorr_ljungbox(residuals.dropna(), lags=20, return_df=True)
print(f"\nLjung-Box test (first 20 lags):")
print(f"Test statistic: {ljung_box['lb_stat'].iloc[-1]:.4f}")
print(f"P-value: {ljung_box['lb_pvalue'].iloc[-1]:.4f}")
if ljung_box['lb_pvalue'].iloc[-1] > 0.05:
print("✓ No significant autocorrelation detected (p > 0.05)")
else:
print("✗ Significant autocorrelation detected (p < 0.05)")
# Shapiro-Wilk test for normality
shapiro_stat, shapiro_p = stats.shapiro(residuals.dropna())
print(f"\nShapiro-Wilk normality test:")
print(f"Test statistic: {shapiro_stat:.4f}")
print(f"P-value: {shapiro_p:.4f}")
if shapiro_p > 0.05:
print("✓ Residuals appear normally distributed (p > 0.05)")
else:
print("✗ Residuals do not appear normally distributed (p < 0.05)")
# Analyze residuals for the Holt-Winters model
analyze_residuals(hw_model, souvenirs_ts, "Holt-Winters Residual Analysis")
Model Comparison
def compare_models(ts, models_dict):
"""
Compare multiple forecasting models using various metrics
"""
results = {}
for name, model in models_dict.items():
fitted_values = model.fittedvalues
residuals = ts - fitted_values
# Calculate metrics
mae = np.mean(np.abs(residuals))
mse = np.mean(residuals**2)
rmse = np.sqrt(mse)
mape = np.mean(np.abs(residuals / ts)) * 100
# AIC and BIC (if available)
aic = getattr(model, 'aic', np.nan)
bic = getattr(model, 'bic', np.nan)
results[name] = {
'MAE': mae,
'MSE': mse,
'RMSE': rmse,
'MAPE': mape,
'AIC': aic,
'BIC': bic
}
# Create comparison DataFrame
comparison_df = pd.DataFrame(results).T
return comparison_df
# Fit multiple models for comparison
simple_exp = SimpleExpSmoothing(souvenirs_ts).fit()
holt_exp = Holt(souvenirs_ts).fit()
hw_add = ExponentialSmoothing(souvenirs_ts, trend='add', seasonal='add', seasonal_periods=12).fit()
hw_mul = ExponentialSmoothing(souvenirs_ts, trend='add', seasonal='mul', seasonal_periods=12).fit()
models = {
'Simple Exponential': simple_exp,
'Holt Linear': holt_exp,
'Holt-Winters Additive': hw_add,
'Holt-Winters Multiplicative': hw_mul
}
comparison = compare_models(souvenirs_ts, models)
print("Model Comparison:")
print("=" * 80)
print(comparison.round(4))
# Find best model by RMSE
best_model = comparison['RMSE'].idxmin()
print(f"\nBest model by RMSE: {best_model}")
7. Advanced Topics
Differencing for Non-Stationary Time Series
def check_stationarity(ts, title="Time Series"):
"""
Check stationarity using visual inspection and statistical tests
"""
from statsmodels.tsa.stattools import adfuller
fig, axes = plt.subplots(2, 2, figsize=(15, 10))
# Original series
ts.plot(ax=axes[0,0], title=f'{title} - Original')
axes[0,0].set_ylabel('Value')
# First difference
ts_diff = ts.diff().dropna()
ts_diff.plot(ax=axes[0,1], title=f'{title} - First Difference')
axes[0,1].set_ylabel('First Difference')
# ACF and PACF
plot_acf(ts.dropna(), ax=axes[1,0], title='ACF - Original')
plot_acf(ts_diff, ax=axes[1,1], title='ACF - First Difference')
for ax in axes.flat:
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# Augmented Dickey-Fuller test
print(f"Stationarity Test Results for {title}:")
print("=" * 50)
# Original series
adf_result = adfuller(ts.dropna())
print(f"Original Series:")
print(f"ADF Statistic: {adf_result[0]:.6f}")
print(f"P-value: {adf_result[1]:.6f}")
print(f"Critical Values: {adf_result[4]}")
if adf_result[1] <= 0.05:
print("✓ Series is stationary (p <= 0.05)")
else:
print("✗ Series is non-stationary (p > 0.05)")
# First difference
adf_diff = adfuller(ts_diff)
print(f"\nFirst Difference:")
print(f"ADF Statistic: {adf_diff[0]:.6f}")
print(f"P-value: {adf_diff[1]:.6f}")
if adf_diff[1] <= 0.05:
print("✓ First difference is stationary (p <= 0.05)")
else:
print("✗ First difference is non-stationary (p > 0.05)")
# Test stationarity on our seasonal data
check_stationarity(souvenirs_ts, "Souvenirs Sales")
Box-Jenkins Methodology (ARIMA)
def auto_arima_analysis(ts, max_p=5, max_d=2, max_q=5, seasonal=True, m=12):
"""
Automatic ARIMA model selection (simplified version)
"""
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.statespace.sarimax import SARIMAX
best_aic = np.inf
best_order = None
best_seasonal_order = None
best_model = None
# Make series stationary
ts_diff = ts.diff().dropna()
print("Searching for best ARIMA model...")
# Grid search for best parameters
for p in range(max_p + 1):
for d in range(max_d + 1):
for q in range(max_q + 1):
try:
if seasonal:
# SARIMA model
model = SARIMAX(ts, order=(p, d, q),
seasonal_order=(1, 1, 1, m))
fitted_model = model.fit(disp=False)
else:
# ARIMA model
model = ARIMA(ts, order=(p, d, q))
fitted_model = model.fit()
if fitted_model.aic < best_aic:
best_aic = fitted_model.aic
best_order = (p, d, q)
if seasonal:
best_seasonal_order = (1, 1, 1, m)
best_model = fitted_model
except:
continue
print(f"Best ARIMA order: {best_order}")
if seasonal:
print(f"Best seasonal order: {best_seasonal_order}")
print(f"Best AIC: {best_aic:.4f}")
return best_model, best_order
# Example ARIMA analysis
try:
best_arima, best_order = auto_arima_analysis(souvenirs_ts, max_p=2, max_d=1, max_q=2)
# Plot ARIMA results
fig, ax = plt.subplots(figsize=(15, 8))
souvenirs_ts.plot(ax=ax, label='Observed', alpha=0.7)
best_arima.fittedvalues.plot(ax=ax, label='ARIMA Fitted', color='red', linewidth=2)
# Generate forecasts
forecast_arima = best_arima.forecast(steps=12)
forecast_dates = pd.date_range(start=souvenirs_ts.index[-1] + pd.DateOffset(months=1),
periods=12, freq='M')
forecast_series = pd.Series(forecast_arima, index=forecast_dates)
forecast_series.plot(ax=ax, label='ARIMA Forecast', color='green', linewidth=2, linestyle='--')
ax.set_title(f'ARIMA{best_order} Model Results')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
except Exception as e:
print(f"ARIMA analysis failed: {e}")
print("This is normal - ARIMA fitting can be sensitive to data characteristics.")
8. Practical Examples and Case Studies
Case Study 1: Economic Data Analysis
# Simulate economic data (GDP, unemployment, etc.)
def create_economic_data():
"""
Create simulated economic time series data
"""
np.random.seed(123)
dates = pd.date_range('2000-01-01', '2023-12-01', freq='Q')
# GDP growth (with trend and cycles)
gdp_trend = 2.5 + 0.3 * np.sin(2 * np.pi * np.arange(len(dates)) / 20)
gdp_cycle = 1.5 * np.sin(2 * np.pi * np.arange(len(dates)) / 40)
gdp_noise = np.random.normal(0, 0.8, len(dates))
gdp = gdp_trend + gdp_cycle + gdp_noise
# Unemployment rate (inverse relationship with GDP)
unemployment = 6.5 - 0.5 * gdp + np.random.normal(0, 0.5, len(dates))
unemployment = np.clip(unemployment, 3.0, 12.0) # Keep realistic bounds
return pd.DataFrame({
'GDP_Growth': gdp,
'Unemployment': unemployment
}, index=dates)
economic_data = create_economic_data()
# Analyze the relationship
fig, axes = plt.subplots(2, 2, figsize=(15, 10))
# Time series plots
economic_data['GDP_Growth'].plot(ax=axes[0,0], title='GDP Growth Rate', color='blue')
axes[0,0].set_ylabel('GDP Growth (%)')
axes[0,0].grid(True, alpha=0.3)
economic_data['Unemployment'].plot(ax=axes[0,1], title='Unemployment Rate', color='red')
axes[0,1].set_ylabel('Unemployment (%)')
axes[0,1].grid(True, alpha=0.3)
# Scatter plot
axes[1,0].scatter(economic_data['GDP_Growth'], economic_data['Unemployment'], alpha=0.6)
axes[1,0].set_xlabel('GDP Growth (%)')
axes[1,0].set_ylabel('Unemployment (%)')
axes[1,0].set_title('GDP vs Unemployment (Okun\'s Law)')
axes[1,0].grid(True, alpha=0.3)
# Add trend line
z = np.polyfit(economic_data['GDP_Growth'], economic_data['Unemployment'], 1)
p = np.poly1d(z)
axes[1,0].plot(economic_data['GDP_Growth'], p(economic_data['GDP_Growth']), "r--", alpha=0.8)
# Cross-correlation
from scipy.signal import correlate
gdp_norm = (economic_data['GDP_Growth'] - economic_data['GDP_Growth'].mean()) / economic_data['GDP_Growth'].std()
unemp_norm = (economic_data['Unemployment'] - economic_data['Unemployment'].mean()) / economic_data['Unemployment'].std()
correlation = correlate(gdp_norm, unemp_norm, mode='full')
lags = np.arange(-len(gdp_norm)+1, len(gdp_norm))
axes[1,1].plot(lags, correlation)
axes[1,1].set_xlabel('Lag')
axes[1,1].set_ylabel('Cross-correlation')
axes[1,1].set_title('GDP-Unemployment Cross-correlation')
axes[1,1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# Calculate correlation coefficient
corr_coef = economic_data['GDP_Growth'].corr(economic_data['Unemployment'])
print(f"Correlation between GDP Growth and Unemployment: {corr_coef:.4f}")
Case Study 2: Sales Forecasting with Multiple Seasonality
def create_sales_data_with_multiple_seasonality():
"""
Create sales data with weekly and monthly patterns
"""
np.random.seed(456)
dates = pd.date_range('2020-01-01', '2023-12-31', freq='D')
# Base trend
trend = 1000 + 50 * np.arange(len(dates)) / 365.25
# Monthly seasonality (higher sales in December, lower in February)
monthly_season = 200 * np.sin(2 * np.pi * (dates.dayofyear - 60) / 365.25)
# Weekly seasonality (higher sales on weekends)
weekly_season = 100 * np.sin(2 * np.pi * dates.dayofweek / 7)
# Random noise
noise = np.random.normal(0, 50, len(dates))
sales = trend + monthly_season + weekly_season + noise
sales = np.maximum(sales, 0) # No negative sales
return pd.Series(sales, index=dates)
daily_sales = create_sales_data_with_multiple_seasonality()
# Aggregate to different frequencies
weekly_sales = daily_sales.resample('W').sum()
monthly_sales = daily_sales.resample('M').sum()
# Create comprehensive analysis
fig, axes = plt.subplots(3, 2, figsize=(18, 15))
# Daily sales
daily_sales.plot(ax=axes[0,0], title='Daily Sales', alpha=0.7)
axes[0,0].set_ylabel('Daily Sales')
# Monthly aggregation
monthly_sales.plot(ax=axes[0,1], title='Monthly Sales', marker='o')
axes[0,1].set_ylabel('Monthly Sales')
# Seasonal decomposition of monthly data
monthly_decomp = seasonal_decompose(monthly_sales, model='additive', period=12)
monthly_decomp.trend.plot(ax=axes[1,0], title='Monthly Trend', color='red')
axes[1,0].set_ylabel('Trend')
monthly_decomp.seasonal.plot(ax=axes[1,1], title='Monthly Seasonal Component', color='green')
axes[1,1].set_ylabel('Seasonal')
# Weekly analysis
weekly_avg_by_day = daily_sales.groupby(daily_sales.index.dayofweek).mean()
day_names = ['Mon', 'Tue', 'Wed', 'Thu', 'Fri', 'Sat', 'Sun']
axes[2,0].bar(day_names, weekly_avg_by_day.values)
axes[2,0].set_title('Average Sales by Day of Week')
axes[2,0].set_ylabel('Average Daily Sales')
axes[2,0].tick_params(axis='x', rotation=45)
# Monthly seasonal pattern
monthly_avg = daily_sales.groupby(daily_sales.index.month).mean()
month_names = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun',
'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
axes[2,1].bar(month_names, monthly_avg.values)
axes[2,1].set_title('Average Sales by Month')
axes[2,1].set_ylabel('Average Daily Sales')
axes[2,1].tick_params(axis='x', rotation=45)
for ax in axes.flat:
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# Forecast monthly sales using Holt-Winters
hw_sales_model = ExponentialSmoothing(monthly_sales, trend='add', seasonal='add', seasonal_periods=12).fit()
sales_forecast = hw_sales_model.forecast(steps=12)
# Plot forecast
fig, ax = plt.subplots(figsize=(15, 8))
monthly_sales.plot(ax=ax, label='Historical Sales')
hw_sales_model.fittedvalues.plot(ax=ax, label='Fitted Values', color='red')
forecast_dates = pd.date_range(start=monthly_sales.index[-1] + pd.DateOffset(months=1),
periods=12, freq='M')
pd.Series(sales_forecast, index=forecast_dates).plot(ax=ax, label='Forecast',
color='green', linestyle='--', marker='o')
ax.set_title('Monthly Sales Forecast using Holt-Winters')
ax.set_ylabel('Monthly Sales')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
print("Sales Forecasting Summary:")
print("=" * 40)
print(f"Historical average monthly sales: ${monthly_sales.mean():,.0f}")
print(f"Forecasted average monthly sales: ${sales_forecast.mean():,.0f}")
print(f"Growth rate implied: {(sales_forecast.mean()/monthly_sales.mean() - 1)*100:.1f}%")
9. Best Practices and Tips
Model Selection Guidelines
def model_selection_guide():
"""
Print guidelines for selecting appropriate time series models
"""
guide = """
TIME SERIES MODEL SELECTION GUIDE
================================
1. DATA CHARACTERISTICS:
- Constant level, no trend, no seasonality → Simple Exponential Smoothing
- Linear trend, no seasonality → Holt's Linear Exponential Smoothing
- Linear trend with seasonality → Holt-Winters Exponential Smoothing
- Non-linear trends → Consider transformation or ARIMA models
- Multiple seasonalities → Consider complex seasonal models or decomposition
2. FORECASTING HORIZON:
- Short-term (< 1 season) → Exponential smoothing methods work well
- Medium-term (1-2 seasons) → Holt-Winters or seasonal ARIMA
- Long-term (> 2 seasons) → Consider structural models or machine learning
3. DATA FREQUENCY:
- Annual data → Simple/Holt's exponential smoothing
- Quarterly data → Consider seasonal patterns (period=4)
- Monthly data → Strong seasonal patterns likely (period=12)
- Daily data → Multiple seasonalities possible (weekly, monthly, yearly)
4. MODEL VALIDATION:
- Always check residuals for patterns
- Use out-of-sample validation when possible
- Consider multiple models and ensemble approaches
- Monitor forecast accuracy over time
5. PRACTICAL CONSIDERATIONS:
- Interpretability vs. accuracy trade-off
- Computational complexity for real-time applications
- Robustness to outliers and structural breaks
- Ability to incorporate external variables
"""
print(guide)
model_selection_guide()
Performance Metrics and Evaluation
def comprehensive_forecast_evaluation(actual, forecast, model_name="Model"):
"""
Comprehensive evaluation of forecast performance
"""
# Align the series
actual = actual.dropna()
forecast = forecast[:len(actual)]
# Calculate metrics
errors = actual - forecast
abs_errors = np.abs(errors)
pct_errors = 100 * abs_errors / np.abs(actual)
metrics = {
'MAE': np.mean(abs_errors),
'MSE': np.mean(errors**2),
'RMSE': np.sqrt(np.mean(errors**2)),
'MAPE': np.mean(pct_errors),
'Bias': np.mean(errors),
'Max Error': np.max(abs_errors),
'Min Error': np.min(abs_errors),
'Std of Errors': np.std(errors)
}
print(f"Forecast Evaluation for {model_name}")
print("=" * 50)
for metric, value in metrics.items():
print(f"{metric:15s}: {value:10.4f}")
# Direction accuracy (for trends)
actual_direction = np.sign(actual.diff().dropna())
forecast_direction = np.sign(pd.Series(forecast).diff().dropna())
direction_accuracy = np.mean(actual_direction == forecast_direction[:len(actual_direction)]) * 100
print(f"{'Direction Acc':15s}: {direction_accuracy:10.1f}%")
# Visualization
fig, axes = plt.subplots(2, 2, figsize=(15, 10))
# Actual vs Forecast
axes[0,0].plot(actual.values, label='Actual', alpha=0.7)
axes[0,0].plot(forecast, label='Forecast', alpha=0.7)
axes[0,0].set_title('Actual vs Forecast')
axes[0,0].legend()
axes[0,0].grid(True, alpha=0.3)
# Residuals
axes[0,1].plot(errors, alpha=0.7)
axes[0,1].axhline(y=0, color='red', linestyle='--')
axes[0,1].set_title('Forecast Errors')
axes[0,1].grid(True, alpha=0.3)
# Error histogram
axes[1,0].hist(errors, bins=20, density=True, alpha=0.7)
axes[1,0].set_title('Distribution of Errors')
axes[1,0].grid(True, alpha=0.3)
# Scatter plot
axes[1,1].scatter(actual, forecast, alpha=0.6)
min_val = min(actual.min(), np.min(forecast))
max_val = max(actual.max(), np.max(forecast))
axes[1,1].plot([min_val, max_val], [min_val, max_val], 'r--', alpha=0.8)
axes[1,1].set_xlabel('Actual')
axes[1,1].set_ylabel('Forecast')
axes[1,1].set_title('Actual vs Forecast Scatter')
axes[1,1].grid(True, alpha=0.3)
plt.suptitle(f'Forecast Evaluation: {model_name}', fontsize=16)
plt.tight_layout()
plt.show()
return metrics
# Example evaluation using our Holt-Winters model
if 'hw_model' in locals():
metrics = comprehensive_forecast_evaluation(souvenirs_ts, hw_model.fittedvalues,
"Holt-Winters Model")
10. Conclusion and Further Reading
This tutorial has covered the essential concepts and techniques for time series analysis in Python, including:
- Data preparation and visualization: Creating and plotting time series data
- Decomposition: Separating trend, seasonal, and irregular components
- Smoothing techniques: Moving averages and exponential smoothing methods
- Forecasting models: Simple, Holt’s, and Holt-Winters exponential smoothing
- Model validation: Residual analysis and performance evaluation
- Advanced topics: Stationarity testing and ARIMA modeling
Key Python Libraries for Time Series
- pandas: Data manipulation and time series functionality
- statsmodels: Statistical modeling and time series analysis
- scikit-learn: Machine learning approaches to forecasting
- matplotlib/seaborn: Visualization
- scipy: Statistical functions and tests
Recommended Next Steps
- Advanced ARIMA modeling: Learn about seasonal ARIMA (SARIMA) and automatic model selection
- Machine learning approaches: Explore Random Forest, XGBoost, and neural networks for forecasting
- Multivariate time series: Vector autoregression (VAR) and cointegration analysis
- State space models: Kalman filters and dynamic linear models
- Real-time forecasting: Stream processing and online learning techniques
Additional Resources
- Books: “Forecasting: Principles and Practice” by Hyndman & Athanasopoulos
- Python packages:
pmdarima,sktime,darts,prophet - Online courses: Time series analysis courses on Coursera, edX
- Documentation: Statsmodels and scikit-learn documentation
Sample Code Repository Structure
time_series_analysis/
├── data/
│ ├── raw/
│ └── processed/
├── notebooks/
│ ├── 01_data_exploration.ipynb
│ ├── 02_decomposition.ipynb
│ ├── 03_forecasting.ipynb
│ └── 04_model_validation.ipynb
├── src/
│ ├── data_processing.py
│ ├── models.py
│ ├── evaluation.py
│ └── visualization.py
├── tests/
└── requirements.txt
This tutorial provides a solid foundation for time series analysis in Python. The key is to practice with real data and understand the underlying assumptions and limitations of each method. Remember that no single model works for all time series – the art lies in selecting and combining the right techniques for your specific problem.
This tutorial is a Python adaptation of “A Little Book of R for Time Series” by Avril Coghlan. All examples have been converted to use Python libraries while maintaining the pedagogical structure of the original work.
Read more in https://amzn.to/4nMy9l1
Leave a Reply