Analyse de la saisonnalité des séries chronologiques des tendances de Google à l'aide de la FFT

Question

J'essaie d'évaluer le spectre d'amplitude de la série chronologique des tendances Google en utilisant une transformation de Fourier rapide. Si vous regardez les données pour "régime" dans les données fournies aquí il présente un schéma saisonnier très fort :

Je pensais pouvoir analyser ce modèle en utilisant une FFT, qui devrait vraisemblablement présenter un pic important pour une période d'un an.

Cependant, lorsque j'applique une FFT comme ceci ( a_gtrend_ham étant la série temporelle multipliée par une fenêtre de Hamming) :

import matplotlib.pyplot as plt
import numpy as np
from numpy.fft import fft, fftshift
import pandas as pd

gtrend = pd.read_csv('multiTimeline.csv',index_col=0)

gtrend.index = pd.to_datetime(gtrend.index, format='%Y-%m')

# Sampling rate
fs = 12 #Points per year

a_gtrend_orig = gtrend['diet: (Worldwide)']
N_gtrend_orig = len(a_gtrend_orig)
length_gtrend_orig = N_gtrend_orig / fs
t_gtrend_orig = np.linspace(0, length_gtrend_orig, num = N_gtrend_orig, endpoint = False)

a_gtrend_sel = a_gtrend_orig.loc['2005-01-01 00:00:00':'2017-12-01 00:00:00']
N_gtrend = len(a_gtrend_sel)
length_gtrend = N_gtrend / fs
t_gtrend = np.linspace(0, length_gtrend, num = N_gtrend, endpoint = False)

a_gtrend_zero_mean = a_gtrend_sel - np.mean(a_gtrend_sel)
ham = np.hamming(len(a_gtrend_zero_mean))
a_gtrend_ham = a_gtrend_zero_mean * ham

N_gtrend = len(a_gtrend_ham)
ampl_gtrend = 1/N_gtrend * abs(fft(a_gtrend_ham))
mag_gtrend = fftshift(ampl_gtrend)
freq_gtrend = np.linspace(-0.5, 0.5, len(ampl_gtrend))
response_gtrend = 20 * np.log10(mag_gtrend)
response_gtrend = np.clip(response_gtrend, -100, 100)

Le spectre d'amplitude qui en résulte ne présente pas de pic dominant :

Où se situe mon incompréhension de la manière d'utiliser la FFT pour obtenir le spectre de la série de données ?

Answer 1

Voici une implémentation propre de ce que je pense que vous essayez d'accomplir. J'ai inclus une sortie graphique et une brève discussion sur ce que cela signifie probablement.

Tout d'abord, nous utilisons le rfft() car les données sont à valeur réelle. Cela permet d'économiser le temps et les efforts (et de réduire le taux de bogues) qui découlent autrement de la génération des fréquences négatives redondantes. Et nous utilisons rfftfreq() pour générer la liste de fréquences (encore une fois, il n'est pas nécessaire de la coder manuellement, et l'utilisation de l'API réduit le taux d'erreurs).

Pour vos données, la fenêtre de Tukey est plus appropriée que la fenêtre de Hamming et les fonctions similaires basées sur cos ou sin. Remarquez également que nous soustrayons la médiane avant de la multiplier par la fonction fenêtre. La médiane() est une estimation assez robuste de la ligne de base, certainement plus que la moyenne().

Dans le graphique, vous pouvez voir que les données chutent rapidement à partir de leur valeur initiale et finissent par être faibles. Les fenêtres de Hamming et similaires échantillonnent le milieu de manière trop étroite pour cela et atténuent inutilement un grand nombre de données utiles.

Pour les graphiques FT, nous sautons la case de la fréquence zéro (le premier point) car elle ne contient que la ligne de base et son omission permet une mise à l'échelle plus commode pour les axes des ordonnées.

Vous remarquerez quelques composantes de haute fréquence dans le graphique de la sortie FT. J'inclus un exemple de code ci-dessous qui illustre une origine possible de ces composantes haute fréquence.

Ok, voici le code :

import matplotlib.pyplot as plt
import numpy as np
from numpy.fft import rfft, rfftfreq
from scipy.signal import tukey

