Fourier Terms in R: Model Complex Seasonality

Fourier terms are pairs of sine and cosine waves that stand in for a repeating seasonal pattern using only a handful of numbers. In R you build them with fourier() and hand them to auto.arima() as extra predictors, which lets a single model capture very long seasonal cycles, and even several overlapping cycles at once, that a seasonal ARIMA cannot. This tutorial uses the forecast and ggplot2 packages for the runnable code, and shows the modern fable version to run in your own R session.

Why do weekly and daily seasonal cycles break a seasonal ARIMA model?

Sales climb every December. Website traffic peaks every Monday. A seasonal ARIMA model copes fine with a short cycle like 12 months, but it starts to struggle when the cycle is long, say 52 weeks, and it cannot handle two cycles at the same time. The quickest way to feel the problem is to build a weekly series and let R try to fit it.

Let's create five years of weekly demand with a clear yearly rhythm baked in, then look at it.

The printout confirms what we built: 260 weekly observations with a season that is 52 steps long. The plot shows the same wave rolling past five times, once per year. This is exactly the kind of series a seasonal model should love, so let's hand it to auto.arima() and see what it does.

Read the label ARIMA(0,0,0)(1,1,0)[52]. The 1 in the middle of the seasonal bracket means R was forced to take one seasonal difference, subtracting each week from the same week a year earlier. That throws away the first 52 weeks before the model can learn anything. Worse, the entire yearly shape is now squeezed into a single seasonal coefficient, sar1. That is a rigid, blunt way to describe a rich pattern, and it only works at all because there is exactly one seasonal period. The moment you have a daily cycle sitting on top of a weekly one, this whole approach falls apart.

What exactly is a Fourier term?

A Fourier term is a single sine or cosine wave that completes a whole number of cycles per season. The first harmonic completes a single cycle per season, and each higher harmonic completes one more cycle than the last, so harmonic number k makes k cycles. The remarkable fact, proved by the mathematician Joseph Fourier, is that adding a few of these waves together can trace almost any repeating shape. We use them in sine-and-cosine pairs so the combined wave can peak at any point in the cycle, not just where a plain sine would.

Here is the same idea written down. For a season of length $m$, the $k$-th Fourier pair is

$$s_k(t) = \sin\left(\frac{2 \pi k t}{m}\right), \qquad c_k(t) = \cos\left(\frac{2 \pi k t}{m}\right)$$

and the seasonal part of the series is approximated by adding $K$ of these pairs on top of a baseline:

$$\hat{y}_t = a + \sum_{k=1}^{K}\left[\alpha_k \sin\left(\frac{2 \pi k t}{m}\right) + \beta_k \cos\left(\frac{2 \pi k t}{m}\right)\right] + N_t$$

Where:

• $t$ = the time index (week 1, 2, 3, and onward)
• $m$ = the season length (52 for weekly data)
• $K$ = how many sine and cosine pairs you keep
• $\alpha_k$ and $\beta_k$ = the wave sizes the model estimates
• $N_t$ = whatever is left over, handled by an ARIMA model

If formulas are not your thing, skip straight to the code, because the fourier() function builds all of this for you. Pass it the series and a value of K, and it returns a matrix of the sine and cosine columns.

With K = 2 you get four columns, one sine and one cosine for each of the first two harmonics. Each row is one week, and each number is where that wave sits at that moment. These four columns are the predictors you will feed to the model. It helps to actually see the waves they represent.

Over two years, S1 and C1 each make two full swings (one per year), while S2 makes four (two per year). The model will scale and add waves like these until their sum matches the real seasonal shape.

How do you fit a dynamic harmonic regression in R?

Now we combine the pieces. You feed the Fourier columns to auto.arima() through its xreg argument (short for external regressors) and set seasonal = FALSE. The Fourier terms carry the season, so auto.arima() is free to model only the short-term wobble that is left, treating it as ordinary, non-seasonal ARIMA errors. This two-part recipe, a regression on Fourier terms with ARIMA errors, is called dynamic harmonic regression.

The label now reads "Regression with ARIMA(0,1,1) errors", with no seasonal bracket at all. The season has moved into the six Fourier coefficients, from S1-52 through C3-52, and the only ARIMA structure left is a tiny non-seasonal (0,1,1). Notice too that this model uses every observation, with no year of data sacrificed to seasonal differencing.

