Forecast Power Demand and Size a Genset

Compare SARIMA, Prophet, and XGBoost to forecast LV demand and size a generator.

Introduction

30 Second Summary

Every building that depends on backup power faces a critical question: how big should the generator be? Oversize it and you waste money on idle capacity. Undersize it and you risk blackouts when demand peaks.

In this project, you will build a complete forecasting pipeline that predicts power demand on a residential low-voltage network, then translates those predictions into a data-driven generator set sizing recommendation.

What You'll Build

You run your final notebook cell and see a printed engineering report: predicted peak kW, converted kVA, safety margin applied, and the recommended standard genset size. Every number traces back to a model you trained and validated against nearly four years of real household power data.

By the end of this project, you'll have:

  • A rigorous time-series diagnostic report showing exactly where and why the power series is non-stationary, seasonal, and heteroscedastic, built from ADF/KPSS tests and STL decomposition.
  • A head-to-head model comparison of SARIMA, Prophet, and XGBoost with quantified RMSE, MAE, and MAPE on held-out test data, so you can see which model tracks peaks most accurately.
  • A genset sizing output that converts your best model's forecast into a kVA recommendation using IEC 60364 and ISO 8528-1 methodology, ready to hand to an engineer.
  • Secret Mission: Add temperature as an external regressor via SARIMAX and Prophet, and quantify whether weather data improves forecast accuracy over calendar-only features.

Are there any prerequisites?

You need a working Python environment with Jupyter Notebook installed. Familiarity with pandas and basic statistics (mean, standard deviation, correlation) will help you follow along comfortably.

Before We Start

Before we dive in, let's get clear on what you're about to build and why it matters.

Set Up and Load the Data

Every forecasting model starts with real data. Without actual power consumption measurements, there is nothing to diagnose, model, or predict.

In this step, you will download nearly 4 years of minute-level household power data from the UCI Machine Learning Repository, load it into a pandas DataFrame, resample it to hourly resolution, and split it into training and test sets ready for time-series modelling.

In this step, get ready to:
  • Download the UCI Individual Household Electric Power Consumption dataset and load it into pandas.
  • Resample to hourly resolution, handle missing values, and select the target variable.
  • Split into train (2007-2009) and test (2010) sets and confirm the shapes.
Download the UCI dataset

The UCI Individual Household Electric Power Consumption dataset contains 2,075,259 one-minute measurements from a single household in Sceaux, France, recorded between December 2006 and November 2010.

  • Confirm the file is in place by running this in a notebook cell:
import os
print(os.path.getsize('household_power_consumption.txt') / 1e6, 'MB')

You should see a number close to 127, confirming the full file is present.

Seeing FileNotFoundError?

Make sure household_power_consumption.txt is in the same directory as your notebook, not inside a subfolder from the zip extraction.

help me fix FileNotFoundError when loading the UCI power consumption dataset

Load the raw data into pandas

The file is semicolon-delimited with separate Date and Time columns. You will parse them into a single datetime index and coerce all measurement columns to numeric values.

  • Load the dataset and print its shape by running this code in a new cell:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# Load dataset
df = pd.read_csv(
    'household_power_consumption.txt',
    sep=';',
    parse_dates={'datetime': ['Date', 'Time']},
    infer_datetime_format=True,
    na_values=['?'],
    low_memory=False
)
df.set_index('datetime', inplace=True)
df = df.apply(pd.to_numeric, errors='coerce')

print('Raw shape:', df.shape)  # (2075259, 7)

What does this code do?

  • sep=';' tells pandas the file uses semicolons as delimiters, not commas.
  • parse_dates={'datetime': ['Date', 'Time']} merges the two separate date and time columns into one datetime index.
  • na_values=['?'] marks the ? characters in the raw data as missing values (the dataset uses ? instead of blanks for some gaps).
  • df.apply(pd.to_numeric, errors='coerce') ensures all 7 measurement columns are proper numeric types, converting any remaining non-numeric entries to NaN.

Before you run the cell, do you expect the shape to show 7 columns or 9? Think about what happened to the Date and Time columns.

  • Run the cell. You should see:
Raw shape: (2075259, 7)

The output confirms all 2,075,259 rows loaded successfully. The Date and Time columns were merged into the index, leaving 7 numeric measurement columns.

Seeing fewer rows or a parsing error?

If your shape shows fewer rows, the file may have been partially downloaded. Re-download and confirm the file is approximately 127 MB.

If you see a date-parsing error, make sure you are using parse_dates={'datetime': ['Date', 'Time']} with curly braces (a dictionary), not a list.

help me debug a parsing error when loading the UCI power dataset with pandas

Resample to hourly and prepare the target variable

Minute-level data is too granular for the forecasting models you will build. Resampling to hourly resolution reduces noise and makes model fitting computationally feasible.

  • Resample the data, handle missing values, and select the target by running this code:
# Resample to hourly
df_hourly = df.resample('H').mean()

# Handle missing values
df_hourly = df_hourly.ffill().interpolate(method='time')

# Select target
y = df_hourly['Global_active_power']

What does this code do?

  • resample('H').mean() groups all 60 minute-level readings into hourly bins and takes their average. This gives you one value per hour.
  • ffill() fills any remaining NaN gaps by carrying the last valid observation forward. This handles short gaps where the meter temporarily stopped.
  • interpolate(method='time') fills any gaps at the start of the series (where forward-fill has no prior value) using time-weighted linear interpolation.
  • Global_active_power is the household's total active power draw in kilowatts. This is the variable you will forecast.
  • Confirm the resampled series has no missing values by running:
print('Hourly shape:', df_hourly.shape)
print('Missing values:', y.isna().sum())

You should see approximately 34,000+ rows and 0 missing values.

Split into train and test sets

You will use three full years (2007-2009) for training and the final year (2010) for testing. This gives the models enough history to learn seasonal patterns while reserving a clean holdout period.

  • Create the train/test split and print the results by running:
# Train/test split
train = y['2007':'2009']
test = y['2010']

