====== 7. Časové řady ======
===== Cíle cvičení =====
* Vyzkoušet si práci s knihovnou Pandas
* Seznámit se s principem klouzavého průměru
* Dekompozice časové řady
===== Doplňkové poklady =====
* Cheat Sheet: The pandas DataFrame Object {{ :courses:b0b37nsi:tutorials:pandas_dataframe_notes.pdf |pdf}}
* Zajímavé datové soubory: [[https://datahub.io/collections|link]]
===== Klouzavý průměr =====
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))
==== Data ====
Pro účely cvičení budeme používat dostupnou časovou řadu, obsahující počet porodů v Kalifornii, v roce 1959.
{{ :courses:b0b37nsi:tutorials: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
==== Načtení a zobrazení dat ====
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()
==== Vyhlazení dat ====
Pro vyhlazení dat použijeme funkci [[https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.Series.rolling.html|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 jako vstup pro predikce ====
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'':
{{:courses:b0b37nsi:tutorials:ma.png?400|}}
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))
==== Přímá predikce ====
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()
===== Dekompozice časové řady =====
==== Data ====
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.
{{ :courses:b0b37nsi:tutorials:airline-passengers.zip |}}
==== Načtení dat ====
from pandas import read_csv
from matplotlib import pyplot
series = read_csv('airline-passengers.csv', header=0, index_col=0)
series.plot()
pyplot.show()
==== Sezónní dekompozice ====
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()
==== Kvalita fitování ====
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()