使用Prophet 模型进行时间序列趋势特征提取和预测（一）

Prophet是由Facebook开发的时间序列预测工具，专为处理具有趋势、季节性和节假日效应的数据而设计。它能够捕捉数据的整体趋势、考虑周期性模式，同时允许用户指定特殊日期的影响。Prophet简单易用，具有自动检测变化点的功能，同时提供可解释的结果。无论是销售数据、股价还是其他时间序列数据，都可以用prophet进行预测，但更多的是通过prophet提取特征。
本文包含模型构成和代码详解。

完整代码和数据可关注gzh’finance褪黑素’回复关键词【1012】免费+无套路 获取！

1 Prophet理论详解

(1) 总体趋势模型：

• Prohpet 的趋势模型由线性和饱和增长两部分组成。
• 其中，$g(t)$表示趋势模型，$k$是饱和增长的速率参数，$m(t)$表示线性趋势。

$$g(t)=\frac{1}{1+\exp \left(-k\left(t-t_{\text{saturation}}\right)\right)}\cdot m(t)$$

(2) 季节性成分：

• 季节性成分包括每日、每周和每年的季节性。
• $s(t)$表示季节性成分。

$$s(t)=\sum _{n=1}^N\left(a_n\cdot \sin \left(\frac{2\pi nt}{T}\right)+b_n\cdot \cos \left(\frac{2\pi nt}{T}\right)\right)$$

其中，$N$是季节性的阶数，$T$是季节性的周期。

(3) 假期效应：

• Prohpet 考虑用户定义的假期效应和内置的一些常见节假日。
• $h(t)$表示假期效应。

$$h(t)=\sum _j\left(k_j\cdot I(t\in {\text{holiday}}_j)\right)$$

其中，$k_j$是节假日效应的强度参数，$I(\cdot )$是指示函数。

(4) 整体模型：

• Prohpet 的整体模型是趋势、季节性和假期效应的组合。

$$y(t)=g(t)+s(t)+h(t)+{\epsilon }_t$$

其中，$y(t)$是时间序列的观测值，${\epsilon }_t$是误差项。

运行环境：python3.8及以下 注意高版本不能使用fbprophet哦！

2 导入数据和相应的库

import baostock as bs
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from datetime import datetime, date

import seaborn as sns
import matplotlib.pyplot as plt
from fbprophet import Prophet
from sklearn.metrics import mean_squared_error, mean_absolute_error
sns.set_style("whitegrid")

pjme = pd.read_csv('PJME_hourly.csv',
                   index_col=[0], parse_dates=[0]) # 设置index为日期

3.数据预处理和探查

create_features从时间序列数据中提取一些时间特征。这些特征包括小时、星期几、季度、月份、年份、一年中的第几天、月中的第几天和一年中的第几周。

def create_features(df, label=None):
    """
    从时间序列中提取时间特征
    """
    df = df.copy()
    df['date'] = df.index
    df['hour'] = df['date'].dt.hour
    df['dayofweek'] = df['date'].dt.dayofweek
    df['quarter'] = df['date'].dt.quarter
    df['month'] = df['date'].dt.month
    df['year'] = df['date'].dt.year
    df['dayofyear'] = df['date'].dt.dayofyear
    df['dayofmonth'] = df['date'].dt.day
    df['weekofyear'] = df['date'].dt.weekofyear
    
    X = df[['hour','dayofweek','quarter','month','year',
           'dayofyear','dayofmonth','weekofyear']]
    if label:
        y = df[label]
        return X, y
    return X

X, y = create_features(pjme, label='PJME_MW')

features_and_target = pd.concat([X, y], axis=1)

使用Seaborn 的 pairplot ，绘制上面特征的变化：

sns.pairplot(features_and_target.dropna(),
             hue='hour',
             x_vars=['hour','dayofweek',
                     'year','weekofyear'],
             y_vars='PJME_MW',
             height=5,
             plot_kws={'alpha':0.15, 'linewidth':0}
            )
plt.suptitle('Power Use MW by Hour, Day of Week, Year and Week of Year')
plt.show()

4.划分数据集并使用prophet进行预测

split_date = '01-Jan-2015'
pjme_train = pjme.loc[pjme.index <= split_date].copy()
pjme_test = pjme.loc[pjme.index > split_date].copy()

# 绘制训练集和测试集的划分情况
pjme_test \
    .rename(columns={'PJME_MW': 'TEST SET'}) \
    .join(pjme_train.rename(columns={'PJME_MW': 'TRAINING SET'}),
          how='outer') \
    .plot(figsize=(15,5), title='PJM East', style='.')
plt.show()

# 格式化数据为ds和y用于进行prophet模型的预测
pjme_train.reset_index() \
    .rename(columns={'Datetime':'ds',
                     'PJME_MW':'y'}).head()

# 训练模型
model = Prophet()
model.fit(pjme_train.reset_index() \
              .rename(columns={'Datetime':'ds',
                               'PJME_MW':'y'}))

# 进行预测
pjme_test_fcst = model.predict(df=pjme_test.reset_index() \
                                   .rename(columns={'Datetime':'ds'}))

# 可视化展示结果
f, ax = plt.subplots(1)
f.set_figheight(5)
f.set_figwidth(15)
fig = model.plot(pjme_test_fcst,
                 ax=ax)
plt.show()

# 绘制模型的组成部分
fig = model.plot_components(pjme_test_fcst)

plot_components 会生成一个包含趋势、季节性、年度和周度成分的图表。这些成分有助于直观地了解时间序列的结构和周期性。
趋势：显示时间序列的总体趋势，可能包括增长或下降。
季节性：显示数据中的季节性模式，即在特定时间段内的重复模式。这可以是每周、每月或每年的季节性。
年度成分：显示一年中的周期性成分，通常与季节性有关，但可能更长。
周度成分： 显示一周中的周期性成分，特别是对于包含每周重复模式的数据。

# 比较预测与实际情况
f, ax = plt.subplots(1)
f.set_figheight(5)
f.set_figwidth(15)
ax.scatter(pjme_test.index, pjme_test['PJME_MW'], color='r')
fig = model.plot(pjme_test_fcst, ax=ax)

# 看看第一个月的预测
f, ax = plt.subplots(1)
f.set_figheight(5)
f.set_figwidth(15)
ax.scatter(pjme_test.index, pjme_test['PJME_MW'], color='r')
fig = model.plot(pjme_test_fcst, ax=ax)
ax.set_xbound(lower='01-01-2015',
              upper='02-01-2015')
ax.set_ylim(0, 60000)
plot = plt.suptitle('January 2015 Forecast vs Actuals')

完整代码和数据可关注gzh’finance褪黑素’回复关键词【1012】免费+无套路 获取！