print('Train shape:', train.shape, '| From:', train.index[0], 'to', train.index[-1])
print('Test shape:', test.shape, '| From:', test.index[0], 'to', test.index[-1])

What does this code do?

  • y['2007':'2009'] uses pandas datetime index slicing to grab all hourly data points between January 2007 and December 2009.
  • y['2010'] grabs all of 2010 as the unseen test window. The models will never train on this data.

What do you expect the train shape to be? Three years of hourly data is roughly 3 x 365 x 24 = 26,280 hours.

  • Run the cell. You should see the train set spanning from early 2007 through end of 2009, and the test set covering 2010.

Seeing an empty test set?

If test.shape shows 0 rows, your dataset may not contain 2010 data. Confirm the raw date range by running print(df.index.min(), df.index.max()). The data should span from December 2006 to November 2010.

help me debug my train/test split for the UCI power consumption time series

✔️ Awesome, I've got everything!

Great. Make sure you have saved your notebook with all cells executed.

ⓧ I'd like to double check the full code

Here is the complete code for this step. Confirm each block is present in your notebook cells:

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# Load dataset
df = pd.read_csv(
    'household_power_consumption.txt',
    sep=';',
    parse_dates={'datetime': ['Date', 'Time']},
    infer_datetime_format=True,
    na_values=['?'],
    low_memory=False
)
df.set_index('datetime', inplace=True)
df = df.apply(pd.to_numeric, errors='coerce')

print('Raw shape:', df.shape)  # (2075259, 7)
# Resample to hourly
df_hourly = df.resample('H').mean()

# Handle missing values
df_hourly = df_hourly.ffill().interpolate(method='time')

# Select target
y = df_hourly['Global_active_power']

# Train/test split
train = y['2007':'2009']
test = y['2010']

print('Train shape:', train.shape, '| From:', train.index[0], 'to', train.index[-1])
print('Test shape:', test.shape, '| From:', test.index[0], 'to', test.index[-1])

Your data pipeline is ready. You have 2M+ raw measurements loaded, cleaned, resampled to hourly, and split into training and test sets. Next up, you will diagnose this series for stationarity, seasonality, and heteroscedasticity to understand its statistical properties before fitting any models.

Explore and Diagnose the Time Series

You have clean hourly power data split into training and test sets. Before fitting any forecasting model, you need to understand the statistical properties of the series.

Skipping diagnostics leads to misspecified models. This step reveals whether the series is stationary, what seasonal patterns dominate, and whether variance is constant over time.

In this step, get ready to:
  • Plot time series profiles to identify visual patterns.
  • Run formal stationarity tests and compute ACF/PACF.
  • Decompose the series and test for heteroscedasticity.
Visualize patterns in the training data

Visual inspection is the first diagnostic. You need to see the overall trend, typical daily usage shape, and weekly rhythm before running any statistical tests.

  • In a new notebook cell, import the diagnostic functions you will use throughout this step by running this code:
from statsmodels.tsa.stattools import adfuller, kpss
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.seasonal import STL
from statsmodels.stats.diagnostic import acorr_ljungbox

What do these imports provide?

  • adfuller and kpss are formal stationarity tests with opposing null hypotheses.
  • plot_acf and plot_pacf visualize autocorrelation at different lags.
  • STL performs seasonal-trend decomposition using LOESS.
  • acorr_ljungbox tests whether a sequence of autocorrelations is significantly different from zero.
  • Run the cell. No output is expected. This just makes the functions available for subsequent cells.
  • In a new cell, generate three exploratory plots by running this code:
# --- Exploratory Plots ---
fig, axes = plt.subplots(3, 1, figsize=(14, 10))
train.plot(ax=axes[0], title='Full Training Time Series')
train.groupby(train.index.hour).mean().plot(ax=axes[1], title='Average Daily Profile (by Hour)')
train.groupby(train.index.dayofweek).mean().plot(ax=axes[2], title='Average Weekly Profile (by Day of Week)')
plt.tight_layout()
plt.show()

What do these plots show?

  • The top plot shows the full training series from 2007 to 2009. You can eyeball any long-term trend or shifts in level.
  • The middle plot groups all hours across all days, averaging them. This reveals the typical daily consumption shape (low overnight, peaks in morning and evening).
  • The bottom plot groups by day of week (0 = Monday, 6 = Sunday). This reveals weekly patterns such as higher weekend usage.
  • Run the cell. You should see three vertically stacked plots showing clear daily and weekly cycles.

Plots not rendering?

  • Make sure your earlier cell that imports matplotlib.pyplot ran successfully. The plt object must be available in your notebook session.
  • If you see a blank output, try adding %matplotlib inline at the top of your notebook.

help me debug why my matplotlib plots aren't showing in Jupyter

Test for stationarity and identify seasonal lags

A time series must be stationary (constant mean and variance over time) for many forecasting models to work correctly. Two complementary tests give you confidence in the diagnosis.

Why two stationarity tests?

The ADF test has a null hypothesis of "the series has a unit root" (non-stationary). The KPSS test has a null hypothesis of "the series is stationary." Running both gives a confirmatory cross-check. If ADF rejects its null AND KPSS fails to reject its null, you have strong evidence of stationarity.

Before you run the tests, what do you expect? Given the clear seasonal cycles you just plotted, do you think the series will test as stationary or non-stationary?

  • Run the ADF and KPSS stationarity tests by adding this code in a new cell:
# --- Stationarity Tests ---
adf_result = adfuller(train.dropna())
print(f'ADF p-value: {adf_result[1]:.6f}')  # expect < 0.05

kpss_result = kpss(train.dropna(), regression='c', nlags='auto')
print(f'KPSS p-value: {kpss_result[1]:.6f}')  # expect > 0.05

