时序分析与预测完全指南

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

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

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

它是静止的吗?

有季节性吗?

目标变量是否自相关?

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

预测未来是困难的

自相关

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

自相关示例

上面是一个自相关的例子。仔细观察,你会发现第一个值和第 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.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)  

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)  

ts_ax.set_title('Time Series Analysis Plots\n Dickey-Fuller: p='.format(p_value))  

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):  

seasonal_order=(param[2], D, param[3], s)).fit(disp=-1)  

except:  

continue  

aic = model.aic  

#Save best model, AIC and parameters  

if aicpan>

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]  

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()pan>

# 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()

我们将只关注氮氧化物浓度。因此,我们删除所有其他不相关的列。

然后,我们导入 Prophet。

Prophet 要求日期列命名为 ds,特征列命名为 y,因此我们进行了适当的更改。

此时,我们的数据如下:

然后,我们定义一个训练集。为此,我们将保留最后 30 个条目进行预测和验证。

之后,我们简单地初始化 Prophet,将模型与数据匹配,并进行预测!

你会看到:

这里,yhat 代表预测值,yhat_lower 和 yhat_upper 分别代表预测值的下限和上限。

Prophet 让你可以轻松绘制预测图,我们得到:

NOx 浓度预测

如你所见,Prophet 只是用一条直线来预测未来的 NOx 浓度。

然后,我们检查时间序列是否具有某些有趣的特性,例如季节性:

在这里,Prophet 没有发现季节性的趋势。

通过计算模型的平均绝对百分误差(MAPE)和平均绝对误差(MAE)来评估模型的性能,我们发现 MAPE 为 13.86%,MAE 为 109.32,这还不错!记住,我们根本没有对模型进行微调。

最后,我们只需绘制预测的上限和下限:

每周 NOx 平均浓度预测

恭喜你达到目的!这篇文章很长,但内容丰富。你学会了如何强有力地分析和建模时间序列,并将你的知识应用到两个不同的项目中。我希望你觉得这篇文章有用。

雷锋网雷锋网雷锋网

  • 发表于:
  • 原文链接https://kuaibao.qq.com/s/20190815A0R5TZ00?refer=cp_1026
  • 腾讯「云+社区」是腾讯内容开放平台帐号(企鹅号)传播渠道之一,根据《腾讯内容开放平台服务协议》转载发布内容。

扫码关注云+社区

领取腾讯云代金券