Lab 4: Forecasting¶

Your fourth graded Python lab. You forecast monthly carbon dioxide (CO₂) levels measured at the Mauna Loa Observatory: turn weekly readings into a monthly series, test it for stationarity, hold back the most recent part of the record, fit a seasonal ARIMA (SARIMA) model to the rest, forecast the held-back months, measure the error with RMSE and MAE, answer questions about what you see, and submit one PDF on GAP.

Adapted from the DATA 300 GAP course, In-Class Module 1, Lab 4: Forecasting. Written in our own words; follow GAP for the official wording.

1. What and why¶

The data is the Mauna Loa CO₂ record that ships with the statsmodels package: weekly readings, in parts per million (ppm), from 1958 to 2001, with a few weeks missing.

2. Route A: Try JupyterLab (as on GAP)¶

Part 1: Open and name a notebook¶

Open Try JupyterLab, start a Python (Pyodide) notebook and rename it with the file name GAP gives (your last name, the course code and "AppliedLab4").

Part 2: Load the data and test for stationarity¶

Paste GAP's first code block into a cell and run it. In outline, it:

1. installs statsmodels, scikit-learn, seaborn and matplotlib into the browser session with %pip;
2. imports pandas, numpy, matplotlib, statsmodels, the adfuller test and scikit-learn's error functions, and hides warnings;
3. loads the CO₂ data from statsmodels, resamples the weekly readings to monthly averages with resample('ME'), and fills the gaps by carrying the previous month forward (ffill);
4. plots the monthly series;
5. runs the Augmented Dickey-Fuller (ADF) test on the whole monthly series and prints its statistic, its p-value and a one-line conclusion.

Part 3: Train/test split in time order¶

Paste and run GAP's second block. It keeps the first 80% of the months as the training set and the last 20% as the test set (no shuffling), prints how many months are in each, and plots the two parts in different colours.

Part 4: Fit SARIMA and measure accuracy¶

Paste and run GAP's third block. It fits a SARIMA model with the orders GAP gives, using the training data only; forecasts as many months as the test set holds; prints the RMSE and MAE of those forecasts against the actual test values, in ppm; and plots the end of the training data, the actual test data and the forecast together. The fit takes a while: see Known pitfalls.

Part 5: Answer the questions¶

Add a Markdown cell at the bottom and answer GAP's four questions in full sentences. They cover: the ADF p-value and why it leads to differencing (d = 1); why a time-ordered train/test split is more rigorous than fitting on all the data, and what risk it reduces; your MAE and whether that error is excellent, acceptable or poor given the level of the series; and whether the forecast follows both the long-run trend and the yearly up-and-down pattern. Read every value from your own output, give units (ppm), and name the output it comes from.

Part 6: Export and submit¶

Run everything from the top, then File → Print → Save as PDF with Background graphics ticked. See Submission below.

3. Route B: offline in Anaconda JupyterLab¶

Use this route if %pip install is slow in Try JupyterLab or the SARIMA fit freezes the browser tab, or if you prefer your own laptop. The CO₂ data ships inside statsmodels, so nothing is downloaded. Create a notebook inside the DATA300_GAP folder, named as in Route A. Paste the block into one cell, or split it at the # --- lines into three cells that match GAP's parts. The questions are the same as in Route A.

Two small differences from GAP's code, neither of which changes any number: the monthly resampling chooses 'ME' or 'M' to suit your pandas version, and the charts use colour-blind-safe colours.

# Route B for Lab 4: run offline in Anaconda JupyterLab.
# The CO2 data ships with statsmodels, so nothing is installed or downloaded.
import warnings

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
from statsmodels.tsa.stattools import adfuller
from sklearn.metrics import mean_squared_error, mean_absolute_error

warnings.filterwarnings("ignore")

# --- Part 2: load, resample to monthly, plot, ADF test
dataset = sm.datasets.co2.load_pandas()
df = dataset.data

# pandas 2.2 renamed the month-end alias from 'M' to 'ME'; older versions only know 'M'
major, minor = (int(x) for x in pd.__version__.split(".")[:2])
MONTH_END = "ME" if (major, minor) >= (2, 2) else "M"
print("pandas", pd.__version__, "- monthly alias:", MONTH_END)

df = df.resample(MONTH_END).mean().ffill()

plt.figure(figsize=(10, 4))
plt.plot(df.index, df["co2"], color="#009E73")
plt.title("Monthly CO2 levels (Mauna Loa Observatory)")
plt.ylabel("CO2 (ppm)")
plt.xlabel("Date")
plt.grid(True, linestyle="--", alpha=0.5)
plt.show()

print("--- Augmented Dickey-Fuller (ADF) test ---")
adf_result = adfuller(df["co2"])
print(f"ADF statistic: {adf_result[0]:.4f}")
print(f"p-value: {adf_result[1]:.4f}")
if adf_result[1] > 0.05:
    print("Conclusion: p-value > 0.05, so the series is NOT stationary.")
