Modeling Variable Seasonal Features with the Fourier Transform | by Florin Andrei | Oct, 2023

Enhance time collection forecast efficiency with a way from sign processing

Modeling time collection knowledge and forecasting it are complicated matters. There are numerous strategies that may very well be used to enhance mannequin efficiency for a forecasting job. We are going to talk about right here a way which will enhance the way in which forecasting fashions take up and make the most of time options, and generalize from them. The principle focus would be the creation of the seasonal options that feed the time collection forecasting mannequin in coaching — there are straightforward good points to be made right here for those who embody the Fourier remodel within the characteristic creation course of.

This text assumes you might be acquainted with the essential facets of time collection forecasting — we won’t talk about this subject basically, however solely a refinement of 1 facet of it. For an introduction to time collection forecasting, see the Kaggle course on this subject — the approach mentioned right here builds on high of their lesson on seasonality.

Take into account the Quarterly Australian Portland Cement manufacturing dataset, displaying the entire quarterly manufacturing of Portland cement in Australia (in hundreds of thousands of tonnes) from 1956:Q1 to 2014:Q1.

df = pd.read_csv(‘Quarterly_Australian_Portland_Cement_production_776_10.csv’, usecols=[‘time’, ‘value’])# convert time from yr float to a correct datetime formatdf[‘time’] = df[‘time’].apply(lambda x: str(int(x)) + ‘-‘ + str(int(1 + 12 * (x % 1))).rjust(2, ‘0’))df[‘time’] = pd.to_datetime(df[‘time’])df = df.set_index(‘time’).to_period()df.rename(columns={‘worth’: ‘manufacturing’}, inplace=True) productiontime 1956Q1 0.4651956Q2 0.5321956Q3 0.5611956Q4 0.5701957Q1 0.529… …2013Q1 2.0492013Q2 2.5282013Q3 2.6372013Q4 2.5652014Q1 2.229

[233 rows x 1 columns]

That is time collection knowledge with a development, seasonal parts, and different attributes. The observations are made quarterly, spanning a number of a long time. Let’s check out the seasonal parts first, utilizing the periodogram plot (all code is included within the companion pocket book, linked on the finish).

The periodogram exhibits the ability density of spectral parts (seasonal parts). The strongest element is the one with a interval equal to 4 quarters, or 1 yr. This confirms the visible impression that the strongest up-and-down variations within the graph occur about 10 occasions per decade. There’s additionally a element with a interval of two quarters — that’s the identical seasonal phenomenon, and it merely means the 4-quarter periodicity just isn’t a easy sine wave, however has a extra complicated form. We are going to ignore parts with durations larger than 4, for simplicity.

Should you attempt to match a mannequin to this knowledge, maybe with a purpose to forecast the cement manufacturing for occasions past the tip of the dataset, it will be a good suggestion to seize this yearly seasonality in a characteristic or set of options, and supply these to the mannequin in coaching. Let’s see an instance.

Seasonal parts might be modeled in a lot of other ways. Typically, they’re represented as time dummies, or as sine-cosine pairs. These are artificial options with a interval equal to the seasonality you’re making an attempt to mannequin, or a submultiple of it.

Yearly time dummies:

seasonal_year = DeterministicProcess(index=df.index, fixed=False, seasonal=True).in_sample()print(seasonal_year) s(1,4) s(2,4) s(3,4) s(4,4)time 1956Q1 1.0 0.0 0.0 0.01956Q2 0.0 1.0 0.0 0.01956Q3 0.0 0.0 1.0 0.01956Q4 0.0 0.0 0.0 1.01957Q1 1.0 0.0 0.0 0.0… … … … …2013Q1 1.0 0.0 0.0 0.02013Q2 0.0 1.0 0.0 0.02013Q3 0.0 0.0 1.0 0.02013Q4 0.0 0.0 0.0 1.02014Q1 1.0 0.0 0.0 0.0

[233 rows x 4 columns]

Yearly sine-cosine pairs:

cfr = CalendarFourier(freq=’Y’, order=2)seasonal_year_trig = DeterministicProcess(index=df.index, seasonal=False, additional_terms=[cfr]).in_sample()with pd.option_context(‘show.max_columns’, None, ‘show.expand_frame_repr’, False):print(seasonal_year_trig) sin(1,freq=A-DEC) cos(1,freq=A-DEC) sin(2,freq=A-DEC) cos(2,freq=A-DEC)time 1956Q1 0.000000 1.000000 0.000000 1.0000001956Q2 0.999963 0.008583 0.017166 -0.9998531956Q3 0.017166 -0.999853 -0.034328 0.9994111956Q4 -0.999963 -0.008583 0.017166 -0.9998531957Q1 0.000000 1.000000 0.000000 1.000000… … … … …2013Q1 0.000000 1.000000 0.000000 1.0000002013Q2 0.999769 0.021516 0.043022 -0.9990742013Q3 0.025818 -0.999667 -0.051620 0.9986672013Q4 -0.999917 -0.012910 0.025818 -0.9996672014Q1 0.000000 1.000000 0.000000 1.000000

[233 rows x 4 columns]

Time dummies can seize very complicated waveforms of the seasonal phenomenon. On the flip aspect, if the interval is N, then you definately want N time dummies — so, if the interval could be very lengthy, you’ll need lots of dummy columns, which might not be fascinating. For non-penalized fashions, usually simply N-1 dummies are used — one is dropped to keep away from collinearity points (we’ll ignore that right here).

Sine/cosine pairs can mannequin arbitrarily very long time durations. Every pair will mannequin a pure sine waveform with some preliminary section. Because the waveform of the seasonal characteristic turns into extra complicated, you’ll need so as to add extra pairs (improve the order of the output from CalendarFourier).

(Facet observe: we use a sine/cosine pair for every interval we wish to mannequin. What we actually need is only one column of A*sin(2*pi*t/T + phi) the place A is the load assigned by the mannequin to the column, t is the time index of the collection, and T is the element interval. However the mannequin can not modify the preliminary section phi whereas becoming the info — the sine values are pre-computed. Fortunately, the formulation above is equal to: A1*sin(2*pi*t/T) + A2*cos(2*pi*t/T) and the mannequin solely wants to seek out the weights A1 and A2.)

TLDR: What these two strategies have in frequent is that they symbolize seasonality by way of a set of repetitive options, with values that cycle as usually as as soon as per the time interval being modeled (time dummies, and the bottom sine/cosine pair), or a number of occasions per interval (larger order sine/cosine). And every one in every of these options has values various inside a hard and fast interval: from 0 to 1, or from -1 to 1. These are all of the options we’ve to mannequin seasonality right here.

Let’s see what occurs after we match a linear mannequin on such seasonal options. However first, we have to add development options to the options assortment used to coach the mannequin.

Let’s create development options after which be a part of them with the seasonal_year time dummies generated above:

trend_order = 2trend_year = DeterministicProcess(index=df.index, fixed=True, order=trend_order).in_sample()X = trend_year.copy()X = X.be a part of(seasonal_year) const development trend_squared s(1,4) s(2,4) s(3,4) s(4,4)time 1956Q1 1.0 1.0 1.0 1.0 0.0 0.0 0.01956Q2 1.0 2.0 4.0 0.0 1.0 0.0 0.01956Q3 1.0 3.0 9.0 0.0 0.0 1.0 0.01956Q4 1.0 4.0 16.0 0.0 0.0 0.0 1.01957Q1 1.0 5.0 25.0 1.0 0.0 0.0 0.0… … … … … … … …2013Q1 1.0 229.0 52441.0 1.0 0.0 0.0 0.02013Q2 1.0 230.0 52900.0 0.0 1.0 0.0 0.02013Q3 1.0 231.0 53361.0 0.0 0.0 1.0 0.02013Q4 1.0 232.0 53824.0 0.0 0.0 0.0 1.02014Q1 1.0 233.0 54289.0 1.0 0.0 0.0 0.0

[233 rows x 7 columns]

That is the X dataframe (options) that will probably be used to coach/validate the mannequin. We’re modeling the info with quadratic development options, plus the 4 time dummies wanted for the yearly seasonality. The y dataframe (goal) will probably be simply the cement manufacturing numbers.

Let’s carve a validation set out of the info, containing the yr 2010 observations. We are going to match a linear mannequin on the X options proven above and the y goal represented by cement manufacturing (the take a look at portion), after which we’ll consider mannequin efficiency in validation. We may even do the entire above with a trend-only mannequin that can solely match the development options however ignores seasonality.

