prophet模型预测时间序列

prophet：facebook开源的的一个时间序列预测算法。
原理：基于时间序列分解和机器学习的拟合来做的。

把时间序列分成不同时间间隔和整体趋势的组合。
间隔：年季度月周日

prophet适用范围

.有至少几个月（最好是一年）的每小时、每天或每周观察的历史数据；
较强的季节性趋势；
事先知道的以不定期的间隔发生的重要节假日；
缺失的历史数据或较大的异常数据的数量在合理范围内；
历史趋势的变化；
数据中蕴含的非线性增长的趋势都有一个自然极限或饱和状态；

prophet使用

输入

输入已知的时间序列的时间戳和相应的值
输入需要预测的时间序列的长度

传入prophet的数据分为两列 ds 和 y ,ds表示时间序列的时间戳，y表示时间序列的取值

输出

输出未来的时间序列走势
输出结果可以提供必要的统计指标，包括拟合曲线，上界和下界等

yhat，表示时间序列的预测值
yhat_lower，表示预测值的下界
yhat_upper，表示预测值的上界

模型数学原理

模型整体由三部分组成：

growth(增长趋势)
seasonality(季节趋势)
holidays(节假日对预测值的影响)

g(t) 表示趋势项，它表示时间序列在非周期上面的变化趋势；
s(t) 表示周期项，或者称为季节项，一般来说是以周或者年为单位；
h(t) 表示节假日项，表示时间序列中那些潜在的具有非固定周期的节假日对预测值造成的影响；
即误差项或者称为剩余项，表示模型未预测到的波动，服从高斯分布；
Prophet 算法就是通过拟合这几项，然后最后把它们累加起来就得到了时间序列的预测值。

趋势项

分为两种：
线性增长：基于分段线性函数。
非线性增长：基于逻辑回归函数。

线性

k 表示增长率

表示增长率的变化量

m 表示偏移量

非线性：

**C(t)**表示承载量:
它是一个随时间变化的函数，限定了所能增长的最大值
在使用 Prophet 的 growth = ‘logistic’ 的时候，需要提前设置好 C(t) 的取值才行

k 表示增长率，在代码中写为 changepoints。
在现实的时间序列中，曲线的走势肯定不会一直保持不变，在某些特定的时候或者有着某种潜在的周期曲线会发生变化，模型定义了增长率k发生变化时对应的点，将其称作changepoints。

如何设置changeppints呢？

人工指定的方式指定变点的位置
算法来自动选择
在默认的函数里面，Prophet 会选择 n_changepoints = 25 个变点，然后设置变点的范围是前 80%
也就是在时间序列的前 80% 的区间内会设置变点。
之后还要看一些边界条件是否合理，例如时间序列的点数是否少于 n_changepoints 等内容；
其次如果边界条件符合，那变点的位置就是均匀分布的

有了变点changeppints之后，比如S是变点，在S个点之内，有一段固定变化的的增长率

m 表示偏移量

线性和非线性区别：

模型调节参数

对于趋势变化，模型需要调节的参数：

changepoint_range：变点的位置
n_changepoint：变点的个数
changepoint_prior_scale：增长的变化率

changepoint_range 指的是百分比，需要在前 changepoint_range 那么长的时间序列中设置变点，
默认的函数中是 changepoint_range = 0.8。

n_changepoint 表示变点的个数，
默认的函数中是 n_changepoint = 25。
changepoint_prior_scale 表示变点增长率的分布情况，

季节项

傅立叶级数来模拟

模型调节参数

seasonality_prior_scale 这个值越大，表示季节的效应越明显；这个值越小，表示季节的效应越不明显。
seasonality_mode 两种模式，分别是加法和乘法，默认是加法的形式。

节假日项

模型调节参数

holidays_prior_scale
默认值是 10，当值越大时，表示节假日对模型的影响越大；当值越小时，表示节假日对模型的效果越小。

介绍一个开源数据库-Quandl 库

注册网站：
https://data.nasdaq.com/
注册完之后，在你的个界面找到API

pip install quandl

从quandl库中免费提取出特斯拉和通用股票的数据

import quandl
import pandas as pd
quandl.ApiConfig.api_key = '这里要写自己的API'
tesla = quandl.get('WIKI/TSLA')
gm = quandl.get('WIKI/GM')
gm.head()

结果如下

这里备注一下（废话可跳过，有卡住的朋友可以看）

本人从quandl网站注册了账号之后，复制API带代码中，一直报403错误
说我的邮件没确认什么的，但是该网站确实也没给我的注册邮箱
哎，只能等过两天看看邮箱有没有收到确认邮件了
所以暂时放弃从quandl上获取数据
果然睡一觉起来，在邮箱里找到确认的信，点了之后，复制的API就可以使用了

