关于机器学习:基于趋势和季节性的时间序列预测

工夫序列预测是基于工夫数据进行预测的工作。它包含建设模型来进行观测，并在诸如天气、工程、经济、金融或商业预测等利用中推动将来的决策。

本文次要介绍工夫序列预测并形容任何工夫序列的两种次要模式(趋势和季节性)。并基于这些模式对工夫序列进行合成。最初应用一个被称为 Holt-Winters 节令办法的预测模型，来预测有趋势和 / 或节令成分的工夫序列数据。

为了涵盖所有这些内容，咱们将应用一个工夫序列数据集，包含 1981 年至 1991 年期间墨尔本 (澳大利亚) 的温度。这个数据集能够从这个 Kaggle 下载，也能够本文最初的 GitHub 下载，其中蕴含本文的数据和代码。该数据集由澳大利亚政府气象局托管，并依据该局的“默认应用条款”(Open Access Licence)取得许可。

首先，导入运行代码所需的下列库。除了最典型的库之外，该代码还基于 statsmomodels 库提供的函数，该库提供了用于预计许多不同统计模型的类和函数，如统计测试和预测模型。

from datetime import datetime, timedelta
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.tsa.stattools import adfuller, kpss
from statsmodels.tsa.api import ExponentialSmoothing
%matplotlib inline

上面是导入数据的代码。数据由两列组成，一列是日期，另一列是 1981 年至 1991 年间墨尔本 (澳大利亚) 的温度。

# date
numdays = 365*10 + 2
base = '2010-01-01'
base = datetime.strptime(base, '%Y-%m-%d')
date_list = [base + timedelta(days=x) for x in range(numdays)]
date_list = np.array(date_list)
print(len(date_list), date_list[0], date_list[-1])

# temp
x = np.linspace(-np.pi, np.pi, 365)
temp_year = (np.sin(x) + 1.0) * 15
x = np.linspace(-np.pi, np.pi, 366)
temp_leap_year = (np.sin(x) + 1.0)
temp_s = []
for i in range(2010, 2020):
    if i == 2010:
        temp_s = temp_year + np.random.rand(365) * 20
    elif i in [2012, 2016]:
        temp_s = np.concatenate((temp_s, temp_leap_year * 15 + np.random.rand(366) * 20 + i % 2010))
    else:
        temp_s = np.concatenate((temp_s, temp_year + np.random.rand(365) * 20 + i % 2010))
print(len(temp_s))

# df
data = np.concatenate((date_list.reshape(-1, 1), temp_s.reshape(-1, 1)), axis=1)
df_orig = pd.DataFrame(data, columns=['date', 'temp'])
df_orig['date'] =  pd.to_datetime(df_orig['date'], format='%Y-%m-%d')
df = df_orig.set_index('date')
df.sort_index(inplace=True)
df

在咱们开始剖析工夫序列的模式之前，让咱们将每个垂直虚线对应于一年开始的数据可视化。

ax = df_orig.plot(x='date', y='temp', figsize=(12,6))
xcoords = ['2010-01-01', '2011-01-01', '2012-01-01', '2013-01-01', '2014-01-01',
           '2015-01-01', '2016-01-01', '2017-01-01', '2018-01-01', '2019-01-01']
for xc in xcoords:
    plt.axvline(x=xc, color='black', linestyle='--')
ax.set_ylabel('temperature')

在进入下一节之前，让咱们花点工夫看一下数据。这些数据仿佛有一个季节性的变动，夏季温度回升，冬季温度降落(南半球)。而且气温仿佛不会随着工夫的推移而减少，因为无论哪一年的平均气温都是雷同的。