from numpy.fft import fft, fftshift
import pandas as pd

gtrend = pd.read_csv('multiTimeline.csv',index_col=0,skiprows=2)
#print(gtrend)

gtrend.index = pd.to_datetime(gtrend.index, format='%Y-%m')
#print(gtrend.index)

a_gtrend_orig = gtrend['diet: (Worldwide)']
t_gtrend_orig = np.linspace( 0, len(a_gtrend_orig)/12, len(a_gtrend_orig), endpoint=False )

a_gtrend_windowed = (a_gtrend_orig-np.median( a_gtrend_orig ))*tukey( len(a_gtrend_orig) )

plt.subplot( 2, 1, 1 )
plt.plot( t_gtrend_orig, a_gtrend_orig, label='raw data'  )
plt.plot( t_gtrend_orig, a_gtrend_windowed, label='windowed data' )
plt.xlabel( 'years' )
plt.legend()

a_gtrend_psd = abs(rfft( a_gtrend_orig ))
a_gtrend_psdtukey = abs(rfft( a_gtrend_windowed ) )

# Notice that we assert the delta-time here,
# It would be better to get it from the data.
a_gtrend_freqs = rfftfreq( len(a_gtrend_orig), d = 1./12. )

# For the PSD graph, we skip the first two points, this brings us more into a useful scale
# those points represent the baseline (or mean), and are usually not relevant to the analysis
plt.subplot( 2, 1, 2 )
plt.plot( a_gtrend_freqs[1:], a_gtrend_psd[1:], label='psd raw data' )
plt.plot( a_gtrend_freqs[1:], a_gtrend_psdtukey[1:], label='windowed psd' )
plt.xlabel( 'frequency ($yr^{-1}$)' )
plt.legend()

plt.tight_layout()
plt.show()

Et voici la sortie affichée graphiquement. Il y a des signaux forts à 1/an et à 0,14 (qui se trouve être la moitié de 1/14 ans), et il y a un ensemble de signaux de plus haute fréquence qui, à première vue, peuvent sembler assez mystérieux.

Nous voyons que la fonction de fenêtrage est en fait assez efficace pour ramener les données à la ligne de base et vous voyez que les forces relatives du signal dans le FT ne sont pas très modifiées par l'application de la fonction de fenêtrage.

Si vous examinez les données de près, il semble y avoir des variations répétées au cours de l'année. Si ces variations se produisent avec une certaine régularité, on peut s'attendre à ce qu'elles apparaissent comme des signaux dans le FT, et en effet la présence ou l'absence de signaux dans le FT est souvent utilisée pour distinguer le signal du bruit. Mais comme nous allons le montrer, il existe une meilleure explication pour les signaux à haute fréquence.

Voici maintenant un exemple de code qui illustre une façon de produire ces composantes haute fréquence. Dans ce code, nous créons un seul ton, puis nous créons un ensemble de pics à la même fréquence que le ton. Ensuite, nous faisons une transformation de Fourier des deux signaux et enfin, nous traçons un graphique des données brutes et de la FT.

import matplotlib.pyplot as plt
import numpy as np
from numpy.fft import rfft, rfftfreq

t = np.linspace( 0, 1, 1000. )

y = np.cos( 50*3.14*t )

y2 = [ 1. if 1.-v < 0.01 else 0. for v in y ]

plt.subplot( 2, 1, 1 )
plt.plot( t, y, label='tone' )
plt.plot( t, y2, label='spikes' )
plt.xlabel('time')

plt.subplot( 2, 1, 2 )
plt.plot( rfftfreq(len(y),d=1/100.), abs( rfft(y) ), label='tone' )
plt.plot( rfftfreq(len(y2),d=1/100.), abs( rfft(y2) ), label='spikes' )
plt.xlabel('frequency')
plt.legend()

plt.tight_layout()
plt.show()

Ok, voici les graphiques du ton, et des pics, et ensuite leurs transformées de Fourier. Remarquez que les pics produisent des composantes de haute fréquence qui sont très similaires à celles de nos données.

En d'autres termes, l'origine des composantes haute fréquence se situe très probablement dans les échelles de temps courtes associées au caractère piquant des signaux dans les données brutes.