安装相应的配置环境

从入门到放弃，安装包就是一大门槛，有时候，就是因为装包，会导致我放弃使用某个模型，这个facebook的prophet，我花费了半天的时间安装包，最终总算找到一条相对简便，通用的方法。

要不是已经花了半天的时间去看这个prophet理解，可能真的就放弃了。

直接使用anaconda了。

这个三角位置点开的话，选择terminal’，可以直接启动相应的虚拟环境

建议直接在下面create一个新的虚拟环境，然后在新的虚拟环境里配置porphet的包

直接刷选conda配置包，将现有包升到最新，
这也是为什么让大家建立一个虚拟环境的原因

conda config --add channels conda-forge
conda update --all highest priority

conda install pystan

conda install fbprophet

conda install plotly -y

安装完成之后就可以开始了

再次明确prophet的预测

输入
输入已知的时间序列的时间戳和相应的值
输入需要预测的时间序列的长度

输出
输出未来的时间序列走势
输出可以提供必要的统计指标，包括拟合曲线，上界和下界

prophet是把时间分为几个部分
季节项 [公式] ，趋势项 [公式] ，剩余项 [公式]

文件要求
csv 文件存储两列，第一列的名字是 ‘ds’ ，表示时间序列的时间戳,；
第二列的名称是 ‘y’，表示时间序列的取值。

prophet 可以计算出 yhat，yhat_lower，yhat_upper，分别表示时间序列的预测值，预测值的下界，预测值的上界。

由于输出太长了，所以截断成两个表输出

Prophet使用注意

Prophet 所需要的两列名称是 ‘ds’ 和 ‘y’，下面需要改一下名字
如果列1是日期，列2是要预测的值

df = df.rename(columns={'列名1（这一列是日期）':'ds', '列名2':'y'})

日期类型，可能需要转化，如果输入数据日期那一列不是日期格式的话

下面是将秒数据，修改成 YYYY-MM-DD hh:mm:ss 的形式

df['ds'] = pd.to_datetime(df['ds'],unit='s')

秒和时间格式转化的理解例子：
下图是秒数和时间互转

当然如果日期是字符串格式也需要转换为日期类型的格式

日期格式的归一化操作

df['y'] = (df['y'] - df['y'].mean()) / (df['y'].std())

代码部分

导包

import pandas as pd
import fbprophet
import matplotlib.pyplot as plt
import numpy as np

从源网站下载数据

import quandl
quandl.ApiConfig.api_key = '我的API，不能告诉你哦！'
tesla = quandl.get('WIKI/TSLA')
gm = quandl.get('WIKI/GM')

查看下载出来的数据

gm.head()

tesla.head()

绘制原始数据的图像用于观察

# 画图查看股票价格
plt.plot(gm.index, gm['Adj. Close'])
plt.title('GM Stock Price')
plt.ylabel('Price ($)');
plt.show()
plt.plot(tesla.index, tesla['Adj. Close'], 'r')
plt.title('Tesla Stock Price')
plt.ylabel('Price ($)');
plt.show();

可以看出从11年到18年，，一共6个年份
自己设置6个每年的平均股份数据
先以tesla为例进行

# 导入年平均的股票数
# 这里e6就是后面6个零，e9就是后面9个零
tesla_shares = {2018: 168e6, 2017: 162e6, 2016: 144e6, 2015: 128e6, 2014: 125e6, 2013: 119e6, 2012: 107e6, 2011: 100e6, 2010: 51e6}

给原有tesla创建一个新的列，设为年份 year

创建另一个新列，设为市值 cap

市值 = 年平均股票数 x 股票价格

# 创建一个新列设为年份
tesla['Year'] = tesla.index.year
#将索引提前出来设为新的列日期
tesla.reset_index(level=0, inplace = True)
# 创建市值列
tesla['cap'] = 0

根据年份将市值cap进行填补

for i, year in enumerate(tesla['Year']):# 将每一年的平均股票数 按 data上year匹配到每个数据shares = tesla_shares.get(year)# 市值价 = 年平均股票数 x 收盘价tesla.loc[i,'cap'] = shares * tesla.loc[i,'Adj. Close']

小tips:

get函数

dict.get(key, default=None)
参数：
key – 字典中要查找的键。
default – 如果指定键的值不存在时，返回该默认值。

get（key）就是从字典中，根据key来获取相应的value值，如果key不存在，也不会报错，就不返回值了。
举例子：

dict = {'Name': 'Runoob', 'Age': 27}
dict.get('Age')
dict.get('Sex', "Not Available")

