How to determine multiple Periodicities present in Timeseries data?

Question

My objective is to detect all kinds of seasonalities and their time periods that are present in a timeseries waveform.

I'm currently using the following dataset: https://www.kaggle.com/rakannimer/air-passengers

At the moment, I've tried the following approaches:

1) Use of FFT:

import pandas as pd
import numpy as np
from statsmodels.tsa.seasonal import seasonal_decompose
 
#https://www.kaggle.com/rakannimer/air-passengers
df=pd.read_csv('AirPassengers.csv')
 
df.head()

frequency_eval_max = 100
A_signal_rfft = scipy.fft.rfft(df['#Passengers'], n=frequency_eval_max)
n = np.shape(A_signal_rfft)[0] # np.size(t)
frequencies_rel = len(A_signal_fft)/frequency_eval_max * np.linspace(0,1,int(n))

fig=plt.figure(3, figsize=(15,6))
plt.clf()
plt.plot(frequencies_rel, np.abs(A_signal_rfft), lw=1.0, c='paleturquoise')
plt.stem(frequencies_rel, np.abs(A_signal_rfft))
plt.xlabel("frequency")
plt.ylabel("amplitude")

This results in the following plot:

But it doesn't result in anything conclusive or comprehensible.

Ideally I wish to see the peaks representing daily, weekly, monthly and yearly seasonality.

Could anyone point out what am I doing wrong?

2) Autocorrelation:

from pandas.plotting import autocorrelation_plot
plt.rcParams.update({'figure.figsize':(10,6), 'figure.dpi':120})
autocorrelation_plot(df['#Passengers'].tolist())

After doing which I get a plot like the following:

But how do I read this plot and how can I derive the presence of the various seasonalities and their periods from this?

3) SLT Decomposition Algorithm

df.set_index('Month',inplace=True)
df.index=pd.to_datetime(df.index)
#drop null values
df.dropna(inplace=True)
df.plot()

result=seasonal_decompose(df['#Passengers'], model='multiplicable', period=12)

result.seasonal.plot()

This gives the following plot:

But here I can only see one kind of seasonality.

So how do we detect all the types of seasonalities and their time periods that are present using this method?

---

Hence, I've tried 3 different approaches but they seem either erroneous or incomplete.

Could anyone please help me out with the most effective approach (even apart from the ones I've tried) to detect all kinds of seasonalities and their time periods for any given timeseries data?

---

---

EDIT: As suggested by Max, I tried out FFT (with 500 samples transformed in FFT) & Autocorrelation on a different dataset having more no. of data points resulting in the following. (Note: this is a square wave having daily and weekly periodicity)

The following is my code for the same:

1. Generation of the waveform:

import json
import sys, os
import numpy as np
import pandas as pd
import glob
import pickle

from statsmodels.tsa.stattools import adfuller, acf, pacf
from scipy.signal import find_peaks, square
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import matplotlib.pyplot as plt

#GENERATION OF A FUNCTION WITH DUAL SEASONALITY & NOISE

def white_noise(mu, sigma, num_pts):
    """ Function to generate Gaussian Normal Noise
    Args:
        sigma: std value
        num_pts: no of points
        mu: mean value

    Returns:
        generated Gaussian Normal Noise
    """
    
    noise = np.random.normal(mu, sigma, num_pts)
    return noise

def signal_line_plot(input_signal: pd.Series, title: str = "", y_label: str = "Signal"):
    """ Function to plot a time series signal
    Args:
        input_signal: time series signal that you want to plot
        title: title on plot
        y_label: label of the signal being plotted
        
    Returns:
        signal plot
    """
    
    plt.plot(input_signal)
    plt.title(title)
    plt.ylabel(y_label)
    plt.show()

t_week = np.linspace(1,480, 480)
t_weekend=np.linspace(1,192,192)
T=96 #Time Period
x_weekday = 10*square(2*np.pi*t_week/T, duty=0.7)+10 + white_noise(0, 1,480)
x_weekend = 2*square(2*np.pi*t_weekend/T, duty=0.7)+2 + white_noise(0,1,192)
x_daily_weekly = np.concatenate((x_weekday, x_weekend)) 
x_daily_weekly_long = np.concatenate((x_daily_weekly,x_daily_weekly,x_daily_weekly,x_daily_weekly,x_daily_weekly,x_daily_weekly,x_daily_weekly,x_daily_weekly,x_daily_weekly,x_daily_weekly))
signal_line_plot(x_daily_weekly_long)
signal_line_plot(x_daily_weekly_long[0:1000])

#x_daily_weekly_long is the final waveform

2. FFT:

frequency_eval_max = 500
A_signal_rfft = scipy.fft.rfft(x_daily_weekly_long, n=frequency_eval_max)
n = np.shape(A_signal_rfft)[0] # np.size(t)
frequencies_rel = len(A_signal_fft)/frequency_eval_max * np.linspace(0,1,int(n))

fig=plt.figure(3, figsize=(15,6))
plt.clf()
plt.plot(frequencies_rel, np.abs(A_signal_rfft), lw=1.0, c='paleturquoise')
plt.stem(frequencies_rel, np.abs(A_signal_rfft))
plt.xlabel("frequency")
plt.ylabel("amplitude")

3) Autocorrelation:

from pandas.plotting import autocorrelation_plot
plt.rcParams.update({'figure.figsize':(10,6), 'figure.dpi':120})
autocorrelation_plot(x_daily_weekly_long)

Answer 1

Your dataset is sampled on a monthly basis. The shortest seasonality you will possibly be able to see is $2$ months. Daily, weekly and monthly are impossible.

While your FFT plot is somewhat clumsy to interpret due to the fact, that you choose to transform 100 samples and to scale the plot considering the DC part (the first FFT bin, $0\text{Hz}$), some things can be garnered from it. The peak around $0.25$ represents the period of $2$ years length, the one at about $0.5$ is the year long period. The peak around $0.95$ stands for $8$ months. The autocorrelation plot clearly depicts the year long period. That's where the peaks at multiples of 12 come from, so does the STL plot (as you concluded yourself. $8$ months is visible here, too.

EDIT

Your synthesized signal has a sampling rate of 96/day (one sample every quarter hour). The length of the signal is 6780 samples, which translates to 10 weeks. The maximum frequency represented by the FFT plot is 48/day. The frequency axis of the FFT plot has then a bin-to-bin distance of $$\frac{48}{(6780/2)}\cdot\frac{1}{\text{day}} \approx 0,014\cdot\frac{1}{\text{day}}$$

This means, that at sample $1/0,014\approx71$ there should be a peak representing the daily periodicity.

Usually, you would go the other way around: for each peak, you calculate the period its bin stands for. For example the high peak in the 5th bin stands for the two week period. $\frac{1}{0,014\cdot5}\approx 14$

Like this, you can interpret each peak in the FFT.

For the autocorrelation it is analog: you look for peaks and dips, and since it is in the time domain, you can read the period directly from the position. E.g. a peak at sample 96 represents a daily period, since the sampling frequency is 96/day.

If you want a somewhat smoother FFT, try sine signals instead of square waves, as these tend to be "noisy" in the frequency domain (a rect function transforms to an si function).