时序分析与预测完全指南

本文作者：skura

2019-08-15 18:45

导语：了解移动平均、指数平滑、平稳性、自相关、时间序列，将这些技术应用于两个项目中。

时序分析与预测完全指南

无论我们是想预测金融市场的趋势还是用电量，时间都是我们模型中必须考虑的一个重要因素。例如，预测一天中什么时候会出现用电高峰是很有趣的，可以以此为依据调整电价或发电量。

输入时间序列。时间序列只是按时间顺序排列的一系列数据点。在时间序列中，时间往往是独立变量，其目标通常是预测未来。

然而，在处理时间序列时，还有一些其他因素会发挥作用。

它是静止的吗？

有季节性吗？

目标变量是否自相关？

在这篇文章中，我将介绍时间序列的不同特征，以及我们如何对它们进行建模才能获得准确的预测。

时序分析与预测完全指南

预测未来是困难的

自相关

通俗地说，自相关是观测值之间的相似度，它是观测值之间时间滞后的函数。

时序分析与预测完全指南

自相关示例

上面是一个自相关的例子。仔细观察，你会发现第一个值和第 24 个值具有很高的自相关性。同样，第 12 个值和第 36 个观测值也高度相关。这意味着我们将在每 24 个时间单位中找到一个非常相似的值。

注意，这个图看起来像正弦函数。这是季节性的征兆，你可以通过在上面的图中找到 24 小时的周期来找到它的价值。

季节性

季节性是指周期性波动。例如，白天的用电量高，晚上的用电量低，或者圣诞节期间的在线销售额增加，节后销售再次放缓。

时序分析与预测完全指南

季节性示例

如你所见，每天都有明显的季节性。每天晚上，你都会看到一个高峰，最低点出现在每天的开始和结束。
记住，如果季节性是满足正弦函数的，它也可以从自相关图中推导出来。简单地看一下周期，它给出了季节的长度。

平稳性

平稳性是时间序列的一个重要特征。如果时间序列的统计性质不随时间变化，则称其为平稳的。换句话说，它有不变的均值和方差，协方差不随时间变化。

时序分析与预测完全指南

平稳过程示例

再看上面的图，我们看到上面的过程是平稳的，平均值和方差不会随时间变化。

通常，股票价格不是一个平稳的过程，因为我们可能会看到一个增长的趋势，或者，其波动性可能会随着时间的推移而增加（这意味着方差正在变化）。

理想情况下，我们需要一个用于建模的固定时间序列。当然，不是所有的都是平稳的，但是我们可以通过做不同的变换，使它们保持平稳。

如何测试过程是否平稳

你可能已经注意到在上图的标题「Dickey-Fuller」。这是我们用来确定时间序列是否稳定的统计测试。

在不讨论 Dickey-Fuller 测试的技术特性的情况下，它检测了单位根是否存在空假设。

如果是，则 P>0，并且过程不是平稳的。

否则，p=0，无效假设被拒绝，过程被认为是平稳的。

例如，下面的过程不是平稳的。请注意为什么平均值不随时间变化。

时序分析与预测完全指南

非平稳过程示例

时间序列建模

有很多方法可以模拟时间序列来进行预测。在此，我将介绍：

移动平均
指数平滑
ARIMA

移动平均

移动平均模型可能是最简单的时间序列建模方法。这个模型简单来说就是，下一个值是所有过去值的平均值。

虽然很简单，但是这个模型的效果可能好到出乎意料，它代表了一个好的起点。

否则，移动平均值可用于识别数据中有趣的趋势。我们可以定义一个窗口来应用移动平均模型来平滑时间序列，并突出不同的趋势。

时序分析与预测完全指南

24 小时窗口上的移动平均值示例

在上面的图中，我们将移动平均模型应用于一个 24 小时窗口。绿线平滑了时间序列，我们可以看到 24 小时内有 2 个峰值。

当然，窗口越长，趋势就越平滑。下面是一个较小窗口上移动平均值的示例。

时序分析与预测完全指南

12 小时窗口上的移动平均值示例

指数平滑