What do these results mean?

  • adfuller returns a tuple where index [1] is the p-value. A p-value below 0.05 means you reject the null of a unit root. The series is stationary.
  • kpss also returns a tuple where index [1] is the p-value. A p-value above 0.05 means you fail to reject the null of stationarity. The series is stationary.
  • Both tests agreeing on stationarity gives you confidence to proceed without differencing.
  • Run the cell. You should see an ADF p-value less than 0.05 and a KPSS p-value greater than 0.05, confirming stationarity.

Seeing an InterpolationWarning from KPSS?

This warning means the test statistic falls outside the precomputed lookup table. The reported p-value is a boundary approximation. If the KPSS p-value prints as 0.100000, it means the actual p-value is at least 0.10, which still indicates you cannot reject stationarity.

help me interpret the KPSS InterpolationWarning I'm seeing

  • Now visualize autocorrelation to identify seasonal lags. Add this code in a new cell:
# --- ACF / PACF ---
fig, axes = plt.subplots(2, 1, figsize=(14, 6))
plot_acf(train.dropna(), lags=200, ax=axes[0], title='ACF')
plot_pacf(train.dropna(), lags=200, ax=axes[1], title='PACF')
plt.tight_layout()
plt.show()

What do ACF and PACF reveal?

The ACF (Autocorrelation Function) shows how correlated the series is with itself at different lag distances. The PACF (Partial Autocorrelation Function) isolates the direct correlation at each lag, removing the effect of intermediate lags.

Look for significant spikes at lag 24 (one day) and lag 168 (one week). These confirm the daily and weekly seasonal periods you spotted in the exploratory plots.

  • Run the cell. You should see significant ACF spikes repeating at multiples of 24, with a prominent spike at lag 168 confirming weekly seasonality.
Decompose the series and test for heteroscedasticity

Seasonal decomposition separates the series into trend, seasonal, and residual components. Once you isolate the residuals, you can test whether variance is constant or changes over time.

  • Apply STL decomposition with a daily period by running this code in a new cell:
# --- STL Decomposition ---
stl = STL(train.dropna(), period=24)
stl_result = stl.fit()
stl_result.plot()
plt.suptitle('STL Decomposition (period=24)')
plt.tight_layout()
plt.show()

What does STL decomposition show?

  • The period=24 parameter tells STL to look for a repeating pattern every 24 hours.
  • The trend component shows the slowly-moving average level of power consumption.
  • The seasonal component shows the repeating 24-hour cycle extracted from the data.
  • The residual component is what remains after removing trend and seasonality. Ideally it should look like random noise with constant spread.
  • Run the cell. You should see a four-panel plot (observed, trend, seasonal, residual) with a clear daily wave in the seasonal panel.

The residuals do not look like uniform noise. They appear to have periods of high variance and periods of low variance. This is called heteroscedasticity. You will confirm this formally with a statistical test.

  • Plot the squared residuals and run the Ljung-Box test on them by adding this code in a new cell:
# --- Heteroscedasticity Check ---
residuals = stl_result.resid
plt.figure(figsize=(12, 4))
plt.plot(residuals**2)
plt.title('Squared STL Residuals (Heteroscedasticity Check)')
plt.show()

lb_test = acorr_ljungbox(residuals**2, lags=[24, 48, 168], return_df=True)
print('Ljung-Box on squared residuals:')
print(lb_test)

What does this test prove?

Squaring the residuals converts them into a measure of variance at each time point. If variance were constant, the squared residuals would show no autocorrelation.

The Ljung-Box test checks whether the autocorrelations of the squared residuals are collectively zero. Very small p-values (near 0.0) at lags 24, 48, and 168 confirm that the variance changes systematically over time. This is heteroscedasticity.

  • Run the cell. You should see a time series plot of squared residuals with visible clustering of high values. The Ljung-Box table should show p-values extremely close to zero at all tested lags, confirming heteroscedasticity.

Getting a ValueError about NaN values?

The acorr_ljungbox function does not handle missing values. Make sure you passed train.dropna() into the STL constructor. If any NaN remains in the residuals, the test will fail.

help me fix NaN errors in my Ljung-Box test on STL residuals

✔️ Awesome, I've got everything!

Great. Confirm that all your cells ran without errors and that you saved your notebook.

ⓧ I'd like to double check the full code

Scroll through your notebook and confirm these cells are present in order:

  • Original data loading and train/test split cell (from Step 1).
  • Imports cell with adfuller, kpss, plot_acf, plot_pacf, STL, acorr_ljungbox.
  • Exploratory plots cell (three subplots: full series, daily profile, weekly profile).
  • Stationarity tests cell (ADF and KPSS prints).
  • ACF/PACF cell (two subplots with 200 lags).
  • STL decomposition cell (period=24).
  • Heteroscedasticity check cell (squared residuals plot and Ljung-Box table).

You now have a complete diagnostic picture of the power consumption series. It is stationary, strongly seasonal at 24h and 168h periods, and heteroscedastic. Next up, you will fit parametric forecasting models and see how well they handle these properties.

Build Parametric Baselines with ETS and SARIMA

In the previous step, you diagnosed the power consumption time series and confirmed it is stationary, has strong daily and weekly seasonality, and shows heteroscedastic residuals.

Now you need to answer a critical question: can classical parametric models capture these patterns on their own? Fitting them first establishes a measurable baseline and reveals exactly where pure autoregressive structure falls short.

In this step, get ready to:
  • Fit an ETS model and evaluate its 7-day forecast accuracy.
  • Fit a SARIMA model and evaluate its 7-day forecast accuracy.
  • Diagnose SARIMA residuals to identify where the model fails.
Add imports and define the test window

Before fitting any model, you need the forecasting classes and evaluation metrics. You also need a consistent 7-day test window that every model will be evaluated against.

  • In your Jupyter Notebook, add a new cell below your existing code and paste the following imports and test window setup:
from statsmodels.tsa.holtwinters import ExponentialSmoothing
from statsmodels.tsa.statespace.sarimax import SARIMAX
from sklearn.metrics import mean_squared_error, mean_absolute_error
import warnings
warnings.filterwarnings('ignore')

