7. Časové řady

• Vyzkoušet si práci s knihovnou Pandas
• Seznámit se s principem klouzavého průměru
• Dekompozice časové řady

• Cheat Sheet: The pandas DataFrame Object pdf
• Zajímavé datové soubory: link

Vyhlazování (smoothing) dat pomocí klouzavého průměru je jednoduchá, ale efektivní technika zpracování časových řad. Cílem vyhlazování je odstranit šum a lépe odhalit původní signál. Klouzavé průměry jsou jednoduchým a běžným typem vyhlazování používaným v analýze časových řad a při předpovídání časových řad.

Výpočet klouzavého průměru spočívá ve vytvoření nové řady, jejíž hodnoty jsou tvořeny průměrem nezpracovaných pozorování v původní časové řadě. Klíčovým parameterem metody je šířka okna, které se posouvá podél časové osy a v jehož rámci jsou vypočítány průměrné hodnoty v nové řadě.

Nejjednodušší metodou vyhlazení je případ, kdy použijeme pro výpočet nových hodnot několik bodů před a několik bodů po středu okna (centered moving averace). Příklad pro velikost okna 3:

center_ma(t) = mean(obs(t-1), obs(t), obs(t+1))

Metoda potřebuje znát budoucí data, takže se většinou používá na analýzu již pořízených časových řad. Metodu lze použít odstranění trendových a sezónních složek z časové řady, pro predikci použít nelze.

O něco lepší je použít pouze historické hodnoty, zde pro okno velikosti 3:

trail_ma(t) = mean(obs(t-2), obs(t-1), obs(t))

Pro účely cvičení budeme používat dostupnou časovou řadu, obsahující počet porodů v Kalifornii, v roce 1959.

daily-total-female-births.zip

Stáhněte si soubor do pracovního adresáře a rozbalte. Prvních několik řádků vypadá následovně:

"Date","Births"
"1959-01-01",35
"1959-01-02",32
"1959-01-03",30
"1959-01-04",31
"1959-01-05",44

from pandas import read_csv
from matplotlib import pyplot
 
series = read_csv('daily-total-female-births.csv', header=0, index_col=0)
 
print(series.head())
series.plot()
pyplot.show()

Pro vyhlazení dat použijeme funkci rolling(). Tato funkce automaticky vytvoří okno o požadované velikosti a použije ho pro výpočet hodnot nové řady. Následující příklad pracuje s oknem velikosti 3, vyzkoušejte i okna jiné velikosti.

# Tail-rolling average transform
rolling = series.rolling(window=3)
rolling_mean = rolling.mean()
print(rolling_mean.head(10))
 
# plot original and transformed dataset
series.plot()
rolling_mean.plot(color='red')
pyplot.show()

Klouzavý průměr lze použít jako zdroj nových informací při modelování předpovědi časové řady (tzv. supervised learning).

V tomto případě se vypočítaný klouzavý průměr přidá jako nová vstupní funkce použitá k předpovědi dalšího časového kroku. Nejprve je třeba posunout kopii řady o jeden časový krok dopředu. Ta bude představovat vstup do předpovědi neboli verzi řady se zpožděním = 1.

X, 	y
NaN,	obs1
obs1,	obs2
obs2,	obs3

Pro lepší představu graficky pro okno velikosti 3:

from pandas import read_csv
from pandas import DataFrame
from pandas import concat
 
series = read_csv('daily-total-female-births.csv', header=0, index_col=0)
 
df = DataFrame(series.values)
 
width = 3
 
lag1 = df.shift(1)
lag3 = df.shift(width - 1)
 
window = lag3.rolling(window=width)
means = window.mean()
 
dataframe = concat([means, lag1, df], axis=1)
dataframe.columns = ['mean', 't-1', 't+1']
 
print(dataframe.head(10))

from numpy import mean
from sklearn.metrics import mean_squared_error
from matplotlib import pyplot
series = read_csv('daily-total-female-births.csv', header=0, index_col=0)
# prepare situation
X = series.values
window = 3
history = [X[i] for i in range(window)]
test = [X[i] for i in range(window, len(X))]
predictions = list()
# walk forward over time steps in test
for t in range(len(test)):
	length = len(history)
	yhat = mean([history[i] for i in range(length-window,length)])
	obs = test[t]
	predictions.append(yhat)
	history.append(obs)
	print('predicted=%f, expected=%f' % (yhat, obs))
error = mean_squared_error(test, predictions)
print('Test MSE: %.3f' % error)
# plot
pyplot.plot(test)
pyplot.plot(predictions, color='red')
pyplot.show()
# zoom plot
pyplot.plot(test[0:100])
pyplot.plot(predictions[0:100], color='red')
pyplot.show()

Pro toto cvičení použijeme časovou řadu s údaji o měsíčním počtu cestujících v letech 1949 až 1960.

airline-passengers.zip

from pandas import read_csv
from matplotlib import pyplot
series = read_csv('airline-passengers.csv', header=0, index_col=0)
series.plot()
pyplot.show()

import pandas as pd
from matplotlib import pyplot
from statsmodels.tsa.seasonal import seasonal_decompose
series = pd.read_csv('airline-passengers.csv', header=0)
s = series.squeeze()
x = list(s['Passengers'])
y = pd.Series(x, index=pd.date_range("1-1-1949", periods=len(x), freq="M"), name="Passengers")
result = seasonal_decompose(y, model='multiplicative')
result.plot()
pyplot.show()

import matplotlib.pyplot as plt
from statsmodels.tsa.seasonal import STL
from statsmodels.datasets import elec_equip as ds
 
elec_equip = ds.load().data
 
def add_stl_plot(fig, res, legend):
    """Add 3 plots from a second STL fit"""
    axs = fig.get_axes()
    comps = ["trend", "seasonal", "resid"]
    for ax, comp in zip(axs[1:], comps):
        series = getattr(res, comp)
        if comp == "resid":
            ax.plot(series, marker="o", linestyle="none")
        else:
            ax.plot(series)
            if comp == "trend":
                ax.legend(legend, frameon=False)
 
stl = STL(elec_equip, period=12, robust=True)
res_robust = stl.fit()
fig = res_robust.plot()
res_non_robust = STL(elec_equip, period=12, robust=False).fit()
add_stl_plot(fig, res_non_robust, ["Robust", "Non-robust"])
plt.show()