Forecasting with Fourier series: part II

Unlocking 1000x speed with better accuracy

In the first post, we explored how to create a forecaster using Fourier series in a similar vain to Prophet by Meta. The resulting forecaster even beat Prophet in our (non-rigorous) test.

In this post, we will make some corrections and improvements that hopefully yield even better results.

It was just a phase

To recap: any periodic signal can be approximated by summing different sine and cosine waves - even a square wave as can be seen in this example:

Example from the Wikipedia page. As the number of terms increases, the fit improves.

We initiated a harmonic series of sine and cosine waves for each seasonality term (weekly, monthly, quarterly, yearly etc.) up to the Fourier order specified for that seasonality. E.g. if the Fourier order for weekly seasonality was set to 3, then we would generate a sine and cosine wave corresponding to the following frequencies: 1/7, 2/7, 3/7 — resulting in 6 waves in total (3 sine, 3 cosine) for the weekly seasonality.

In the first post, we then used gradient descent to fit both the amplitude and phase of our sine waves. Here lied the issue: whilst we were fitting these very well, it is actually totally unnecessary to fit the phase — a complete oversight on my part. This is because the Fourier series expansion will automatically converge on the target signal as more harmonics are added.

This meant that the gradient descent algorithm would take slightly longer to converge on the target signal as it must also explore deviations in phase rather than amplitude alone (as well as bias, trend and exogenous regressors). This also risks failure to converge or being stuck in some local optimum as if the phase is shifted too much in one direction it may end up back to where it started (as it is periodically repeating) or may overlap too much with its corresponding sine/cosine wave (as one is just shifted a quarter cycle of the other).

Removing gradient descent

Only having to fit the amplitude — or, in other words, weight — of each wave component means we can resolve our fit with a linear regression closed-form solution. This allows us to find the optimal solution more precisely, more eloquently and at barely a fraction of the original time it took to fit (e.g. milliseconds rather than seconds).

• The derivation of the closed-form solution is discussed in both the following: Efficient leave-one-out cross-validation and Linear regression with multiple targets.
• But to summarize, we can replace the entire gradient descent with the following:

np.linalg.inv(X.T @ X) @ X.T @ y

Where X is a matrix of the wave components, bias, trend and any exogenous regressors and y is the timeseries we are using to fit.

Updated results

To see if this has resulted in any improvements, we will use the same example as the previous post: train and bus ridership numbers from the Chigago Transport Authority. We will predict total ridership for one year at a time between 2003 to 2018 training on the previous two years of data for a total of 16 forecasts. The average MAPE score across the 16 forecasts will be taken for each forecast for the model evaluation.

In the original post: Prophet scored an average MAPE of 7.9% and our model achieved 7.4%.

FourierForecast scores: [0.06757175 0.06847599 0.07174342 0.06946523 0.07235714 0.08831864
 0.09115687 0.06550169 0.08012167 0.06707774 0.08542078 0.07360266
 0.07506197 0.07625026 0.07252445 0.06323971]
Prophet scores: [0.06361311 0.08916864 0.07080911 0.0656165  0.0788751  0.11362799
 0.06367285 0.06227564 0.13398795 0.07868203 0.06452716 0.07208917
 0.08931415 0.08308844 0.07541953 0.05568152]

FourierForecast average MAPE: 7.4%
Prophet average MAPE: 7.9%

The improved version (GitHub updated accordingly) scores 7.1%, which is a marginal improvement but coupled with ~1000x speed up and generally better practise from a methodological standpoint, this can be seen as a significant win over the first iteration!

[0.06462474 0.06586173 0.0673727  0.06767039 0.06862962 0.08440217
 0.09051008 0.06334214 0.07686107 0.06499376 0.08170539 0.07122229
 0.07231773 0.07286749 0.06890983 0.05964524]

FourierForecast version II average MAPE: 0.071