指数平滑使用与移动平均相似的逻辑，但这次，对每个观测值分配了不同的递减权重。换言之，离现在的时间距离越远，观察结果的重要性就越低。

在数学上，指数平滑表示为：

时序分析与预测完全指南

指数平滑表达式

这里，alpha 是一个平滑因子，它的值介于 0 和 1 之间。它决定了之前观测值的权重下降的速度。

时序分析与预测完全指南

指数平滑示例

在上面的图中，深蓝色线表示时间序列的指数平滑，平滑系数为 0.3，而橙色线表示平滑系数为 0.05。

如你所见，平滑因子越小，时间序列就越平滑。这是有意义的，因为当平滑因子接近 0 时，我们接近移动平均模型。

双指数平滑

当时间序列中存在趋势时，使用双指数平滑。在这种情况下，我们使用这种技术，它只是指数平滑的两次递归使用。

数学公式为：

时序分析与预测完全指南

双指数平滑表达式

这里，beta 是趋势平滑因子，它的值介于 0 和 1 之间。

下面，你可以看到 alpha 和 beta 的不同值如何影响时间序列的形状。

时序分析与预测完全指南

双指数平滑示例

三指数平滑

该方法通过添加季节平滑因子来扩展双指数平滑。当然，如果你注意到时间序列中的季节性，这很有用。

在数学上，三指数平滑表示为：

时序分析与预测完全指南

三指数平滑表达式

其中 gamma 是季节平滑因子，L 是季节长度。

季节性差分自回归滑动平均模型（SARIMA）

SARIMA 实际上是简单模型的组合，可以生成一个复杂的模型，该模型可以模拟具有非平稳特性和季节性的时间序列。

首先，我们得到了自回归模型 AR(p)。这基本上是时间序列对自身的回归。在这里，我们假设当前值依赖于它以前的值，并且有一定的滞后。它采用一个表示最大滞后的参数 p。为了找到它，我们查看了部分自相关图，在此之后大部分滞后并不显著。

在下面的例子中，p 的值是 4。

时序分析与预测完全指南

部分自相关图示例

然后，我们添加移动平均模型 MA(q)。这需要一个参数 q，它代表自相关图上那些滞后不显著的最大滞后。

下图中，q 为 4。

时序分析与预测完全指南

自相关图示例

之后，我们添加整合顺序 I(d)。参数 d 表示使序列平稳所需的差异数。

最后，我们添加最后一部分：季节性 S(P, D, Q, s)，其中 S 只是季节的长度。此外，这里要求参数 P 和 Q 与 p 和 q 相同，但用于季节部分。最后，D 是季节整合的顺序，表示从系列中删除季节性所需的差异数量。

综合起来，我们得到了 SARIMA(p, d, q)(P, D, Q, s) 模型。

要注意的是：在用 SARIMA 建模之前，我们必须对时间序列进行转换，以消除季节性和任何非平稳行为。

这是一个很好的理论！让我们在第一个项目中应用上面讨论的技术。

我们将想办法预测一家公司的股票价格。现在，预测股票价格几乎是不可能的。然而，这仍然是一个有趣的练习，它将是一个很好的来实践我们所学到知识的方法。

项目 1：股票价格预测

我们将利用 New Germany Fund（GF）的历史股价来预测未来五个交易日的收盘价。

你可以在这里获取数据集和资料。

像往常一样，我强烈推荐你动手编码！启动你的笔记本，我们开始吧！

时序分析与预测完全指南

你绝不会因为这个项目而发财

导入数据

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()

from sklearn.metrics import r2_score, median_absolute_error, mean_absolute_error
from sklearn.metrics import median_absolute_error, mean_squared_error, mean_squared_log_error

from scipy.optimize import minimize
import statsmodels.tsa.api as smt
import statsmodels.api as sm

from tqdm import tqdm_notebook

from itertools import product

def mean_absolute_percentage_error(y_true, y_pred):
    return np.mean(np.abs((y_true - y_pred) / y_true)) * 100

import warnings
warnings.filterwarnings('ignore')

%matplotlib inline

DATAPATH = 'data/stock_prices_sample.csv'