# 7-day test window (168 hours)
test_7d = test.iloc[:168]

What does this code do?

  • ExponentialSmoothing is the ETS (Error-Trend-Seasonal) model class from statsmodels.
  • SARIMAX is the seasonal ARIMA class that handles both pure SARIMA and models with exogenous regressors.
  • test_7d = test.iloc[:168] slices exactly 168 hours (7 days) from the test set. Every model predicts the same 168 values, making comparison fair.
  • Run the cell.

No output appears. The cell executes without errors, confirming the imports resolved correctly.

Fit an ETS model with additive seasonality

Exponential Smoothing (ETS) is one of the simplest parametric forecasting methods. It decomposes a time series into error, trend, and seasonal components and forecasts by exponentially weighting recent observations.

  • Add a new cell and paste the following ETS model code:
# --- ETS Model ---
ets_model = ExponentialSmoothing(
    train,
    trend='add',
    seasonal='add',
    seasonal_periods=24
).fit(optimized=True)

ets_forecast = ets_model.forecast(168)
ets_rmse = np.sqrt(mean_squared_error(test_7d, ets_forecast))
ets_mae = mean_absolute_error(test_7d, ets_forecast)
ets_metrics = {'RMSE': ets_rmse, 'MAE': ets_mae}
print(f'ETS RMSE: {ets_rmse:.4f} | MAE: {ets_mae:.4f}')

What does this code do?

  • trend='add' and seasonal='add' specify additive trend and seasonal components. Additive means the seasonal pattern has a constant amplitude regardless of the trend level.
  • seasonal_periods=24 tells the model to look for a repeating pattern every 24 hours (one daily cycle).
  • ets_model.forecast(168) generates predictions for the next 168 hourly timestamps, covering exactly 7 days.
  • RMSE penalizes large errors heavily. MAE gives the average absolute deviation. Together they quantify how far off the forecast is.
  • Run the cell.

You should see a line printing the ETS RMSE and MAE values. These are your first baseline numbers.

Cell taking a long time?

ETS fitting on ~26,000 hourly training points can take 30-60 seconds. If the cell hangs beyond 2 minutes, interrupt the kernel and re-run. Make sure train contains only the Global_active_power series, not the full DataFrame.

Help me debug my ETS cell

Fit SARIMA and diagnose its residuals

SARIMA (Seasonal AutoRegressive Integrated Moving Average) is a more flexible parametric model. It captures both short-term autocorrelation and seasonal patterns through explicit AR and MA terms at both regular and seasonal lags.

  • Add a new cell and paste the following SARIMA model code:
# --- SARIMA Model ---
sarima_model = SARIMAX(
    train,
    order=(1, 0, 1),
    seasonal_order=(1, 0, 1, 24),
    enforce_stationarity=False,
    enforce_invertibility=False
).fit(disp=False)

print(sarima_model.summary())

sarima_forecast = sarima_model.forecast(168)
sarima_rmse = np.sqrt(mean_squared_error(test_7d, sarima_forecast))
sarima_mae = mean_absolute_error(test_7d, sarima_forecast)
sarima_metrics = {'RMSE': sarima_rmse, 'MAE': sarima_mae}
print(f'SARIMA RMSE: {sarima_rmse:.4f} | MAE: {sarima_mae:.4f}')

What do these SARIMA parameters mean?

  • order=(1, 0, 1) means 1 autoregressive term, 0 differencing, 1 moving average term at the non-seasonal level.
  • seasonal_order=(1, 0, 1, 24) adds 1 seasonal AR term, 0 seasonal differencing, 1 seasonal MA term, with a seasonal period of 24 hours.
  • The seasonal period is set to 24 (daily). Setting it to 168 (weekly) would be computationally prohibitive for statsmodels. This is a known limitation you will observe in the residual diagnostics below.
  • Run the cell.

You should see the full SARIMA model summary table, including the AIC value and coefficient estimates. Below it, the RMSE and MAE values print.

Seeing a convergence warning?

SARIMA fitting on 26,000+ points may show a convergence warning. The enforce_stationarity=False and enforce_invertibility=False flags help avoid hard failures. If you see a warning, the model still produces forecasts. Check that the printed RMSE is a finite number.

Help me troubleshoot my SARIMA model

Now comes the critical diagnostic question. The SARIMA model has a seasonal period of 24 hours. Your ACF/PACF analysis from Step 2 showed significant autocorrelation at both lag 24 and lag 168.

Before you plot the residual diagnostics, think: do you expect the SARIMA residuals to be clean white noise, or do you expect leftover structure at the weekly lag?

  • Add a new cell and paste the following residual diagnostics code:
# --- SARIMA Residual Diagnostics ---
sarima_resid = sarima_model.resid

fig, axes = plt.subplots(1, 2, figsize=(14, 5))
import scipy.stats as stats
stats.probplot(sarima_resid.dropna(), dist='norm', plot=axes[0])
axes[0].set_title('SARIMA Residuals QQ-Plot')
plot_acf(sarima_resid.dropna(), lags=200, ax=axes[1], title='ACF of SARIMA Residuals')
plt.tight_layout()
plt.show()

What are you looking for in these plots?

  • The QQ-plot compares your residuals against a perfect normal distribution. Points falling on the diagonal line mean the residuals are normally distributed. Deviations at the tails mean the model produces occasional large errors it cannot explain.
  • The ACF of residuals should ideally show no significant spikes. Any spike means the model left predictable structure in the data that it failed to capture.
  • Run the cell.

You should see two plots side by side. The QQ-plot tails deviate heavily from the diagonal, confirming heavy-tailed residuals. The ACF plot shows significant spikes at lag 168.

This is the key finding. The SARIMA model with seasonal_order=(1, 0, 1, 24) captures daily seasonality but completely misses the weekly cycle at lag 168. Those spikes are predictable structure left on the table.

Why can't SARIMA handle weekly seasonality here?