Figure 1: The dynamic harmonic regression workflow: turn the season into Fourier columns, fit ARIMA on the rest, then forecast with future Fourier terms.

How do you choose K, the number of Fourier pairs?

K is the one dial you control, and it sets how detailed the seasonal curve can be. A small K gives a smooth, simple season. A large K gives a wiggly one that can start chasing noise. The clean way to pick it is to fit a handful of values and keep the one with the lowest AICc, a score that rewards good fit and penalises extra parameters.

The AICc drops sharply, hits its lowest point at K = 3, then creeps back up. That dip is the sweet spot: enough waves to trace the yearly shape, not so many that the model fits the random noise. A picture makes the smoothness trade-off obvious.

The K = 1 line is a single smooth swell, too plain to follow the real peaks and dips. The K = 3 line hugs the actual series closely. Push K much higher and that line would start wobbling around every little bump, which is overfitting.

How do you forecast future values with Fourier terms?

Forecasting has one twist worth understanding. To predict ahead, auto.arima() needs the Fourier terms for the future weeks too, and here is the good news: those terms depend only on the calendar, never on anything unknown. You generate them the same way you generated the training terms, this time asking for a horizon h. Call fourier(x, K, h) to build the terms for the next h steps, then pass them to forecast() as xreg.

Each row is one future week: a point forecast plus a 95 percent interval showing the range the value is likely to fall in. The point forecasts rise week by week, tracing the start of the next yearly climb. Plotting the whole forecast makes the continuation clear.

The shaded fan is the forecast interval, widening as it reaches further ahead. The central line carries the yearly wave smoothly forward from where the data ends.

How do you model multiple seasonal periods at once?

This is where Fourier terms leave seasonal ARIMA behind for good. High-frequency data often repeats on more than one clock at the same time. Half-hourly electricity demand has a daily rhythm (48 half-hours) and a weekly rhythm (336 half-hours) layered on top of each other. A seasonal ARIMA can track only one of them. With Fourier terms you simply add a set of waves for each period.

The recipe has two parts. First store the series as an msts object (a multiple seasonal time series) that records every period. Then give fourier() a vector of K values, one per period. The built-in taylor series (half-hourly electricity demand in England and Wales) already arrives in that form. We fit it here with tslm(), a time series version of lm() that is fast and does the job.

The K = c(10, 10) gives ten pairs to the daily cycle and ten to the weekly cycle. You might expect 41 coefficients (an intercept plus 2 times 10 daily and 2 times 10 weekly columns), but because the weekly period is exactly seven daily periods (336 = 7 times 48), one weekly harmonic falls right on top of a daily one, so fourier() drops that duplicate pair and 39 remain. That single model now carries both seasons together. Forecasting works just as before: build future Fourier terms for both periods and pass them along.

Those four numbers are the predicted demand for the first four half-hours of the next day. The forecast rolls the daily peak-and-trough forward while also respecting where we are in the week, something no single-season model could manage.

How do you do this the modern way with fable?

The forecast package is the classic toolkit. The newer tidyverts family (the tsibble and fable packages) does the same modelling in tidy, pipe-friendly form, and it is what Hyndman's textbook now teaches. The idea is identical: a fourier() term for the season and PDQ(0, 0, 0) to switch off seasonal ARIMA so the Fourier terms own the season. The only real difference is that the data lives in a tsibble, a time-aware version of a data frame, instead of a ts object.

Because tsibble and fable are not part of the in-browser toolkit, run the block below in your own R session after installing them with install.packages(c("fable", "tsibble", "dplyr")).

library(dplyr)
library(tsibble)
library(fable)

demand_tsbl <- tsibble(week = 1:260,
                       demand = as.numeric(weekly_demand),
                       index = week)

fit_fbl <- demand_tsbl |>
  model(dhr = ARIMA(demand ~ fourier(period = 52, K = 3) + PDQ(0, 0, 0)))