27
Not Available

get使用要注意！不能加 default= 这几个字！

dict={}
dict.get('python') #返回None
dict.get('python', default=None) #直接报错，这个错误在某些情况没有任何提示！让我找了好久bug
dict.get('python', default='HHHH'）#直接报错
dict.get('python', 'HHHH') 返回HHHH

同样逻辑处理通用汽车

# 同样逻辑处理通用汽车的数据
gm_shares = {2018: 1.42e9, 2017: 1.50e9, 2016: 1.54e9, 2015: 1.59e9, 2014: 1.61e9, 2013: 1.39e9, 2012: 1.57e9, 2011: 1.54e9, 2010:1.50e9}
gm['Year'] = gm.index.year
gm.reset_index(level=0, inplace = True)
gm['cap'] = 0
for j, year in enumerate(gm['Year']):shares = gm_shares.get(year)gm.loc[j, 'cap'] = shares * gm.loc[j, 'Adj. Close']

查看一下处理后，新增的两列数据，cap市值，和year年份

将tesla公司和gm通用公司的市值变化绘制到同一张图

下面想将tesla公司和gm通用公司的市值变化绘制到同一张图
首先是根据时间，将tesla和gm的两个表合并
使用的.merge（）函数

# 合并tesla和gm的数据，通过merge关联，交集为两个公司相同的时间段的数据
cars = gm.merge(tesla, how='inner', on='Date')

将合并后的表的列名/字段进行改名，以便区分

# 将car_x的列名/字段 改为gm_cap
# 将tesla_cap的列名/字段 改为tesla_cap
cars.rename(columns={'cap_x': 'gm_cap', 'cap_y': 'tesla_cap'}, inplace=True)

# 只提取需要的列
cars = cars.loc[:,['Date', 'gm_cap', 'tesla_cap']]

缩放市值的倍数

# 将市值列同步缩小e9即 9个10的倍数
cars['gm_cap'] = cars['gm_cap'] / 1e9
cars['tesla_cap'] = cars['tesla_cap'] / 1e9
cars.head()

将市值绘制出来进行对比

# 将市值绘制出来进行对比
plt.figure(figsize = (10,5),dpi=100) # 设置画布
plt.plot(cars.Date, cars['gm_cap'],label = 'gm') # 绘制gm市值变化
plt.plot(cars.Date, cars['tesla_cap'],label = 'tesla') # 绘制tesla市值
plt.legend() # 添加小插板
plt.title('Market Cap of GM and Tesla') # 设置标题
plt.ylabel('MarketCap(Billions$)') # 设置纵坐标标题
plt.grid(linestyle='-.') # 添加网格
plt.show()

采用prophet进行预测

prophet封装比较严重
数据分为两列
日期列名为ds
值列名为y
因此使用prophet模型首先需要修改列名

后面由于值这一列数组太大，采用全部缩小9个0的倍数

# prophet需要 日期ds 值 y
gm = gm.rename(columns={'Date': 'ds', 'cap': 'y'})
# 缩小y值
gm['y'] = gm['y']/1e9
gm = gm.loc[:,['ds','y']]

训练模型

# 通过propeht训练 gm的数据并进行预测gm_prophet = fbprophet.Prophet(changepoint_prior_scale = 0.15)
gm_prophet.fit(gm)

设置需要预测的未来的时间长度

# 设置预测的未来的时间
gm_forecast = gm_prophet.make_future_dataframe(periods=365 * 2, freq='D')
# 导入数据进行预测
gm_forecast = gm_prophet.predict(gm_forecast)

绘制预测图

# 使用ophet 的绘图函数来可视化预测
gm_prophet.plot(gm_forecast, xlabel = 'Date', ylabel = 'Market Cap (billions $)')
plt.title('Market Cap of GM') # 设置标题

额外拓展-周期变化点

可以查看模型的周期变化点

gm_prophet.changepoints[:10]

查看预测输出的模型的结果

合并两个公司未来的预测的数据为一个表

# 将两个公司的数据合并
# 遍历 预测出来数据的每一行，在前面添加gm_存入list中保存
gm_names = ['gm_%s' % column for column in gm_forecast.columns]
tesla_names = ['tesla_%s' % column for column in tesla_forecast.columns]# 创建备份数据进行合并
# 合并的数据框
merge_gm_forecast = gm_forecast.copy()
merge_tesla_forecast = tesla_forecast.copy()# 根据list更新列名
merge_gm_forecast.columns = gm_names
merge_tesla_forecast.columns = tesla_names# 合并两组数据
forecast = pd.merge(merge_gm_forecast,merge_tesla_forecast,how = 'inner',left_on = 'gm_ds', right_on = 'tesla_ds')# 将日期列更名（删除gm_ds和tesla_ds，将其中一个改名为Date）
forecast = forecast.rename(columns={'gm_ds':'Date'})
freocast = forecast.drop('tesla_ds',axis=1)

forecast

绘制未来变化趋势图进行对比

# 绘制GM和Tesla两个公司趋势对比
plt.figure(figsize=(10,8),dpi=80)
plt.plot(forecast['Date'],forecast['gm_trend'],'b-',label='GM')
plt.plot(forecast['Date'],forecast['tesla_trend'],'r-',label='Tesla')
plt.legend()
plt.grid()
plt.xlabel('Date')
plt.ylabel('Market Cap($)')
plt.title('GM vs Tesla Trend')

# 绘制预测
plt.figure(figsize=(10, 8))
plt.plot(forecast['Date'], forecast['gm_yhat'],'b-',label = 'GM')
plt.plot(forecast['Date'], forecast['tesla_yhat'], 'r-',label = 'Tesla')
plt.legend()
plt.grid()
plt.xlabel('Date')
plt.ylabel('Market Cap (billions $)')
plt.title('GM vs Tesla Estimate')

找到tesla超过GM的首个日期点

# 从预测的曲线上找到Tesla超过GM的日期
overtake_date = min(forecast.loc[forecast['tesla_yhat'] > forecast['gm_yhat'],'Date'])
# overtake_date = forecast.loc[forecast['tesla_yhat'] > forecast['gm_yhat'],'Date']
print(overtake_date)

绘制出预测可能存在的波动性

# 将预测可能存在的波动性情况也画出来
fig, ax = plt.subplots(1, 1, figsize=(10, 8)) # 设置画布大小，设置一个图# 绘制第一个图
ax.plot(forecast['Date'], forecast['gm_yhat'], label = 'gm prediction')
# 绘制一个区间变化的图线，
ax.fill_between(forecast['Date'].dt.to_pydatetime(), forecast['gm_yhat_upper'], forecast['gm_yhat_lower'], alpha=0.6, edgecolor = 'k')# 绘制第二个图
ax.plot(forecast['Date'], forecast['tesla_yhat'], 'r', label = 'tesla prediction');
ax.fill_between(forecast['Date'].dt.to_pydatetime(), forecast['tesla_yhat_upper'], forecast['tesla_yhat_lower'], alpha=0.6, edgecolor = 'k')plt.legend() # 显示右上角小桌板
plt.xlabel('Date')  # 设置x轴 名
plt.ylabel('Billions $') # 设置y轴 名
plt.title('Market Cap Prediction for GM and Tesla') # 设置 标题名

小tip：

plt.figure()参数详解

matpltlib.pyplot.figure(
num = None,               # 设定figure名称。系统默认按数字升序命名的figure_num（透视表输出窗口）e.g. “figure1”。可自行设定figure名称，名称或是INT，或是str类型；
figsize=None,             # 设定figure尺寸。系统默认命令是rcParams["figure.fig.size"] = [6.4, 4.8]，即figure长宽为6.4 * 4.8；
dpi=None,                 # 设定figure像素密度。系统默命令是rcParams["sigure.dpi"] = 100；
facecolor=None,           # 设定figure背景色。系统默认命令是rcParams["figure.facecolor"] = 'w'，即白色white；
edgecolor=None, frameon=True,    # 设定要不要绘制轮廓&轮廓颜色。系统默认绘制轮廓，轮廓染色rcParams["figure.edgecolor"]='w',即白色white；
FigureClass=<class 'matplotlib.figure.Figure'>,   # 设定使不使用一个figure模板。系统默认不使用；
clear=False,                     # 设定当同名figure存在时，是否替换它。系统默认False，即不替换。
**kwargs)

fig,axes 往往是用来绘制在一个图上排列多个图

举例子：

设置一个画布，，分为2x2，就是四块

fig, axes = plt.subplots(2, 2)  # 此处是一个2*2的图

位置分别是
( 0，0) | (0，1)
(1，0) | (1，1)

给图传入数据

data = pd.Series(np.random.rand(16), index=list('abcdefghijklmnop'))

设置子图的位置

data.plot.bar(ax=axes[1,1], color='b', alpha=0.5)  # ax=[1,1] 即位置是第2行、第二列。（python从0开始计数，所以“1”代表第2的）data.plot.barh(ax=axes[0,1], color='k', alpha=0.5) # alpha：设定图表的透明度；

例子2：
导包

import numpy as np
import matplotlib.pyplot as plt

创建画布，设置子图
区分点s
fig, ax1 = plt.subplot()
fig, ax1 = plt.subplots()

 # 做1*1个子图，等价于 " fig, ax1 = plt.subplot() "，等价于 " fig, ax1 = plt.subplots() "
fig, ax1 = plt.subplots(1, 1)

 # 让2个子图的x轴一样，同时创建副坐标轴。ax2 = ax1.twinx()

画图

# 作y=sin(x)函数
x1 = np.linspace(0, 4 * np.pi, 100)
y1 = np.sin(x1)
ax1.plot(x1, y1)
#  作y = cos(x)函数
x2 = np.linspace(0, 4 * np.pi, 100)       # 表示在区间[0, 4π]之间取100个点（作为横坐标，“线段是有无数多个点组成的”）。
y2 = np.cos(x2)
ax2.plot(x2, y2)

保存图

plt.savefig('路径xxx //sin_cos_2.png')               # 将窗口另存为png格式图片

查看GM通用的趋势和规律

plot函数创建了实际/预测的图表，plot_components提供了趋势/季节性的图表。

gm_prophet.plot_components(gm_forecast)

增长趋势模型（trend）

季节性模型（weekly和yearly）

如果在每年趋势”yearly”中你认为当前的效果过拟合了，
那么就可以调解seasonality_prior_scale这个参数，值越小，这里的季节性波动就越小。

调节参数建议

如果预测结果的误差很大，考虑选取的模型是否准确；
尝试调整增长率模型（growth）的参数；
在必要的情况下也需要调整季节性（seasonality）参数。

.如果对历史数据进行仿真预测时发现，
从一个截点到下一个截点误差急剧的增加，
说明在两个截点期间数据的产生过程发生了较大的变化，此时两个截点之间应该增加一个”changepoint”，来对这期间的不同阶段分别建模。

Prophet设置的参数

Capacity：在增量函数是逻辑回归函数的时候，需要设置的容量值。

季节性和节假日：可以根据实际的业务需求来指定相应的节假日。

光滑参数

growth：增长趋势模型。分为两种：”linear”与”logistic”，分别代表线性与非线性的增长，默认值：”linear”.

选择logistic后需要设置 capacity参数
Capacity：在增量函数是逻辑回归函数的时候，需要设置的容量值,表示非线性增长趋势中限定的最大值，预测值将在该点达到饱和.

Change Points：可以通过 n_changepoints 和 changepoint_range 来进行等距的变点设置（模型自己修正），
也可以通过人工设置的方式来指定时间序列的变点，默认值：“None”.

**n_changepoints：**用户指定潜在的”changepoint”的个数，默认值：25。

changepoint_prior_scale：增长趋势模型的灵活度。调节”changepoint”选择的灵活度，
值越大，选择的”changepoint”越多，从而使模型对历史数据的拟合程度变强，然而也增加了过拟合的风险。默认值：0.05。

seasonality_prior_scale（seasonality模型中的）：调节季节性组件的强度。
值越大，模型将适应更强的季节性波动，值越小，越抑制季节性波动，默认值：10.0.

holidays_prior_scale（holidays模型中的）：调节节假日模型组件的强度。
值越大，该节假日对模型的影响越大，值越小，节假日的影响越小，默认值：10.0。

freq：数据中时间的统计单位（频率），默认为”D”，按天统计。

periods：需要预测的未来时间的个数。例如按天统计的数据，想要预测未来一年时间内的情况，则需填写365。

mcmc_samples：mcmc采样，用于获得预测未来的不确定性。
若大于0，将做mcmc样本的全贝叶斯推理，如果为0，将做最大后验估计，默认值：0。

interval_width：衡量未来时间内趋势改变的程度。表示预测未来时使用的趋势间隔出现的频率和幅度与历史数据的相似度，
值越大越相似，默认值：0.80。
当mcmc_samples = 0时，该参数仅用于增长趋势模型的改变程度，
当mcmc_samples > 0时，该参数也包括了季节性趋势改变的程度。

uncertainty_samples：用于估计未来时间的增长趋势间隔的仿真绘制数，默认值：1000。

参考文章：

https://www.cnblogs.com/lfri/p/12245352.html
https://zhuanlan.zhihu.com/p/92615532
https://www.jianshu.com/p/afe4e7c0a7eb
https://www.jianshu.com/p/0c06ad7bccaa
https://facebook.github.io/prophet/docs/quick_start.html#python-api
https://zhuanlan.zhihu.com/p/52330017
https://www.runoob.com/python/att-dictionary-get.html
https://blog.csdn.net/htuhxf/article/details/82986440
https://www.jianshu.com/p/218757bee516