Setting the seasonal period to 168 (weekly) in statsmodels is computationally prohibitive. The state-space representation grows with the seasonal period, making s=168 impractical for 26,000+ observations.

This is the exact gap that Prophet and SARIMAX with exogenous calendar features will fix in the next step. Prophet handles multiple seasonalities natively, and exogenous day-of-week features let SARIMAX encode weekly patterns without expanding the seasonal period.

✔️ Awesome, I've got everything!

Great. Make sure your notebook is saved. You now have ETS and SARIMA baselines with quantified error metrics and a clear diagnostic showing where SARIMA falls short.

ⓧ I'd like to double check the full code

Scroll through your notebook and confirm you have these cells in order:

  • Your original imports and data loading cell (from Steps 1-2).
  • Your exploratory plots and diagnostics cells (from Step 2).
  • The new imports cell with ExponentialSmoothing, SARIMAX, metrics, warnings, and the test_7d definition.
  • The ETS model cell printing RMSE and MAE.
  • The SARIMA model cell printing the summary table and RMSE/MAE.
  • The SARIMA residual diagnostics cell showing the QQ-plot and ACF plot.

You now have two parametric baselines with hard numbers and a clear diagnosis of where they fail. The SARIMA model cannot capture weekly seasonality with its daily seasonal period, leaving significant predictable structure in the residuals. Next up, you will fix this gap by bringing in Prophet and SARIMAX with exogenous calendar features.

Enrich with Exogenous Features Using Prophet and SARIMAX

In Step 3, you fitted a SARIMA model and saw its residual diagnostics reveal two clear failures: heavy tails in the QQ-plot and significant ACF spikes at lag 168.

Those lag-168 spikes tell you the model is missing weekly seasonality. SARIMA with s=24 can only model a daily cycle. To capture the weekly pattern, you need to feed that information in from outside. Prophet handles multiple seasonalities natively, and SARIMAX lets you inject calendar features as exogenous regressors.

In this step, get ready to:
  • Engineer cyclical calendar features (hour-of-day, day-of-week, is_weekend).
  • Fit Prophet with daily and weekly seasonality and evaluate its 7-day forecast.
  • Fit SARIMAX with calendar exogenous features and compare residual ACF to SARIMA.
Engineer calendar exogenous features

Raw hour and day-of-week integers create an artificial gap between the highest and lowest values (e.g., hour 23 and hour 0 are adjacent in time but 23 apart numerically). Cyclical sin/cos encoding solves this by placing adjacent hours close together on a unit circle.

  • First, add the Prophet import to your imports cell. Find this line at the top of your notebook:
from sklearn.metrics import mean_squared_error, mean_absolute_error
  • Add this import directly below it:
from prophet import Prophet
  • Run the imports cell to reload all dependencies.
  • Create a new cell below your existing SARIMA residual diagnostics. Build the calendar feature function by pasting this code:
# --- Calendar Feature Engineering ---
def make_calendar_features(index):
    """Create cyclical time features from a DatetimeIndex."""
    df_feat = pd.DataFrame(index=index)
    # Cyclical hour encoding (0-23 mapped to unit circle)
    df_feat['hour_sin'] = np.sin(2 * np.pi * index.hour / 24)
    df_feat['hour_cos'] = np.cos(2 * np.pi * index.hour / 24)
    # Cyclical day-of-week encoding (0-6 mapped to unit circle)
    df_feat['dow_sin'] = np.sin(2 * np.pi * index.dayofweek / 7)
    df_feat['dow_cos'] = np.cos(2 * np.pi * index.dayofweek / 7)
    # Binary weekend flag
    df_feat['is_weekend'] = (index.dayofweek >= 5).astype(int)
    return df_feat

train_exog = make_calendar_features(train.index)
test_7d_exog = make_calendar_features(test_7d.index)
print(train_exog.head(8))

What does this code do?

The make_calendar_features function takes a DatetimeIndex and returns a DataFrame with five columns. Each column gives the model a different angle on "when" a measurement was taken.

  • hour_sin and hour_cos encode the 24-hour cycle so midnight (hour 0) and 11 PM (hour 23) sit next to each other.
  • dow_sin and dow_cos do the same for the 7-day week, capturing the weekly pattern SARIMA missed.
  • is_weekend is a simple binary flag (1 for Saturday/Sunday, 0 otherwise) because residential power use often differs on weekends.
  • Run the cell.

You should see a table with 8 rows showing the five feature columns. Notice that hour_sin starts at 0.0 for midnight and hour_cos starts at 1.0, confirming the cyclical encoding is working.

Seeing a NameError for train or test_7d?

Make sure you ran all previous cells in order. The train and test_7d variables come from earlier cells. Restart the kernel and run all cells from the top if needed.

help me debug a NameError when running my calendar feature engineering cell

Fit Prophet with daily and weekly seasonality

Prophet expects a DataFrame with two columns: ds (datetime) and y (values). It decomposes the series into trend, seasonalities, and holidays internally. You enable daily and weekly seasonality so it captures both the 24-hour and 168-hour patterns your SARIMA model could not.

  • Create a new cell and fit the Prophet model by pasting this code:
# --- Prophet Model ---
# Reshape training data into Prophet's expected format
train_prophet = train.reset_index().rename(
    columns={'datetime': 'ds', 'Global_active_power': 'y'}
)

# Fit with explicit daily + weekly seasonality
prophet_model = Prophet(
    daily_seasonality=True,
    weekly_seasonality=True,
    yearly_seasonality=False
)
prophet_model.fit(train_prophet)

Why disable yearly seasonality?

Your test window is only 7 days. Yearly seasonality adds complexity without improving a short-horizon forecast. Keeping only daily and weekly seasonality focuses the model on the patterns that matter for this prediction window.

  • Run the cell. Prophet will print fitting status messages as it optimizes. This may take 10-30 seconds.

You should see output confirming the model fit completed with no convergence warnings.

  • Create a new cell to generate the 7-day forecast and compute metrics:
