Fpppy section 12.1 MSTL code example suggestion

[presmerats]

Hi,
I've been looking at the Section 12.1 of the Python version and the final forecast visualization in Figure 12.4 doesn't look quite right. It seems that the forecast is in another time scale, the five days of forecast data seem to be squashed inside an interval of 3.5 days.

If we compare to the same section on the R version of the book, the same visualization (now figure 12.3) shows the forecasted values placed in the correct dates with the same time scale as the original data (although the date range of the visualisation only shows the forecast of 3 days)

Here is a python code to solve the problem and create a visualization of the forecast that is in the same time scale as the data.

First I will post the way I loaded the data (so that you can check if this looks coherent to you).

#load the data 
bank_calls = pd.read_csv('../data/bank_calls.csv', parse_dates=["ds"])

# let's first reindex the data 
bank_calls_or = bank_calls.reset_index(drop=True, ).reset_index(drop=False)
bank_calls_or.drop('ds', axis=1,inplace=True)
bank_calls_or.rename(columns={
    "index":"ds",
    }, 
    inplace=True)

# remove periods with 0 calls
mask_y_not_0 = bank_calls_or["y"] != 0 
bank_calls_or = bank_calls_or.loc[mask_y_not_0,:].copy()

# every data point is 5min from the previous one, expressed in units of seconds:  x60x5 
# this would result in an Unix Epoch of 0 which corresponds to 1970-01-01
# we move the start of the date range to 2003-03-03 -> + 1046674800
bank_calls_or['ds'] = bank_calls_or['ds']*5*60 + 1046674800
bank_calls_or['ds'] = bank_calls_or['ds'].astype("datetime64[s]")

# ensure the number of calls is a float from which we can compute the sqrt later on
bank_calls_or["y"] = bank_calls_or["y"].astype("float64")

Then I proceed to run the code for the MSTL model and its decomposition visualization that is indicated in the book. Copying it here for reference:

sf = StatsForecast(models=[MSTL(season_length=[169, 169 * 5])], freq="5min")

sf = sf.fit(bank_calls_or.assign(y = lambda df: np.sqrt(df["y"])))
dcmp = sf.fitted_[0, 0].model_

fig, axes = plt.subplots(nrows=5, ncols=1, sharex=True)
sns.lineplot(data=dcmp, x=dcmp.index, y="data", ax=axes[0], color="black")
sns.lineplot(data=dcmp, x=dcmp.index, y="trend", ax=axes[1], color="black")
sns.lineplot(data=dcmp, x=dcmp.index, y="seasonal169", ax=axes[2], color="black")
sns.lineplot(data=dcmp, x=dcmp.index, y="seasonal845", ax=axes[3], color="black")
sns.lineplot(data=dcmp, x=dcmp.index, y="remainder", ax=axes[4], color="black")
fig.suptitle("STL decomposition")
fig.text(0.5, 0.94, "sqrt(Calls) = trend + seasonal169 + seasonal845 + remainder", 
         ha='center')
plt.xlabel("")
plt.show()

Then finally comes the code that forecasts 5 days of data and generates date times for them, that's the part that I suggest to be changed:

from matplotlib.ticker import MaxNLocator

# We predict 5 days of 169 datapoints each,
fc = sf.predict(h=169 * 5, level=[80, 95])
fc[fc.select_dtypes(include=np.number).columns] = fc.select_dtypes(include=np.number).transform(lambda x: x ** 2)


# We generate 7 days of datapoints at a 5 minute interval
df_dates = pd.DataFrame({
    'ds': pd.date_range(start="2003-10-27 07:00:00", periods = 7*24*60//5, freq="5min")
})

# we keep desired timestamps from the date range:
# - only weekdays 
# - from 7am to 9:05pm 
df_dates = df_dates[df_dates["ds"].dt.weekday <= 4]

mask_hour1 = df_dates["ds"].dt.hour >= 7
mask_hour2 = df_dates["ds"].dt.hour < 21
mask_desired1 = mask_hour1 & mask_hour2 

mask_9pm = df_dates["ds"].dt.hour >= 21 
mask_10pm = df_dates["ds"].dt.hour < 22 
mask_9_05 = df_dates["ds"].dt.minute <= 5
mask_desired2 = mask_9pm & mask_10pm & mask_9_05

mask_desired = mask_desired1 | mask_desired2
df_dates = df_dates.loc[mask_desired,:]


# make sure both indices are the same so that they can be joined together
df_dates.reset_index(drop=True, inplace=True)
fc.reset_index(drop=True, inplace=True)

# Join the datasets so now fc has a ds column with the proper timestamps
fc = df_dates.join(fc,how="left", rsuffix="_")
fc.drop(labels=['ds_',], axis=1, inplace=True)


# finally fill the gaps
fc = fill_gaps(fc, freq="5min")


# remove periods with 0 calls and fill gaps in the original dataframe
mask_y_not_0 = bank_calls["y"] != 0 
bank_calls = bank_calls.loc[mask_y_not_0,:].copy()
bank_calls = fill_gaps(bank_calls, freq="5min")

_ = plot_series(bank_calls, fc, level=[80, 95], 
            max_insample_length=169 * (14 + 5), 
            xlabel="DateTime [5m]", ylabel="Calls",
            title="Five-minute call volume to bank", 
            rm_legend=True,
            #palette="turbo"
            )

ax = _.axes[0]
ax.xaxis.set_major_locator(MaxNLocator(integer=True))
ax.legend(bbox_to_anchor=(1.04, 1), loc="upper left")
_

And this is the output obtained (sorry I couldn't get the ggplot color flavors):

I hope this helps to correct it. Thanks anyway

[mergenthaler]

Thanks @presmerats for noticing that issue and providing a posible solution. We will take a close look and come back with an answer soon. Sorry that we did not answer before and PLEASE keep raising issues or suggestions in case you find more.