Python Code

# import libraries
import pandas as pd
from statsforecast import StatsForecast
from statsforecast.models import AutoARIMA
import matplotlib.pyplot as plt

# read data
data_raw = pd.read_csv("../posts/2024-10-02-ts-fundamentals-whats-a-time-series/example_ts_data.csv")

data_raw = (
    # select columns
    data_raw[["Country", "Product", "Date", "Revenue"]]
    # change data types
    .assign(
        Date = pd.to_datetime(data_raw["Date"]), 
        Revenue = pd.to_numeric(data_raw["Revenue"])
    )
)

# print the first few rows
print(data_raw.head())

# filter on specific series
us_ic_raw = data_raw[(data_raw["Country"] == "United States") & (data_raw["Product"] == "Ice Cream")]

# create unique id
us_ic_raw["unique_id"] = us_ic_raw["Country"] + "_" + us_ic_raw["Product"]

# convert date to datetime
us_ic_raw["Date"] = pd.to_datetime(us_ic_raw["Date"])

# plot the data
plt.figure(figsize=(10, 6))
plt.plot(us_ic_raw.index, us_ic_raw["Revenue"], label="Ice Cream Revenue")
plt.xlabel("Date")
plt.ylabel("Revenue")
plt.title("Ice Cream Revenue in United States")
plt.legend()

# get final data for forecasting
us_ic_clean = us_ic_raw[["unique_id", "Date", "Revenue"]].copy()

# set up models to train
sf = StatsForecast(
    models=[AutoARIMA(season_length=12)],
    freq='ME',
)

# fit the model and forecast for 12 months ahead
Y_hat_df = sf.forecast(df = us_ic_clean, 
                       time_col = "Date", 
                       target_col = "Revenue", 
                       id_col = "unique_id", 
                       h=12, 
                       level=[95], 
                       fitted=True)

print(Y_hat_df.head())

# convert date to be first of the month
Y_hat_df["Date"] = Y_hat_df["Date"].dt.to_period("M").dt.to_timestamp()

# concat both df together
future_fcst_df = pd.concat([us_ic_clean, Y_hat_df], axis=0)

# make date the index
future_fcst_df.set_index("Date", inplace=True)

print(future_fcst_df.tail())

# plot the future fcst data
plt.figure(figsize=(10, 6))

# plot the original revenue data as line and forecast as dotted line
plt.plot(future_fcst_df.index, future_fcst_df["Revenue"], label="Actual Revenue")
plt.plot(future_fcst_df.index, future_fcst_df["AutoARIMA"], label="ARIMA Forecast", linestyle='dotted')

# plot the prediction intervals as shaded areas
plt.fill_between(future_fcst_df.index, 
                 future_fcst_df["AutoARIMA-lo-95"], 
                 future_fcst_df["AutoARIMA-hi-95"], 
                 color='gray', alpha=0.2, label='95% Prediction Interval')

# chart formatting
plt.xlabel("Date")
plt.ylabel("Revenue")
plt.title("Forecast Results for US Ice Cream Revenue")
plt.legend()

# save the plot
# plt.savefig("chart1", dpi = 300, bbox_inches = "tight")

# get fitted values of the historical data
# Note: The fitted values are the predicted values for the training data
residual_values = sf.forecast_fitted_values()

print(residual_values.head())

# plot the historical fitted values
plt.figure(figsize=(10, 6))

# plot the original revenue data as line and forecast as dotted line
plt.plot(residual_values.index, residual_values["Revenue"], label="Actual Revenue")
plt.plot(residual_values.index, residual_values["AutoARIMA"], label="Forecasted Revenue", linestyle='dotted')

# plot the prediction intervals as shaded areas
plt.fill_between(residual_values.index, 
                 residual_values["AutoARIMA-lo-95"], 
                 residual_values["AutoARIMA-hi-95"], 
                 color='gray', alpha=0.2, label='95% Prediction Interval')

# chart formatting
plt.xlabel("Date")
plt.ylabel("Revenue")
plt.title("ARIMA Residuals for US Ice Cream Revenue")
plt.legend()

# save the plot
# plt.savefig("chart2", dpi = 300, bbox_inches = "tight")

# calculate residuals and plot the residuals directly as a line chart 
residuals = residual_values["Revenue"] - residual_values["AutoARIMA"]

plt.figure(figsize=(10, 6))
plt.plot(residuals.index, residuals, label="Residuals", color='orange')
plt.axhline(y=0, color='red', linestyle='--')
plt.xlabel("Date")
plt.ylabel("Residuals")
plt.title("Residuals Over Time")

# save the plot
# plt.savefig("chart3", dpi = 300, bbox_inches = "tight")

# create histogram of residuals
plt.figure(figsize=(10, 6))
plt.hist(residuals, bins=6, color='blue', alpha=0.7)
plt.axvline(x=0, color='red', linestyle='--')
plt.xlabel("Residuals")
plt.ylabel("Frequency")
plt.title("Histogram of Residuals")

# save the plot
# plt.savefig("chart4", dpi = 300, bbox_inches = "tight")

# create chart that puts the actual values on one side and the fitted values on the other as a scatter plot
plt.figure(figsize=(10, 6))
# plot the actual values
plt.scatter(residual_values["Revenue"], residual_values["AutoARIMA"], label="Fitted Values", alpha=0.5)
# plot the 45 degree line
plt.plot([residual_values["Revenue"].min(), residual_values["Revenue"].max()], 
         [residual_values["Revenue"].min(), residual_values["Revenue"].max()], 
         color='red', linestyle='--', label="45 Degree Line")
# chart formatting
plt.xlabel("Actual Revenue")
plt.ylabel("Fitted Revenue")
plt.title("Q-Q Plot of Actual vs Fitted Values")

# save the plot
# plt.savefig("chart5", dpi = 300, bbox_inches = "tight")