data = pd.read_csv(DATAPATH, index_col=['DATE'], parse_dates=['DATE'])
data.head(10)

首先，我们导入一些库，这些库将在整个分析过程中都会用到。此外，我们用平均百分比误差（MAPE）作为我们的误差度量。

然后，我们导入数据集，排在前十的是：

时序分析与预测完全指南

数据集的前 10 个条目

正如你所看到的，我们有一些关于 New Germany Fund (GF) 不同股票的数据。此外，我们还有一个关于当天信息的数据，但我们只需要当天结束（EOD）时的股票信息。

数据清洗

data = data[data.TICKER != 'GEF']
data = data[data.TYPE != 'Intraday']

drop_cols = ['SPLIT_RATIO', 'EX_DIVIDEND', 'ADJ_FACTOR', 'ADJ_VOLUME', 'ADJ_CLOSE', 'ADJ_LOW', 'ADJ_HIGH', 'ADJ_OPEN', 'VOLUME', 'FREQUENCY', 'TYPE', 'FIGI']

data.drop(drop_cols, axis=1, inplace=True)

data.head()

首先，我们删除不需要的条目。

然后，我们删除不需要的列，因为我们只想关注股票的收盘价。

如果预览数据集，则应该看到的是：

时序分析与预测完全指南

清洗后的数据集

令人惊叹！我们准备好进行探索性数据分析了！

探索性数据分析（EDA）

# Plot closing price

plt.figure(figsize=(17, 8))
plt.plot(data.CLOSE)
plt.title('Closing price of New Germany Fund Inc (GF)')
plt.ylabel('Closing price ($)')
plt.xlabel('Trading day')
plt.grid(False)
plt.show()

我们绘制数据集整个时间段的收盘价。

你将会得到：

时序分析与预测完全指南

New Germany Fund （GF）收盘价

很明显，你看到的不是一个平稳的过程，很难判断是否有某种季节性。

移动平均

让我们使用移动平均模型来平滑我们的时间序列。为此，我们将使用一个辅助函数，该函数将在指定的时间窗口上运行移动平均模型，并绘制结果平滑曲线：

def plot_moving_average(series, window, plot_intervals=False, scale=1.96):

    rolling_mean = series.rolling(window=window).mean()

    plt.figure(figsize=(17,8))
    plt.title('Moving average\n window size = {}'.format(window))
    plt.plot(rolling_mean, 'g', label='Rolling mean trend')

    #Plot confidence intervals for smoothed values
    if plot_intervals:
        mae = mean_absolute_error(series[window:], rolling_mean[window:])
        deviation = np.std(series[window:] - rolling_mean[window:])
       lower_bound = rolling_mean - (mae + scale * deviation)
       upper_bound = rolling_mean + (mae + scale * deviation)
       plt.plot(upper_bound, 'r--', label='Upper bound / Lower bound')
       plt.plot(lower_bound, 'r--')

    plt.plot(series[window:], label='Actual values')
    plt.legend(loc='best')
    plt.grid(True)

#Smooth by the previous 5 days (by week)
plot_moving_average(data.CLOSE, 5)

#Smooth by the previous month (30 days)
plot_moving_average(data.CLOSE, 30)

#Smooth by previous quarter (90 days)
plot_moving_average(data.CLOSE, 90, plot_intervals=True)

使用5天的时间窗口，我们得到：

时序分析与预测完全指南

上一个交易周的平滑曲线

如你所见，我们几乎看不到趋势，因为它太接近实际曲线。让我们看看上个月和上个季度的平滑处理结果。

时序分析与预测完全指南

上个月（30 天前）的平滑曲线

时序分析与预测完全指南

按上一季度（90 天）平滑

现在更容易发现趋势。请注意，30 天和 90 天的趋势图在末尾显示一条向下的曲线。这意味着股票可能在接下来的几天内会下跌。

指数平滑

现在，让我们用指数平滑来看看它是否能获得更好的趋势。

def exponential_smoothing(series, alpha):

result = [series[0]] # first value is same as series
for n in range(1, len(series)):
result.append(alpha * series[n] + (1 - alpha) * result[n-1])
return result