# Generate 168-hour (7-day) forecast
future = prophet_model.make_future_dataframe(periods=168, freq='H')
prophet_pred = prophet_model.predict(future)
prophet_forecast = prophet_pred.set_index('ds')['yhat'].loc[test_7d.index]

# Compute error metrics
prophet_rmse = np.sqrt(mean_squared_error(test_7d, prophet_forecast))
prophet_mae = mean_absolute_error(test_7d, prophet_forecast)
prophet_mape = np.mean(np.abs(
    (test_7d.values - prophet_forecast.values) / test_7d.values
)) * 100
prophet_metrics = {'RMSE': prophet_rmse, 'MAE': prophet_mae, 'MAPE': prophet_mape}
print(f'Prophet RMSE: {prophet_rmse:.4f} | MAE: {prophet_mae:.4f} | MAPE: {prophet_mape:.2f}%')

What does this code do?

  • make_future_dataframe(periods=168, freq='H') creates a DataFrame extending 168 hours beyond the training data (covering your 7-day test window).
  • prophet_pred.set_index('ds')['yhat'].loc[test_7d.index] extracts just the predicted values (yhat) aligned to your test timestamps.
  • MAPE (Mean Absolute Percentage Error) gives you a percentage-based accuracy measure that is easy to interpret across different scales.
  • Run the cell.

Compare the Prophet RMSE to the SARIMA RMSE from Step 3. Prophet should show a lower RMSE, confirming that handling multiple seasonalities natively improves accuracy.

Prophet forecast values are all NaN?

This happens when the test_7d.index timestamps do not align with the Prophet prediction timestamps. Check that your train series index has a proper DatetimeIndex named datetime. Run train.index.name to verify. If it shows None, re-run your data loading cells.

help me fix NaN values in my Prophet forecast alignment

Fit SARIMAX with calendar features and compare residuals

SARIMAX is the same SARIMA model but with an exog parameter that accepts your external feature matrix. By feeding in the calendar features you engineered, you give the model explicit information about weekly patterns that s=24 alone cannot capture.

  • Create a new cell and fit the SARIMAX model by pasting this code:
# --- SARIMAX Model ---
sarimax_model = SARIMAX(
    train,
    exog=train_exog,
    order=(1, 0, 1),
    seasonal_order=(1, 0, 1, 24),
    enforce_stationarity=False,
    enforce_invertibility=False
).fit(disp=False)

# Generate forecast with exogenous features for the test period
sarimax_forecast = sarimax_model.forecast(168, exog=test_7d_exog)
sarimax_rmse = np.sqrt(mean_squared_error(test_7d, sarimax_forecast))
sarimax_mae = mean_absolute_error(test_7d, sarimax_forecast)
sarimax_mape = np.mean(np.abs(
    (test_7d.values - sarimax_forecast.values) / test_7d.values
)) * 100
sarimax_metrics = {'RMSE': sarimax_rmse, 'MAE': sarimax_mae, 'MAPE': sarimax_mape}
print(f'SARIMAX RMSE: {sarimax_rmse:.4f} | MAE: {sarimax_mae:.4f} | MAPE: {sarimax_mape:.2f}%')

What does this code do?

  • exog=train_exog passes the 5-column calendar feature matrix alongside the time series during fitting.
  • sarimax_model.forecast(168, exog=test_7d_exog) requires you to pass the same features for the forecast period. The model needs to know what hour and day each future point falls on.
  • The model order (1,0,1)(1,0,1,24) is identical to the pure SARIMA. The only change is adding exogenous regressors, so any improvement is purely from the calendar information.
  • Run the cell. Fitting SARIMAX with exogenous variables takes longer than pure SARIMA (1-3 minutes depending on your machine).

You should see the SARIMAX metrics printed. Compare its RMSE to both SARIMA and Prophet.

Before you plot the residual comparison, what do you expect the SARIMAX residual ACF at lag 168 to look like compared to the pure SARIMA from Step 3?

  • Create a new cell and plot the side-by-side residual ACF comparison:
# --- Compare SARIMAX vs SARIMA Residual ACF ---
sarimax_resid = sarimax_model.resid
fig, axes = plt.subplots(1, 2, figsize=(14, 4))
plot_acf(sarima_model.resid.dropna(), lags=200, ax=axes[0], title='SARIMA Residuals ACF')
plot_acf(sarimax_resid.dropna(), lags=200, ax=axes[1], title='SARIMAX Residuals ACF')
plt.tight_layout()
plt.show()
  • Run the cell.

You should see two ACF plots side-by-side. The left plot (SARIMA) shows the significant spikes at lag 168 from Step 3. The right plot (SARIMAX) should show those lag-168 spikes reduced or completely within the confidence bands. This confirms that injecting weekly calendar features solved the weekly seasonality leakage.

SARIMAX residual ACF still shows lag-168 spikes?

Check that train_exog has the same length as train by running len(train_exog) == len(train). If they differ, re-run the calendar feature cell. Also verify that dow_sin and dow_cos columns are not all zeros by checking train_exog.describe().

help me figure out why SARIMAX residuals still show weekly seasonality

✔️ Awesome, I've got everything!

Great. Make sure all cells are saved. You now have Prophet and SARIMAX metrics ready for the final model comparison in the next step.

ⓧ I'd like to double check the full code

Your notebook should now contain these sections in order. Scroll through and confirm each one is present:

  • Imports cell (including from prophet import Prophet).
  • Data loading and resampling cell.
  • Diagnostics cells (ADF, KPSS, STL from Step 2).
  • ETS model cell with ets_metrics.
  • SARIMA model cell with sarima_metrics and residual diagnostics.
  • Calendar feature engineering cell (make_calendar_features function, train_exog, test_7d_exog).
  • Prophet model fitting cell.
  • Prophet forecast and metrics cell (prophet_metrics).
  • SARIMAX model and metrics cell (sarimax_metrics).
  • SARIMAX vs SARIMA residual ACF comparison cell.