report(fit_fbl)
#> Series: demand 
#> Model: LM w/ ARIMA(1,1,1) errors 
#> 
#> Coefficients:
#>          ar1      ma1  fourier(period = 52, K = 3)C1_52  fourier(period = 52, K = 3)S1_52
#>       0.0116  -0.9100                            5.9218                           13.6375
#> s.e.  0.0654   0.0208                            0.3169                            0.3282
#>       fourier(period = 52, K = 3)C2_52  fourier(period = 52, K = 3)S2_52
#>                                 3.7825                            4.3746
#> s.e.                            0.2682                            0.2715
#>       fourier(period = 52, K = 3)C3_52  fourier(period = 52, K = 3)S3_52
#>                                 0.3312                            2.8872
#> s.e.                            0.2580                            0.2595
#> 
#> sigma^2 estimated as 8.984:  log likelihood=-648.62
#> AIC=1315.25   AICc=1315.97   BIC=1347.26

Look at the Fourier coefficients: S1_52 is 13.64 and C1_52 is 5.92, matching the forecast-package fit almost to the decimal. Both packages estimate the same model; they just wrap it in different syntax.

Practice Exercises

These exercises combine the ideas above. Try each before opening the solution.

Exercise 1: Let AICc pick K for you

Instead of reading the AICc table by eye, write code that fits K = 1 through 6 and returns the K with the lowest AICc automatically. Save it to best_k.

Exercise 2: A full multiple-seasonality forecast

Put the whole multiple-seasonality workflow together on the taylor series: fit a harmonic regression with K = c(10, 10), forecast the next full day (48 half-hours), and print the very first forecast value. Save the forecast object to my_fc.

Exercise 3: Check the model is trustworthy

A good forecast model leaves behind residuals that look like random noise. Run checkresiduals() on the K = 3 model dhr_fit from earlier and read the Ljung-Box p-value from the printout. A p-value above 0.05 means the leftovers are indistinguishable from noise, which is what you want.

Frequently Asked Questions

How are Fourier terms different from seasonal dummy variables?

Seasonal dummy variables use one indicator column per step in the season, so a weekly season needs 51 dummies and a daily one needs 364. Fourier terms compress the same shape into a few smooth waves, so the model stays small and the fitted season stays smooth instead of jagged.

What value of K should I start with?

Start small. Try K = 1 up to about 6 or 10, fitting each model and keeping the K with the lowest AICc. Note that K can never exceed m/2 (half the season length), because beyond that point the extra waves simply repeat ones you already have.

Can I add other predictors alongside the Fourier terms?

Yes. The Fourier columns are just regressors, so you can cbind() them with anything else, such as a temperature series or a holiday flag, and pass the combined matrix as xreg. This is how forecasters model electricity demand from the weather and the calendar together.

When should I reach for TBATS or Prophet instead?

Dynamic harmonic regression assumes the seasonal shape stays roughly constant over time. If the seasonality clearly drifts from year to year, TBATS (also in the forecast package) lets it evolve, and Prophet is a popular choice for business series with irregular holidays. Fourier-term regression is the simplest and most transparent place to start.

Do Fourier terms work for daily data with a yearly cycle?

Yes, and it is one of their best uses. A yearly cycle in daily data has a period of about 365, far too long for seasonal ARIMA, yet a handful of Fourier pairs handle it with ease. Set the period to 365.25 (or use an msts object) and pick K by AICc as usual.

Summary

Fourier terms turn a long or complicated season into a few sine and cosine columns that drop straight into a regression. That one trick handles the cases a seasonal ARIMA cannot, namely very long periods and several overlapping cycles, using a model that stays small and forecasts cleanly.

With those five moves you can model weekly sales or half-hourly demand without ever fighting a giant seasonal ARIMA.

References

1. Hyndman, R.J., and Athanasopoulos, G. Forecasting: Principles and Practice (3rd ed), Section 10.5, Dynamic harmonic regression. Link
2. Hyndman, R.J., and Athanasopoulos, G. Forecasting: Principles and Practice (3rd ed), Section 12.1, Complex seasonality. Link
3. forecast package documentation, fourier() reference. Link
4. fable package documentation, ARIMA() model reference. Link
5. forecast package documentation, tslm() reference. Link
6. Hyndman, R.J., and Athanasopoulos, G. Forecasting: Principles and Practice (3rd ed), book home. Link

Continue Learning

• SARIMA in R: How to Fit Seasonal ARIMA Models. The seasonal model that Fourier terms extend to long periods. Start here if the (P,D,Q)[m] notation is new to you.
• ARIMAX in R: ARIMA with External Regressors. Fourier terms are one kind of external regressor; this shows the general xreg idea with real predictors.
• Dynamic Regression in R: Regression with ARIMA Errors. The broader family of models that dynamic harmonic regression belongs to.