def plot_exponential_smoothing(series, alphas):

plt.figure(figsize=(17, 8))
for alpha in alphas:
plt.plot(exponential_smoothing(series, alpha), label="Alpha {}".format(alpha))
plt.plot(series.values, "c", label = "Actual")
plt.legend(loc="best")
plt.axis('tight')
plt.title("Exponential Smoothing")
plt.grid(True);

plot_exponential_smoothing(data.CLOSE, [0.05, 0.3])

这里，我们使用 0.05 和 0.3 作为平滑因子的值。当然你也可以尝试其他值，看看结果如何。

时序分析与预测完全指南

指数平滑

如您所见，alpha 值 0.05 平滑了曲线，同时剔除了大部分向上和向下的趋势。

现在，让我们使用双指数平滑。

双指数平滑

def double_exponential_smoothing(series, alpha, beta):

result = [series[0]]
for n in range(1, len(series)+1):
if n == 1:
level, trend = series[0], series[1] - series[0]
if n >= len(series): # forecasting
value = result[-1]
else:
value = series[n]
last_level, level = level, alpha * value + (1 - alpha) * (level + trend)
trend = beta * (level - last_level) + (1 - beta) * trend
result.append(level + trend)
return result

def plot_double_exponential_smoothing(series, alphas, betas):

plt.figure(figsize=(17, 8))
for alpha in alphas:
for beta in betas:
plt.plot(double_exponential_smoothing(series, alpha, beta), label="Alpha {}, beta {}".format(alpha, beta))
plt.plot(series.values, label = "Actual")
plt.legend(loc="best")
plt.axis('tight')
plt.title("Double Exponential Smoothing")
plt.grid(True)

plot_double_exponential_smoothing(data.CLOSE, alphas=[0.9, 0.02], betas=[0.9, 0.02])

你将得到：

时序分析与预测完全指南

双指数平滑

同样，用不同的 α 和 β 组合进行实验，以获得更好的曲线。

建模

如前所述，我们必须将序列转换为一个平稳的过程，以便对其进行建模。因此，让我们应用 Dickey-Fuller 测试来看看它是否是一个平稳的过程：

def tsplot(y, lags=None, figsize=(12, 7), syle='bmh'):

    if not isinstance(y, pd.Series):
        y = pd.Series(y)

   with plt.style.context(style='bmh'):
       fig = plt.figure(figsize=figsize)
       layout = (2,2)
       ts_ax = plt.subplot2grid(layout, (0,0), colspan=2)
       acf_ax = plt.subplot2grid(layout, (1,0))
       pacf_ax = plt.subplot2grid(layout, (1,1))

       y.plot(ax=ts_ax)
       p_value = sm.tsa.stattools.adfuller(y)[1]
       ts_ax.set_title('Time Series Analysis Plots\n Dickey-Fuller: p={0:.5f}'.format(p_value))
       smt.graphics.plot_acf(y, lags=lags, ax=acf_ax)
       smt.graphics.plot_pacf(y, lags=lags, ax=pacf_ax)
       plt.tight_layout()

tsplot(data.CLOSE, lags=30)

# Take the first difference to remove to make the process stationary
data_diff = data.CLOSE - data.CLOSE.shift(1)

tsplot(data_diff[1:], lags=30)

你将看到：

时序分析与预测完全指南

通过 DickeyFuller 测试，时间序列是非平稳的。另外，从自相关图来看，我们发现它似乎没有明显的季节性。

因此，为了消除高度自相关并使过程稳定，让我们取第一个差异（代码块中的第 23 行）。我们简单地用一天的滞后时间减去时间序列，得到：

时序分析与预测完全指南

令人惊叹的！我们的序列现在是平稳的，可以开始建模了！

SARIMA

#Set initial values and some bounds
ps = range(0, 5)
d = 1
qs = range(0, 5)
Ps = range(0, 5)
D = 1
Qs = range(0, 5)
s = 5

#Create a list with all possible combinations of parameters
parameters = product(ps, qs, Ps, Qs)
parameters_list = list(parameters)
len(parameters_list)