You have now proven that adding exogenous calendar information fixes the weekly seasonality gap SARIMA left behind. Next, you will train an XGBoost model with even richer features and run a head-to-head comparison of all four models.

Forecast with XGBoost and Compare All Models

The parametric models from the previous steps revealed their limits. SARIMA could not capture weekly seasonality without prohibitive computation, and Prophet and SARIMAX improved by adding calendar features.

Tree-based ML models can capture non-linear interactions between features that parametric models cannot. XGBoost with engineered features is the proven top performer for tabular time-series forecasting. In this step, you train an XGBoost model and run a final head-to-head comparison of all five models.

In this step, get ready to:
  • Engineer lag, rolling, and calendar features for XGBoost.
  • Train an XGBRegressor and compute RMSE, MAE, and MAPE on the 7-day test set.
  • Build a final model comparison table and identify what drives peak demand prediction.
Engineer features for XGBoost

XGBoost does not understand time natively. You need to hand-engineer temporal features that encode what parametric models get implicitly from their structure: recent values (lags), trend context (rolling statistics), and cyclical calendar patterns.

  • Create a new cell in your notebook and add the XGBoost feature engineering code:
from xgboost import XGBRegressor

# --- XGBoost Feature Engineering ---
def make_xgb_features(series):
    df_feat = pd.DataFrame({'y': series})
    df_feat['lag_1'] = df_feat['y'].shift(1)
    df_feat['lag_24'] = df_feat['y'].shift(24)
    df_feat['lag_168'] = df_feat['y'].shift(168)
    df_feat['rolling_mean_24'] = df_feat['y'].shift(1).rolling(24).mean()
    df_feat['rolling_std_24'] = df_feat['y'].shift(1).rolling(24).std()
    df_feat['hour_sin'] = np.sin(2 * np.pi * series.index.hour / 24)
    df_feat['hour_cos'] = np.cos(2 * np.pi * series.index.hour / 24)
    df_feat['dow_sin'] = np.sin(2 * np.pi * series.index.dayofweek / 7)
    df_feat['dow_cos'] = np.cos(2 * np.pi * series.index.dayofweek / 7)
    df_feat['is_weekend'] = (series.index.dayofweek >= 5).astype(int)
    return df_feat.drop(columns=['y'])

full_series = pd.concat([train, test_7d])
full_features = make_xgb_features(full_series)

train_features = full_features.loc[train.index].dropna()
y_train = train.loc[train_features.index]

test_7d_features = full_features.loc[test_7d.index].dropna()
y_test_7d = test_7d.loc[test_7d_features.index]

What does this code do?

The make_xgb_features function builds 10 features from the raw time series:

  • Lag features (lag_1, lag_24, lag_168) give the model access to the value 1 hour, 1 day, and 1 week ago.
  • Rolling statistics (rolling_mean_24, rolling_std_24) summarize the last 24 hours of demand as context for trend and volatility.
  • Cyclical encodings (sin/cos of hour and day-of-week) let the model learn smooth periodic patterns without treating hours 23 and 0 as far apart.
  • The is_weekend flag captures the distinct consumption profile of weekends.

You concatenate train and test before computing features so lag values that span the train/test boundary are computed correctly, then split back into separate matrices.

  • Run the cell.
  • In a new cell, type train_features.shape and run it. You should see a tuple showing roughly 26,000 rows and 10 columns.

Seeing a NameError?

If train or test_7d is not defined, re-run the earlier cells that load the data and create the train/test split. Help me fix this NameError in my XGBoost feature cell.

Train XGBoost and compute metrics

With the feature matrices ready, you can train the model. XGBoost uses gradient-boosted decision trees that learn sequentially, with each tree correcting the previous one's errors.

  • Add a new cell and train the XGBRegressor:
# Train XGBoost
xgb_model = XGBRegressor(n_estimators=500, learning_rate=0.05, max_depth=6, random_state=42)
xgb_model.fit(train_features, y_train)

xgb_predictions = xgb_model.predict(test_7d_features)
xgb_rmse = np.sqrt(mean_squared_error(y_test_7d, xgb_predictions))
xgb_mae = mean_absolute_error(y_test_7d, xgb_predictions)
xgb_mape = np.mean(np.abs((y_test_7d.values - xgb_predictions) / y_test_7d.values)) * 100
xgb_metrics = {'RMSE': xgb_rmse, 'MAE': xgb_mae, 'MAPE': xgb_mape}

# Best model forecast
best_forecast = xgb_predictions

What does this code do?

  • The model trains with 500 boosting rounds, a learning rate of 0.05, and max tree depth of 6. These are reasonable defaults for hourly demand data.
  • After predicting on the 7-day test features, the code computes the same three error metrics (RMSE, MAE, MAPE) used for all prior models so you can compare fairly.
  • The best_forecast variable stores the XGBoost predictions for use in the generator sizing step later.
  • Run the cell. Training completes silently (XGBoost does not print progress by default without an evaluation set).
  • In a new cell, type xgb_metrics and run it. You should see a dictionary with RMSE, MAE, and MAPE values.

Getting a ValueError about feature shape?

Make sure you ran the feature engineering cell above first. The train_features and test_7d_features DataFrames must exist before training. Also confirm that dropna() ran successfully. Lag features create NaN rows at the start of the series that must be removed. Help me debug this XGBoost training error.

Compare all models and identify what drives peak demand

You now have metrics from five models. Time to put them side by side and see which one wins on the held-out test window.

  • Add a new cell to build the comparison table:
# --- Model Comparison Table ---
comparison = pd.DataFrame({
    'ETS': ets_metrics,
    'SARIMA': sarima_metrics,
    'Prophet': prophet_metrics,
    'SARIMAX': sarimax_metrics,
    'XGBoost': xgb_metrics
}).T
print('\n=== Model Comparison ===')
print(comparison.to_string())

What does this table show?