else:
    print("Conclusion: p-value <= 0.05, so the series IS stationary.")

# --- Part 3: sequential train/test split (first 80% train, last 20% test)
train_size = int(len(df) * 0.8)
train_data, test_data = df.iloc[:train_size], df.iloc[train_size:]

print(f"Total observations: {len(df)}")
print(f"Training horizon:   {len(train_data)} months")
print(f"Testing horizon:    {len(test_data)} months")

plt.figure(figsize=(10, 4))
plt.plot(train_data.index, train_data["co2"], label="Training data (80%)", color="#0072B2")
plt.plot(test_data.index, test_data["co2"], label="Test data (20%), held back", color="#E69F00")
plt.title("Sequential train/test split")
plt.ylabel("CO2 (ppm)")
plt.xlabel("Date")
plt.legend(loc="upper left")
plt.grid(True, linestyle="--", alpha=0.5)
plt.show()

# --- Part 4: fit SARIMA(1,1,1)(1,1,1,12) on the training data only, forecast the test period
print("Fitting the SARIMA model (this can take 10-60 seconds)...")
model = sm.tsa.statespace.SARIMAX(train_data["co2"], order=(1, 1, 1),
                                  seasonal_order=(1, 1, 1, 12))
sarima_result = model.fit(disp=False)

predictions = sarima_result.get_forecast(steps=len(test_data))
pred_mean = predictions.predicted_mean

rmse = np.sqrt(mean_squared_error(test_data["co2"], pred_mean))
mae = mean_absolute_error(test_data["co2"], pred_mean)
print("\n--- Out-of-sample accuracy ---")
print(f"Root mean squared error (RMSE): {rmse:.2f} ppm")
print(f"Mean absolute error (MAE):      {mae:.2f} ppm")

plt.figure(figsize=(12, 6))
plt.plot(train_data.index[-50:], train_data["co2"].iloc[-50:], label="Recent training data", color="#0072B2")
plt.plot(test_data.index, test_data["co2"], label="Actual test data", color="black")
plt.plot(test_data.index, pred_mean, label="SARIMA forecast", color="#D55E00", linestyle="--")
plt.title("SARIMA forecast against the held-back test period")
plt.ylabel("CO2 (ppm)")
plt.xlabel("Date")
plt.legend(loc="upper left")
plt.grid(True, linestyle="--", alpha=0.5)
plt.show()

Export with Kernel → Restart Kernel and Run All Cells, then File → Print → Save as PDF.

4. Checks you can make¶

These tell you whether your output is right without giving any answers away.

• The monthly series starts in 1958 and ends in 2001, and its plot shows a rising line with a regular yearly wiggle. The y-axis is in ppm.
• The ADF output prints two numbers (a statistic and a p-value) and a conclusion that matches the p-value and the 0.05 rule.
• Training months + testing months = total observations, and the training months are about four times the testing months.
• On the split chart, the training data comes first and the test data after, with no gap and no overlap.
• The forecast covers exactly the test period: the dashed line starts where the training data ends.
• RMSE and MAE are positive, in ppm, and RMSE is at least as large as MAE (it always is).
• Your MAE answer compares the error with the level of the series (for example as a percentage of a typical value), not just as a bare number.
• Your answer on differencing states the ADF test's null hypothesis and connects the p-value to d = 1.

5. Known pitfalls and fixes¶

6. Beyond the lab¶

Optional, not graded. Try these in the same notebook after you have answered GAP's questions:

1. Run the ADF test on df["co2"].diff().dropna() and on df["co2"].diff().diff(12).dropna(). How do the conclusions change, and which part of SARIMA's orders does each difference correspond to?
2. Compute a seasonal naive baseline for the test period (each month's forecast is the value 12 months before; you will need the last 12 training months to start) and its MAE. Does SARIMA beat it, and by how much?
3. Add the MAPE (np.mean(np.abs((actual - forecast) / actual)) * 100) to your error measures. Why is it easy to explain to a policymaker?
4. Draw the forecast's 95% interval with predictions.conf_int() and plt.fill_between. Does it widen across the test period? How many actual months fall outside it?
5. Fit a simpler SARIMA(0,1,1)(0,1,1,12) and compare its RMSE, MAE and fit time with GAP's model. Which would you recommend, and why?
6. Write two sentences for a policymaker that give your forecast's accuracy in words, with units and a comparison to the level of the series, and one assumption behind it.

7. Submission¶

1. Run Kernel → Restart Kernel and Run All Cells (Route B) or run every cell from the top (Route A), and check that every output and chart is visible.
2. Export with File → Print → Save as PDF, Background graphics ticked. If you wrote answers in a separate document instead, combine the PDFs offline as described in Combining PDFs without internet.
3. Open the PDF and check every part is present and readable.
4. Upload it to Applied Lab 4 on GAP and submit before Thursday 12 November, 23:59 China Standard Time (how to submit).