# Train many SARIMA models to find the best set of parameters
def optimize_SARIMA(parameters_list, d, D, s):
"""
Return dataframe with parameters and corresponding AIC

parameters_list - list with (p, q, P, Q) tuples
d - integration order
D - seasonal integration order
s - length of season
"""

results = []
best_aic = float('inf')

for param in tqdm_notebook(parameters_list):
try: model = sm.tsa.statespace.SARIMAX(data.CLOSE, order=(param[0], d, param[1]),
seasonal_order=(param[2], D, param[3], s)).fit(disp=-1)
except:
continue

aic = model.aic

#Save best model, AIC and parameters
if aic < best_aic:
best_model = model
best_aic = aic
best_param = param
results.append([param, model.aic])

result_table = pd.DataFrame(results)
result_table.columns = ['parameters', 'aic']
#Sort in ascending order, lower AIC is better
result_table = result_table.sort_values(by='aic', ascending=True).reset_index(drop=True)

return result_table

result_table = optimize_SARIMA(parameters_list, d, D, s)

#Set parameters that give the lowest AIC (Akaike Information Criteria)
p, q, P, Q = result_table.parameters[0]

best_model = sm.tsa.statespace.SARIMAX(data.CLOSE, order=(p, d, q),
seasonal_order=(P, D, Q, s)).fit(disp=-1)

print(best_model.summary())

现在，对于 SARIMA，我们首先定义一些参数值的范围，以生成 p, q, d, P, Q, D, s 的所有可能组合的列表。

现在，在上面的代码单元中，我们有 625 种不同的组合！我们将尝试每种组合，并训练 SARIMA，以便找到性能最佳的模型。这可能需要一些时间，具体多长时间取决于计算机的处理能力。

完成后，我们将输出最佳模型的摘要，你将看到：

时序分析与预测完全指南

令人惊叹！最后，我们预测未来五个交易日的收盘价，并评估模型的 MAPE。

在这种情况下，有一个 0.79% 的 MAPE，这是非常好的！

将预测价格与实际数据进行比较

# Make a dataframe containing actual and predicted prices
comparison = pd.DataFrame({'actual': [18.93, 19.23, 19.08, 19.17, 19.11, 19.12],
'predicted': [18.96, 18.97, 18.96, 18.92, 18.94, 18.92]},
index = pd.date_range(start='2018-06-05', periods=6,))

#Plot predicted vs actual price

plt.figure(figsize=(17, 8))
plt.plot(comparison.actual)
plt.plot(comparison.predicted)
plt.title('Predicted closing price of New Germany Fund Inc (GF)')
plt.ylabel('Closing price ($)')
plt.xlabel('Trading day')
plt.legend(loc='best')
plt.grid(False)
plt.show()

现在，为了将我们的预测与实际数据进行比较，我们从雅虎财务（YahooFinance）获取财务数据并创建一个数据框架。

然后，我们绘出曲线，看看我们与实际收盘价的差距有多大：

时序分析与预测完全指南

预计值和实际收盘价比较

我们的预测似乎有点偏离。事实上，预测价格很平稳，这意味着我们的模型可能表现不佳。

当然，这不是因为我们的程序，而是因为预测股票价格基本上是不可能的。

从第一个项目开始，我们学习了在使用 SARIMA 建模之前平滑时间序列的整个过程。

现在，让我们介绍一下 Facebook 的 Prophet。它是一个在 python 和 r 中都可用的预测工具。该工具帮助生成高质量的预测。

让我们看看如何在第二个项目中使用它！

项目2-使用 Prophet 预测空气质量

标题说明了一切：我们将使用 Prophet 来帮助我们预测空气质量！

时序分析与预测完全指南

导入数据

import warnings
warnings.filterwarnings('ignore')

import numpy as np
import pandas as pd
from scipy import stats
import statsmodels.api as sm
import matplotlib.pyplot as plt

%matplotlib inline

DATAPATH = 'data/AirQualityUCI.csv'

data = pd.read_csv(DATAPATH, sep=';')
data.head()

打印出前五行：