工夫序列预测模型应用数学方程 (s) 在一系列历史数据中找到模式。而后应用这些方程将数据[中的历史工夫模式投射到将来。

有四种类型的工夫序列模式:

趋势: 数据的长期增减。趋势能够是任何函数，如线性或指数，并能够随工夫改变方向。

季节性: 以固定的频率 (一天中的小时、星期、月、年等) 在系列中反复的周期。节令模式存在一个固定的已知周期

周期性: 当数据涨跌时产生，但没有固定的频率和持续时间，例如由经济情况引起。

乐音: 系列中的随机变动。

大多数工夫序列数据将蕴含一个或多个模式，但可能不是全副。这里有一些例子，咱们能够辨认这些工夫序列模式:

维基百科年度受众(左图): 在这张图中，咱们能够确定一个减少的趋势，受众每年线性减少。

美国用电量季节性图(中图): 每条线对应的是一年，因而咱们能够察看到每年的用电量反复呈现的季节性。

ibex35 的日收盘价(右图): 这个工夫序列随着工夫的推移有一个减少的趋势，以及一个周期性的模式，因为有一些期间 ibex35 因为经济起因而降落。

如果咱们假如对这些模式进行加法合成，咱们能够这样写:

Y[t] = t [t] + S[t] + e[t]

其中 Y[t]为数据，t [t]为趋势周期重量，S[t]为节令重量，e[t]为噪声，t 为工夫周期。

另一方面，乘法合成能够写成:

Y[t] = t [t] S[t] e[t]

当节令稳定不随工夫序列程度变动时，加法合成是最合适的办法。相同，当节令成分的变动与工夫序列程度成正比时，则采纳乘法合成更为适合。

安稳工夫序列被定义为其不依赖于察看到该序列的工夫。因而具备趋势或季节性的工夫序列不是安稳的，而白噪声序列是安稳的。从数学意义上讲，如果一个工夫序列的均值和方差不变，且协方差与工夫无关，那么这个工夫序列就是安稳的。有不同的例子来比拟安稳和非安稳工夫序列。一般来说，安稳工夫序列不会有长期可预测的模式。

为什么稳定性很重要呢?

平稳性曾经成为工夫序列剖析中许多实际和工具的常见假如。其中包含趋势预计、预测和因果推断等。因而，在许多状况下，须要确定数据是否是由固定过程生成的，并将其转换为具备该过程生成的样本的属性。

如何测验工夫序列的平稳性呢?

咱们能够用两种办法来测验。一方面，咱们能够通过查看工夫序列的均值和方差来手动查看。另一方面，咱们能够应用测试函数来评估平稳性。

有些状况可能会让人感到困惑。例如一个没有趋势和季节性但具备周期行为的工夫序列是安稳的，因为周期的长度不是固定的。

查看趋势

为了剖析工夫序列的趋势，咱们首先应用有 30 天窗口的滚动均值办法剖析随时间推移的平均值。

def analyze_stationarity(timeseries, title):
    fig, ax = plt.subplots(2, 1, figsize=(16, 8))

    rolmean = pd.Series(timeseries).rolling(window=30).mean() 
    rolstd = pd.Series(timeseries).rolling(window=30).std()
    ax[0].plot(timeseries, label= title)
    ax[0].plot(rolmean, label='rolling mean');
    ax[0].plot(rolstd, label='rolling std (x10)');
    ax[0].set_title('30-day window')
    ax[0].legend()
    
    rolmean = pd.Series(timeseries).rolling(window=365).mean() 
    rolstd = pd.Series(timeseries).rolling(window=365).std()
    ax[1].plot(timeseries, label= title)
    ax[1].plot(rolmean, label='rolling mean');
    ax[1].plot(rolstd, label='rolling std (x10)');
    ax[1].set_title('365-day window')
    
pd.options.display.float_format = '{:.8f}'.format
analyze_stationarity(df['temp'], 'raw data')
    ax[1].legend()

在上图中，咱们能够看到应用 30 天窗口时滚动均值是如何随工夫稳定的，这是由数据的季节性模式引起的。此外，当应用 365 天窗口时，滚动平均值随工夫减少，表明随工夫略有减少的趋势。

这也能够通过一些测试来评估，如 Dickey-Fuller (ADF)和 Kwiatkowski, Phillips, Schmidt 和 Shin (KPSS):

ADF 测验的后果 (p 值低于 0.05) 表明，存在的原假如能够在 95% 的置信水平上被回绝。因而，如果 p 值低于 0.05，则工夫序列是安稳的。

def ADF_test(timeseries):
    print("Results of Dickey-Fuller Test:")
    dftest = adfuller(timeseries, autolag="AIC")
    dfoutput = pd.Series(dftest[0:4],
        index=[
            "Test Statistic",
            "p-value",
            "Lags Used",
            "Number of Observations Used",
        ],
    )
    for key, value in dftest[4].items():
        dfoutput["Critical Value (%s)" % key] = value
    print(dfoutput)

ADF_test(df)

Results of Dickey-Fuller Test:
Test Statistic                  -3.69171446
p-value                          0.00423122
Lags Used                       30.00000000
Number of Observations Used   3621.00000000
Critical Value (1%)             -3.43215722
Critical Value (5%)             -2.86233853
Critical Value (10%)            -2.56719507
dtype: float64

KPSS 测验的后果 (p 值高于 0.05) 表明，在 95% 的置信水平下，不能回绝的零假如。因而如果 p 值低于 0.05，则工夫序列不是安稳的。

def KPSS_test(timeseries):
    print("Results of KPSS Test:")
    kpsstest = kpss(timeseries.dropna(), regression="c", nlags="auto")    
    kpss_output = pd.Series(kpsstest[0:3], index=["Test Statistic", "p-value", "Lags Used"]
    )
    for key, value in kpsstest[3].items():
        kpss_output["Critical Value (%s)" % key] = value
    print(kpss_output)

KPSS_test(df)

Results of KPSS Test:
Test Statistic           1.04843270
p-value                  0.01000000
Lags Used               37.00000000
Critical Value (10%)     0.34700000
Critical Value (5%)      0.46300000
Critical Value (2.5%)    0.57400000
Critical Value (1%)      0.73900000
dtype: float64

尽管这些测验仿佛是为了检查数据的平稳性，但这些测验对于剖析工夫序列的趋势而不是季节性是有用的。

统计后果还显示了工夫序列的平稳性的影响。尽管两个测验的零假如是相同的。ADF 测验表明工夫序列是安稳的(p 值 > 0.05)，而 KPSS 测验表明工夫序列不是安稳的(p 值 > 0.05)。但这个数据集创立时带有轻微的趋势，因而结果表明，KPSS 测试对于剖析这个数据集更精确。

为了缩小数据集的趋势，咱们能够应用以下办法打消趋势:

df_detrend = (df - df.rolling(window=365).mean()) / df.rolling(window=365).std()

analyze_stationarity(df_detrend['temp'].dropna(), 'detrended data')
ADF_test(df_detrend.dropna())

查看季节性

正如在之前从滑动窗口中察看到的，在咱们的工夫序列中有一个节令模式。因而应该采纳差分办法来去除工夫序列中潜在的节令或周期模式。因为样本数据集具备 12 个月的季节性，我应用了 365 个滞后差值:

df_365lag =  df - df.shift(365)

analyze_stationarity(df_365lag['temp'].dropna(), '12 lag differenced data')
ADF_test(df_365lag.dropna())

当初，滑动均值和标准差随着工夫的推移或多或少放弃不变，所以咱们有一个安稳的工夫序列。

下面办法合在一起的代码如下：

df_365lag_detrend =  df_detrend - df_detrend.shift(365)

analyze_stationarity(df_365lag_detrend['temp'].dropna(), '12 lag differenced de-trended data')
ADF_test(df_365lag_detrend.dropna())

合成模式

基于上述模式的合成能够通过 ’ statmodels ‘ 包中的一个有用的 Python 函数 seasonal_decomposition 来实现:

def seasonal_decompose (df):
    decomposition = sm.tsa.seasonal_decompose(df, model='additive', freq=365)

    trend = decomposition.trend
    seasonal = decomposition.seasonal
    residual = decomposition.resid

    fig = decomposition.plot()
    fig.set_size_inches(14, 7)
    plt.show()

    return trend, seasonal, residual

seasonal_decompose(df)

在看了合成图的四个局部后，能够说，在咱们的工夫序列中有很强的年度季节性成分，以及随时间推移的减少趋势模式。

工夫序列数据的适当模型将取决于数据的特定特色，例如，数据集是否具备总体趋势或季节性。请务必抉择最适数据的模型。

咱们能够应用上面这些模型：

Autoregression (AR)
Moving Average (MA)
Autoregressive Moving Average (ARMA)
Autoregressive Integrated Moving Average (ARIMA)
Seasonal Autoregressive Integrated Moving-Average (SARIMA)
Seasonal Autoregressive Integrated Moving-Average with Exogenous Regressors (SARIMAX)
Vector Autoregression (VAR)
Vector Autoregression Moving-Average (VARMA)
Vector Autoregression Moving-Average with Exogenous Regressors (VARMAX)
Simple Exponential Smoothing (SES)
Holt Winter’s Exponential Smoothing (HWES)

因为咱们的数据中存在季节性，因而抉择 HWES，因为它实用于具备趋势和 / 或节令成分的工夫序列数据。

这种办法应用指数平滑来编码大量的过来的值，并应用它们来预测当初和将来的“典型”值。指数平滑指的是应用指数加权挪动均匀(EWMA)“平滑”一个工夫序列。

在实现它之前，让咱们创立训练和测试数据集:

y = df['temp'].astype(float)
y_to_train = y[:'2017-12-31']
y_to_val = y['2018-01-01':]
predict_date = len(y) - len(y[:'2017-12-31'])

上面是应用均方根误差 (RMSE) 作为评估模型误差的度量的实现。

def holt_win_sea(y, y_to_train, y_to_test, seasonal_period, predict_date):

    fit1 = ExponentialSmoothing(y_to_train, seasonal_periods=seasonal_period, trend='add', seasonal='add').fit(use_boxcox=True)
    fcast1 = fit1.forecast(predict_date).rename('Additive')
    mse1 = ((fcast1 - y_to_test.values) ** 2).mean()
    print('The Root Mean Squared Error of additive trend, additive seasonal of'+ 
          'period season_length={} and a Box-Cox transformation {}'.format(seasonal_period,round(np.sqrt(mse1), 2)))

    y.plot(marker='o', color='black', legend=True, figsize=(10, 5))
    fit1.fittedvalues.plot(style='--', color='red', label='train')
    fcast1.plot(style='--', color='green', label='test')
    plt.ylabel('temp')
    plt.title('Additive trend and seasonal')
    plt.legend()
    plt.show()

holt_win_sea(y, y_to_train, y_to_val, 365, predict_date)

The Root Mean Squared Error of additive trend, additive seasonal of period season_length=365 and a Box-Cox transformation 6.27

从图中咱们能够察看到模型是如何捕获工夫序列的季节性和趋势的，在异样值的预测上则存在一些误差。

在本文中，咱们通过一个基于温度数据集的理论示例来介绍趋势和季节性。除了查看趋势和季节性之外，咱们还看到了如何升高它，以及如何创立一个根本模型，利用这些模式来推断将来几天的温度。

理解次要的工夫序列模式和学习如何实现工夫序列预测模型是至关重要的，因为它们有许多利用。

本文数据和代码：

https://avoid.overfit.cn/post/51c2316b0237445fbb3dbf6228ea3a52

作者：Javier Fernandez

关于机器学习:基于趋势和季节性的时间序列预测

导入库和数据

可视化数据集

工夫序列模式

合成数据

时序建模

总结

Just My Socks（注册教程内含优惠码）

关于机器学习:基于趋势和季节性的时间序列预测

导入库和数据

可视化数据集

工夫序列模式

合成数据

时序建模

总结

Just My Socks（注册教程 内含优惠码）

Just My Socks（注册教程内含优惠码）