def do_forecast(X, index_train, index_test, trend_order):X_train = X.loc[index_train]X_test = X.loc[index_test]

y_train = df[‘production’].loc[index_train]y_test = df[‘production’].loc[index_test]

mannequin = LinearRegression(fit_intercept=False)_ = mannequin.match(X_train, y_train)y_fore = pd.Collection(mannequin.predict(X_test), index=index_test)y_past = pd.Collection(mannequin.predict(X_train), index=index_train)

trend_columns = X_train.columns.to_list()[0 : trend_order + 1]model_trend = LinearRegression(fit_intercept=False)_ = model_trend.match(X_train[trend_columns], y_train)y_trend_fore = pd.Collection(model_trend.predict(X_test[trend_columns]), index=index_test)y_trend_past = pd.Collection(model_trend.predict(X_train[trend_columns]), index=index_train)

RMSLE = mean_squared_log_error(y_test, y_fore, squared=False)print(f’RMSLE: {RMSLE}’)

ax = df.plot(**plot_params, title=’AUS Cement Manufacturing – Forecast’)ax = y_past.plot(coloration=’C0′, label=’Backcast’)ax = y_fore.plot(coloration=’C3′, label=’Forecast’)ax = y_trend_past.plot(ax=ax, coloration=’C0′, linewidth=3, alpha=0.333, label=’Pattern Previous’)ax = y_trend_fore.plot(ax=ax, coloration=’C3′, linewidth=3, alpha=0.333, label=’Pattern Future’)_ = ax.legend()

do_forecast(X, index_train, index_test, trend_order)

RMSLE: 0.03846449744356434

Blue is prepare, pink is validation. The total mannequin is the sharp, skinny line. The trend-only mannequin is the vast, fuzzy line.

This isn’t dangerous, however there may be one obvious concern: the mannequin has realized a constant-amplitude yearly seasonality. Though the yearly variation truly will increase in time, the mannequin can solely persist with fixed-amplitude variations. Clearly, it is because we gave the mannequin solely fixed-amplitude seasonal options, and there is nothing else within the options (goal lags, and many others) to assist it overcome this concern.

Let’s dig deeper into the sign (the info) to see if there’s something there that might assist with the fixed-amplitude concern.

The periodogram will spotlight all spectral parts within the sign (all seasonal parts within the knowledge), and can present an outline of their general energy, but it surely’s an mixture, a sum over the entire time interval. It says nothing about how the amplitude of every seasonal element could range in time throughout the dataset.

To seize that variation, it’s important to use the Fourier spectrogram as an alternative. It’s just like the periodogram, however carried out repeatedly over many time home windows throughout the whole knowledge set. Each periodogram and spectrogram can be found as strategies within the scipy library.

Let’s plot the spectrogram for the seasonal parts with durations of two and 4 quarters, talked about above. As all the time, the total code is within the companion pocket book linked on the finish.

spectrum = compute_spectrum(df[‘production’], 4, 0.1)plot_spectrogram(spectrum, figsize_x=10)

What this diagram exhibits is the amplitude, the “energy” of the 2-quarter and 4-quarter parts over time. They’re fairly weak initially, however grow to be very sturdy round 2010, which matches the variations you may see within the knowledge set plot on the high of this text.

What if, as an alternative of feeding the mannequin fixed-amplitude seasonal options, we modify the amplitude of those options in time, matching what the spectrogram tells us? Would the ultimate mannequin carry out higher in validation?

Let’s choose the 4-quarter seasonal element. We may match a linear mannequin (referred to as the envelope mannequin) on the order=2 development of this element, on the prepare knowledge (mannequin.match()), and lengthen that development into validation / testing (mannequin.predict()):

envelope_features = DeterministicProcess(index=X.index, fixed=True, order=2).in_sample()

spec4_train = compute_spectrum(df[‘production’].loc[index_train], max_period=4)spec4_train

spec4_model = LinearRegression()spec4_model.match(envelope_features.loc[spec4_train.index], spec4_train[‘4.0’])spec4_regress = pd.Collection(spec4_model.predict(envelope_features), index=X.index)