时序分析与预测完全指南

如你所见，数据集包含有关不同气体浓度的信息。每天隔一个小时记录一次。

如果更深入地研究数据集，会发现有许多值 -200 的实例。当然，负浓度是没有意义的，所以我们需要在建模前清洗数据。

数据清洗与特征工程

# Make dates actual dates
data['Date'] = pd.to_datetime(data['Date'])

# Convert measurements to floats
for col in data.iloc[:,2:].columns:
if data[col].dtypes == object:
data[col] = data[col].str.replace(',', '.').astype('float')

# Compute the average considering only the positive values
def positive_average(num):
return num[num > -200].mean()

# Aggregate data
daily_data = data.drop('Time', axis=1).groupby('Date').apply(positive_average)

# Drop columns with more than 8 NaN
daily_data = daily_data.iloc[:,(daily_data.isna().sum() <= 8).values]

# Remove rows containing NaN values
daily_data = daily_data.dropna()

# Aggregate data by week
weekly_data = daily_data.resample('W').mean()

# Plot the weekly concentration of each gas
def plot_data(col):
plt.figure(figsize=(17, 8))
plt.plot(weekly_data[col])
plt.xlabel('Time')
plt.ylabel(col)
plt.grid(False)
plt.show()

for col in weekly_data.columns:
plot_data(col)

在这里，我们首先分析日期列，将其转换为日期类型。

然后，我们把所有的测量值转换成浮点数。

之后，我们用每天的平均值来汇总数据。

我们还有一些需要删除的 NAN。

最后，我们按周汇总数据，这将提供一个更平滑的分析趋势。

我们可以画出每种化学物质浓度的趋势。这里，我们展示了 NOx。

时序分析与预测完全指南

NOx 浓度

氮氧化物是非常有害的，因为它们会形成烟雾和酸雨，同时也会形成细颗粒和臭氧。这些都会对健康产生不利影响，因此氮氧化物的浓度是空气质量的一个关键特征。

建模

# Drop irrelevant columns
cols_to_drop = ['PT08.S1(CO)', 'C6H6(GT)', 'PT08.S2(NMHC)', 'PT08.S4(NO2)', 'PT08.S5(O3)', 'T', 'RH', 'AH']

weekly_data = weekly_data.drop(cols_to_drop, axis=1)

# Import Prophet
from fbprophet import Prophet
import logging

logging.getLogger().setLevel(logging.ERROR)

# Change the column names according to Prophet's guidelines
df = weekly_data.reset_index()
df.columns = ['ds', 'y']
df.head()

# Split into a train/test set
prediction_size = 30
train_df = df[:-prediction_size]

# Initialize and train a model
m = Prophet()
m.fit(train_df)

# Make predictions
future = m.make_future_dataframe(periods=prediction_size)
forecast = m.predict(future)
forecast.head()

# Plot forecast
m.plot(forecast)

# Plot forecast's components
m.plot_components(forecast)

# Evaluate the model
def make_comparison_dataframe(historical, forecast):
return forecast.set_index('ds')[['yhat', 'yhat_lower', 'yhat_upper']].join(historical.set_index('ds'))

cmp_df = make_comparison_dataframe(df, forecast)
cmp_df.head()

def calculate_forecast_errors(df, prediction_size):

df = df.copy()

df['e'] = df['y'] - df['yhat']
df['p'] = 100 * df['e'] / df['y']

predicted_part = df[-prediction_size:]

error_mean = lambda error_name: np.mean(np.abs(predicted_part[error_name]))

return {'MAPE': error_mean('p'), 'MAE': error_mean('e')}

for err_name, err_value in calculate_forecast_errors(cmp_df, prediction_size).items():
print(err_name, err_value)

# Plot forecast with upper and lower bounds
plt.figure(figsize=(17, 8))
plt.plot(cmp_df['yhat'])
plt.plot(cmp_df['yhat_lower'])
plt.plot(cmp_df['yhat_upper'])
plt.plot(cmp_df['y'])
plt.xlabel('Time')
plt.ylabel('Average Weekly NOx Concentration')
plt.grid(False)
plt.show()