Each row is a model. Each column is an error metric (lower is better). MAPE is the most interpretable because it is percentage-based and scale-independent. XGBoost should achieve the lowest MAPE, confirming that non-linear feature interactions outperform pure autoregressive structure for this dataset.

  • Run the cell. You should see a formatted table with five rows (ETS, SARIMA, Prophet, SARIMAX, XGBoost) and three columns (RMSE, MAE, MAPE).
  • Confirm that XGBoost has the lowest MAPE value in the table.

Next, visualize which features XGBoost relied on most for its predictions.

  • Add a new cell for the feature importance plot:
# --- Feature Importance ---
feat_imp = pd.Series(xgb_model.feature_importances_, index=train_features.columns).sort_values(ascending=False)
feat_imp.plot(kind='bar', figsize=(10, 4), title='XGBoost Feature Importance')
plt.tight_layout()
plt.show()

What should the feature importance show?

You should see lag_24 and the hour encodings (hour_sin, hour_cos) dominating the top of the chart. This confirms that yesterday's same-hour demand and the time of day are the strongest predictors of current demand. The weekly lag (lag_168) also ranks highly, validating the weekly seasonality pattern you identified in the diagnostic step.

  • Run the cell. You should see a horizontal bar chart with features ranked by importance.

Before you run the final plot, think about what the XGBoost forecast line should look like compared to the actual demand.

Do you expect XGBoost to track the peaks more closely than the parametric models did, or will it smooth them out?

  • Add a new cell for the predicted vs actual visualization:
# --- Predicted vs Actual ---
plt.figure(figsize=(14, 5))
plt.plot(y_test_7d.values, label='Actual', color='black')
plt.plot(xgb_predictions, label='XGBoost Forecast', color='blue')
plt.title('XGBoost: Predicted vs Actual (7-Day Test Window)')
plt.legend()
plt.show()
  • Run the cell. You should see the blue XGBoost forecast line tracking the black actual demand line closely, including the peaks.

Why does XGBoost track peaks better?

Parametric models like SARIMA generate forecasts from a weighted combination of past values and errors. They tend to regress toward the mean during peak hours. XGBoost learns non-linear splits that can respond sharply when lag features signal a high-demand hour approaching, producing forecasts that match peak magnitudes more faithfully.

XGBoost line looks flat or wrong?

If the prediction line is flat or wildly off, check that test_7d_features has the correct date range by running test_7d_features.index in a new cell. The dates should start at January 1, 2010. Also verify xgb_predictions has 168 values by checking its length. Help me debug my XGBoost forecast plot.

✔️ Awesome, I've got everything!

Great. Make sure you have saved your notebook with all five cells from this step.

ⓧ I'd like to double check the full code

Scroll through your notebook and confirm these cells are present in order:

  • XGBoost Feature Engineering cell: defines make_xgb_features and builds train_features, y_train, test_7d_features, y_test_7d.
  • Train XGBoost cell: creates xgb_model, computes xgb_metrics, and assigns best_forecast.
  • Model Comparison Table cell: prints the 5-model comparison DataFrame.
  • Feature Importance cell: bar chart of XGBoost feature importances.
  • Predicted vs Actual cell: line plot comparing XGBoost forecast to actual demand.

You have a winning model and you know exactly which features drive its accuracy. Next up, you translate these demand predictions into an engineering decision: sizing a generator set using IEC/ISO methodology.

Secret mission

Add Weather Regressors and Size the Generator Set

Your models use calendar features alone. Real power demand correlates strongly with outdoor temperature. Can adding weather data improve your forecasts? And what generator does the winning prediction actually call for? In this extension, you inject temperature as an external regressor, measure the accuracy impact, then translate your best forecast into an engineering-grade genset sizing recommendation.

Clean Up Your Resources

Clean Up Your Resources

Decide whether to keep your resources running, pause them to come back later, or delete them entirely. This project runs entirely locally with no paid services, so there are no ongoing costs.

Resources you used:

  • notebook.ipynb - your Jupyter Notebook containing the full forecasting analysis and genset sizing recommendation.
  • household_power_consumption.txt - the UCI dataset file (~127 MB) downloaded in Step 1.

Keep everything running

No action needed. Choose this if you plan to continue experimenting with additional models, feature engineering approaches, or different forecast horizons.

  • Your notebook retains the complete analysis pipeline, model comparisons, and genset sizing output.
  • The dataset remains available for resampling at different intervals or exploring sub-metering columns.

Pause - I'll come back to this later

Shut down your Jupyter Notebook server to free up memory. Your files stay on disk for when you return.

  • In the terminal where Jupyter is running, press Ctrl+C and confirm the shutdown.
  • Both notebook.ipynb and household_power_consumption.txt remain on your Desktop, ready to resume.

Delete - I don't want to use this again

Remove all project files to reclaim disk space. The dataset alone is ~127 MB.

  • In the terminal where Jupyter is running, press Ctrl+C and confirm the shutdown.
  • Delete the project folder by running this command on macOS/Linux:
rm -rf ~/Desktop/power-forecasting

On Windows, run:

rmdir /s /q %USERPROFILE%\Desktop\power-forecasting
  • Confirm the folder is gone by checking your Desktop. Both notebook.ipynb and household_power_consumption.txt are removed.

Nice Work!

Nice Work!

You just built a complete demand forecasting pipeline that bridges data science with power systems engineering. From raw minute-level household data to a justified generator set recommendation, every step was grounded in real statistical diagnostics and model comparison.

You've learned how to:

  • Diagnose a real-world LV power consumption time series for stationarity, seasonality, and heteroscedasticity using ADF/KPSS tests and STL decomposition.
  • Build and compare four forecasting models (ETS, SARIMA, Prophet, XGBoost) with rigorous train/test evaluation and quantified error metrics.
  • Translate ML demand predictions into a standard genset sizing recommendation using IEC 60364 and ISO 8528-1 methodology, including power factor conversion, safety margin, and load factor verification.
  • Secret Mission: Quantified the value of weather exogenous variables on forecast accuracy by adding temperature as an external regressor to both SARIMAX and Prophet.

Ready to quiz yourself?