ax = spec4_train[‘4.0′].plot(label=’element envelope’, coloration=’grey’)spec4_regress.loc[spec4_train.index].plot(ax=ax, coloration=’C0′, label=’envelope regression: previous’)spec4_regress.loc[index_test].plot(ax=ax, coloration=’C3′, label=’envelope regression: future’)_ = ax.legend()

The blue line is the energy of the variation of the 4-quarter element within the prepare knowledge, fitted as a quadratic development (order=2). The pink line is similar factor, prolonged (predicted) over the validation knowledge.

We’ve got modeled the variation in time of the 4-quarter seasonal element. We will take the output from the envelope mannequin, and multiply it by the point dummies akin to the 4-quarter seasonal element:

spec4_regress = spec4_regress / spec4_regress.imply()

season_columns = [‘s(1,4)’, ‘s(2,4)’, ‘s(3,4)’, ‘s(4,4)’]for c in season_columns:X[c] = X[c] * spec4_regressprint(X)

const development trend_squared s(1,4) s(2,4) s(3,4) s(4,4)time 1956Q1 1.0 1.0 1.0 0.179989 0.000000 0.000000 0.0000001956Q2 1.0 2.0 4.0 0.000000 0.181109 0.000000 0.0000001956Q3 1.0 3.0 9.0 0.000000 0.000000 0.182306 0.0000001956Q4 1.0 4.0 16.0 0.000000 0.000000 0.000000 0.1835811957Q1 1.0 5.0 25.0 0.184932 0.000000 0.000000 0.000000… … … … … … … …2013Q1 1.0 229.0 52441.0 2.434701 0.000000 0.000000 0.0000002013Q2 1.0 230.0 52900.0 0.000000 2.453436 0.000000 0.0000002013Q3 1.0 231.0 53361.0 0.000000 0.000000 2.472249 0.0000002013Q4 1.0 232.0 53824.0 0.000000 0.000000 0.000000 2.4911392014Q1 1.0 233.0 54289.0 2.510106 0.000000 0.000000 0.000000

[233 rows x 7 columns]

The 4 time dummies usually are not both 0 or 1 anymore. They’ve been multiplied by the element envelope, and now their amplitude varies in time, identical to the envelope.

Let’s prepare the primary mannequin once more, now utilizing the modified time dummies.

do_forecast(X, index_train, index_test, trend_order)RMSLE: 0.02546321729737165

We’re not aiming for an ideal match right here, since that is only a toy mannequin, but it surely’s apparent the mannequin performs higher, it’s extra carefully following the 4-quarter variations within the goal. The efficiency metric has improved additionally, by 51%, which isn’t dangerous in any respect.

Enhancing mannequin efficiency is an enormous subject. The habits of any mannequin doesn’t rely on a single characteristic, or on a single approach. Nonetheless, for those who’re trying to extract all you may out of a given mannequin, it is best to most likely feed it significant options. Time dummies, or sine/cosine Fourier pairs, are extra significant once they replicate the variation in time of the seasonality they’re modeling.

Adjusting the seasonal element’s envelope by way of the Fourier remodel is especially efficient for linear fashions. Boosted timber don’t profit as a lot, however you may nonetheless see enhancements nearly it doesn’t matter what mannequin you utilize. Should you’re utilizing a hybrid mannequin, the place a linear mannequin handles deterministic options (calendar, and many others), whereas a boosted timber mannequin handles extra complicated options, it will be important that the linear mannequin “will get it proper”, subsequently leaving much less work to do for the opposite mannequin.

It is usually essential that the envelope fashions you utilize to regulate seasonal options are solely skilled on the prepare knowledge, and they don’t see any testing knowledge whereas in coaching, identical to every other mannequin. When you modify the envelope, the time dummies include info from the goal — they don’t seem to be purely deterministic options anymore, that may be computed forward of time over any arbitrary forecast horizon. So at this level info may leak from validation/testing again into coaching knowledge for those who’re not cautious.

The information set used on this article is obtainable right here below the Public Area (CC0) license:

The code used on this article might be discovered right here:

A pocket book submitted to the Retailer Gross sales — Time Collection Forecasting competitors on Kaggle, utilizing concepts described on this article:

A GitHub repo containing the event model of the pocket book submitted to Kaggle is right here:

All photographs and code used on this article are created by the writer.