0
  • 聊天消息
  • 系统消息
  • 评论与回复
登录后你可以
  • 下载海量资料
  • 学习在线课程
  • 观看技术视频
  • 写文章/发帖/加入社区
创作中心

完善资料让更多小伙伴认识你,还能领取20积分哦,立即完善>

3天内不再提示

时间序列分析和预测基础理论知识

数据分析与开发 来源:机器学习研习院 作者:小猴子 2022-03-16 14:05 次阅读

今天给大家带来一篇实战案例,本案例旨在运用之前学习的时间序列分析和预测基础理论知识,用一个基于交通数据的实际案例数据演示这些方法是如何被应用的。

本文主要内容

★ 首先使用探索性数据分析,从不同时间维度探索分析交通系统乘客数量。

★ 创建一个函数来检查时间序列数据的平稳性,通过一阶差分将非平稳性数据转化为平稳性数据。

★ 然后将数据分为训练集和验证集,简单介绍了并应用多个时间序列预测技术,如朴素法、移动平均方法、简单指数平滑、霍尔特线性趋势法、霍尔特-温特法、ARIMA和SARIMAX模型。

★ 最后使用SARIMAX模型预测未来7个月的流量,因为它有最小的RMSE。如下图所示,蓝色线是训练数据,黄色线是验证数据,红色是使用SARIMAX模型预测的数据。

8e6ef822-a3e7-11ec-952b-dac502259ad0.png

虽然拟合最好的SARIMAX模型,但似乎也没那么棒,当然会有更好的方法来预测该数据。而本文重点是介绍这些基于统计学的经典时间序列预测技术在实际案例中的应用。

8e811d36-a3e7-11ec-952b-dac502259ad0.png

导入相关模块

importpandasaspd
importnumpyasnp
importmatplotlib.pyplotasplt
fromdatetimeimportdatetime
frompandasimportSeries
fromsklearn.metricsimportmean_squared_error
frommathimportsqrt
fromstatsmodels.tsa.seasonalimportseasonal_decompose
importstatsmodels
importstatsmodels.apiassm
fromstatsmodels.tsa.arima_modelimportARIMA

数据集准备

去直接读取用pandas读取csv文本文件,并拷贝一份以备用。

train=pd.read_csv("Train.csv")
test=pd.read_csv("Test.csv")

train_org=train.copy()
test_org=test.copy()

查看数据的列名

train.columns,test.columns
(Index(['ID', 'Datetime', 'Count'], dtype='object'),
 Index(['ID', 'Datetime'], dtype='object'))

查看数据类型

train.dtypes,test.dtypes
(ID           int64
 Datetime    object
 Count        int64
 dtype: object,
 ID           int64
 Datetime    object
 dtype: object)

查看数据大小

test.shape,train.shape
((5112, 2), (18288, 3))

查看数据样貌

train.head()
8ee61df8-a3e7-11ec-952b-dac502259ad0.png

解析日期格式

train['Datetime']=pd.to_datetime(train.Datetime,format='%d-%m-%Y%H:%M')
test['Datetime']=pd.to_datetime(test.Datetime,format='%d-%m-%Y%H:%M')
test_org['Datetime']=pd.to_datetime(test_org.Datetime,format='%d-%m-%Y%H:%M')
train_org['Datetime']=pd.to_datetime(train_org.Datetime,format='%d-%m-%Y%H:%M')

时间日期格式解析结束后,记得查看下结果。

train.dtypes
ID                   int64
Datetime    datetime64[ns]
Count                int64
dtype: object
train.head()
8efc8ec6-a3e7-11ec-952b-dac502259ad0.png

时间序列数据的特征工程

时间序列的特征工程一般可以分为以下几类。本次案例我们根据实际情况,选用时间戳衍生时间特征。

8f1368da-a3e7-11ec-952b-dac502259ad0.png

时间戳虽然只有一列,但是也可以根据这个就衍生出很多很多变量了,具体可以分为三大类:时间特征、布尔特征,时间差特征

本案例首先对日期时间进行时间特征处理,而时间特征包括年、季度、月、周、天(一年、一月、一周的第几天)、小时、分钟...

因为需要对test, train, test_org, train_org四个数据框进行同样的处理,直接使用for循环批量提取年月日小时等特征。

foriin(test,train,test_org,train_org):
i['Year']=i.Datetime.dt.year
i['Month']=i.Datetime.dt.month
i['day']=i.Datetime.dt.day
i['Hour']=i.Datetime.dt.hour
#i["dayoftheweek"]=i.Datetime.dt.dayofweek
test.head()
8f320538-a3e7-11ec-952b-dac502259ad0.png

时间戳衍生中,另一常用的方法为布尔特征,即:

是否年初/年末

是否月初/月末

是否周末

是否节假日

是否特殊日期

是否早上/中午/晚上

等等

下面判断是否是周末,进行特征衍生的布尔特征转换。首先提取出日期时间的星期几。

train['dayoftheweek']=train.Datetime.dt.dayofweek
#返回给定日期时间的星期几
train.head()
8f50c914-a3e7-11ec-952b-dac502259ad0.png

再判断day of the week是否是周末(星期六和星期日),是则返回1 ,否则返回0

defapplyer(row):
ifrow.dayofweek==5orrow.dayofweek==6:
return1
else:
return0
temp=train['Datetime']
temp2=train.Datetime.apply(applyer)
train['weekend']=temp2
train.index=train['Datetime']

对年月乘客总数统计后可视化,看看总体变化趋势。

df=train.drop('ID',1)
ts=df['Count']
plt.plot(ts,label='Passengercount')

8f66c5e8-a3e7-11ec-952b-dac502259ad0.png

探索性数据分析

首先使用探索性数据分析,从不同时间维度探索分析交通系统乘客数量。

对年进行聚合,求所有数据中按年计算的每日平均客流量,从图中可以看出,随着时间的增长,每日平均客流量增长迅速。

train.groupby('Year')['Count'].mean().plot.bar()
8fae5d22-a3e7-11ec-952b-dac502259ad0.png

对月份进行聚合,求所有数据中按月计算的每日平均客流量,从图中可以看出,春夏季客流量每月攀升,而秋冬季客流量骤减。

train.groupby('Month')['Count'].mean().plot.bar()
8fd031e0-a3e7-11ec-952b-dac502259ad0.png

年月

对年月份进行聚合,求所有数据中按年月计算的每日平均客流量,从图可知道,几本是按照平滑指数上升的趋势。

temp=train.groupby(['Year','Month'])['Count'].mean()
temp.plot()#乘客人数(每月)
8feac4a6-a3e7-11ec-952b-dac502259ad0.png

对日进行聚合,求所有数据中每月中的每日平均客流量。从图中可大致看出,在5、11、24分别出现三个峰值,该峰值代表了上中旬的高峰期。

train.groupby('day')['Count'].mean(
).plot.bar(figsize=(15,5))
9006a428-a3e7-11ec-952b-dac502259ad0.png

小时

对小时进行聚合,求所有数据中一天内按小时计算的平均客流量,得到了在中(12)晚(19)分别出现两个峰值,该峰值代表了每日的高峰期。

train.groupby('Hour')['Count'].mean().plot.bar()
901769d4-a3e7-11ec-952b-dac502259ad0.png

是否周末

对是否是周末进行聚合,求所有数据中按是否周末计算的平均客流量,发现工作日比周末客流量客流量多近一倍,果然大家都是周末都喜欢宅在家里。

train.groupby('weekend')['Count'].mean().plot.bar()
90344478-a3e7-11ec-952b-dac502259ad0.png

对星期进行聚合统计,求所有数据中按是周计算的平均客流量。

train.groupby('dayoftheweek')['Count'].mean().plot.bar()
9041e830-a3e7-11ec-952b-dac502259ad0.png

时间重采样

重采样(resampling)指的是将时间序列从一个频率转换到另一个频率的处理过程;
◎ 将高频率数据聚合到低频率称为降采样(downsampling);
◎ 将低频率数据转换到高频率则称为升采样(unsampling);

train.head()
9063722a-a3e7-11ec-952b-dac502259ad0.png

Pandas中的resample,重新采样,是对原样本重新处理的一个方法,是一个对常规时间序列数据重新采样和频率转换的便捷的方法。

Resample方法的主要参数,如需要了解详情,可以戳这里了解更多。

参数 说明
freq 表示重采样频率,例如'M'、'5min'、Second(15)
how='mean' 用于产生聚合值的函数名或数组函数,例如'mean'、'ohlc'、np.max等,默认是'mean',其他常用的值由:'first'、'last'、'median'、'max'、'min'
axis=0 默认是纵轴,横轴设置axis=1

接下来对训练数据进行对月、周、日及小时多重采样。其实我们分月份进行采样,然后求月内的均值。事实上重采样,就相当于groupby,只不过是根据月份这个period进行分组。

train=train.drop('ID',1)
train.timestamp=pd.to_datetime(train.Datetime,format='%d-%m-%Y%H:%M')
train.index=train.timestamp

#每小时的时间序列
hourly=train.resample('H').mean()
#换算成日平均值
daily=train.resample('D').mean()
#换算成周平均值
weekly=train.resample('W').mean()
#换算成月平均值
monthly=train.resample('M').mean()

重采样后对其进行可视化,直观地看看其变化趋势。

90728954-a3e7-11ec-952b-dac502259ad0.png

对测试数据也进行相同的时间重采样处理。

test.Timestamp=pd.to_datetime(test.Datetime,format='%d-%m-%Y%H:%M')
test.index=test.Timestamp
#换算成日平均值
test=test.resample('D').mean()
train.Timestamp=pd.to_datetime(train.Datetime,format='%d-%m-%Y%H:%M')
train.index=train.Timestamp

#C换算成日平均值
train=train.resample('D').mean()

划分训练集和验证集

到目前为止,我们有训练集和测试集,实际上,我们还需要一个验证集,用来实时验证和调整训练模型。下面直接用索引切片的方式做处理。

Train=train.loc['2012-08-25':'2014-06-24']
valid=train['2014-06-25':'2014-09-25']

划分好数据集后,绘制折线图将训练集和验证集进行可视化。

9088ae00-a3e7-11ec-952b-dac502259ad0.png

模型建立

数据准备好了,就到了模型建立阶段,这里我们应用多个时间序列预测技术,如朴素法、移动平均方法、简单指数平滑、霍尔特线性趋势法、霍尔特-温特法、ARIMA和SARIMAX模型。

朴素预测法

如果数据集在一段时间内都很稳定,我们想预测第二天的价格,可以取前面一天的价格,预测第二天的值。这种假设第一个预测点和上一个观察点相等的预测方法就叫朴素预测法(Naive Forecast),即。

因为朴素预测法用最近的观测值作为预测值,因此他最简单的预测方法。虽然朴素预测法并不是一个很好的预测方法,但是它可以为其他预测方法提供一个基准。

dd=np.asarray(Train.Count)
#将结构数据转化为ndarray
y_hat=valid.copy()
y_hat['naive']=dd[len(dd)-1]
plt.plot(Train.index,Train['Count'],label='Train')
plt.plot(valid.index,valid['Count'],label='Valid')
plt.plot(y_hat.index,y_hat['naive'],label='NaiveForecast')
90c604bc-a3e7-11ec-952b-dac502259ad0.png

模型评价

用RMSE检验朴素法的的准确率

rms=sqrt(mean_squared_error(valid.Count,y_hat.naive))
print(rms)
111.79050467496724

移动平均值法

移动平均法也叫滑动平均法,取前面n个点的平均值作为预测值。

计算移动平均值涉及到一个有时被称为"滑动窗口"的大小值。使用简单的移动平均模型,我们可以根据之前数值的固定有限数的平均值预测某个时序中的下一个值。利用一个简单的移动平均模型,我们预测一个时间序列中的下一个值是基于先前值的固定有限个数“p”的平均值。

这样,对于所有的

#最近10次观测的移动平均值,即滑动窗口大小为P=10
y_hat_avg=valid.copy()
y_hat_avg['moving_avg_forecast']=Train['Count'].rolling(10).mean().iloc[-1]

#最近20次观测的移动平均值
y_hat_avg=valid.copy()
y_hat_avg['moving_avg_forecast']=Train['Count'].rolling(20).mean().iloc[-1]

#最近30次观测的移动平均值
y_hat_avg=valid.copy()
y_hat_avg['moving_avg_forecast']=Train['Count'].rolling(50).mean().iloc[-1]
plt.plot(Train['Count'],label='Train')
plt.plot(valid['Count'],label='Valid')
plt.plot(y_hat_avg['moving_avg_forecast'],
label='MovingAverageForecastusing50observations')
90e1cd3c-a3e7-11ec-952b-dac502259ad0.png90f81f6a-a3e7-11ec-952b-dac502259ad0.png911e9f78-a3e7-11ec-952b-dac502259ad0.png

简单指数平滑法

介绍这个之前,需要知道什么是简单平均法(Simple Average),该方法预测的期望值等于所有先前观测点的平均值。

物品价格会随机上涨和下跌,平均价格会保持一致。我们经常会遇到一些数据集,虽然在一定时期内出现小幅变动,但每个时间段的平均值确实保持不变。这种情况下,我们可以认为第二天的价格大致和过去的平均价格值一致。

简单平均法和加权移动平均法在选取时间点的思路上存在较大的差异:

简单平均法将过去数据一个不漏地全部加以同等利用;

移动平均法则不考虑较远期的数据,并在加权移动平均法中给予近期更大的权重。

我们就需要在这两种方法之间取一个折中的方法,在将所有数据考虑在内的同时也能给数据赋予不同非权重。

简单指数平滑法 (Simple Exponential Smoothing)相比更早时期内的观测值,越近的观测值会被赋予更大的权重,而时间越久远的权重越小。它通过加权平均值计算出预测值,其中权重随着观测值从早期到晚期的变化呈指数级下降,最小的权重和最早的观测值相关:

其中是平滑参数。

fromstatsmodels.tsa.apiimportExponentialSmoothing,SimpleExpSmoothing,Holt
y_hat_avg=valid.copy()
fit2=SimpleExpSmoothing(np.asarray(Train['Count'])).fit(smoothing_level=0.6,optimized=False)
y_hat_avg['SES']=fit2.forecast(len(valid))
plt.figure(figsize=(16,8))
plt.plot(Train['Count'],label='Train')
plt.plot(valid['Count'],label='Valid')
plt.plot(y_hat_avg['SES'],label='SES')
plt.legend(loc='best')
plt.show()
91301258-a3e7-11ec-952b-dac502259ad0.png

模型评价

用RMSE检验朴素法的的准确率

rms=sqrt(mean_squared_error(valid.Count,y_hat_avg.SES))
print(rms)
113.43708111884514

霍尔特线性趋势法

Holts线性趋势模型,该方法考虑了数据集的趋势,即序列的增加或减少性质。

尽管这些方法中的每一种都可以应用趋势:简单平均法会假设最后两点之间的趋势保持不变,或者我们可以平均所有点之间的所有斜率以获得平均趋势,使用移动趋势平均值或应用指数平滑。但我们需要一种无需任何假设就能准确绘制趋势图的方法。这种考虑数据集趋势的方法称为霍尔特线性趋势法,或者霍尔特指数平滑法。

y_hat_avg=valid.copy()
fit1=Holt(np.asarray(Train['Count'])
).fit(smoothing_level=0.3,smoothing_slope=0.1)
y_hat_avg['Holt_linear']=fit1.forecast(len(valid))

plt.plot(Train['Count'],label='Train')
plt.plot(valid['Count'],label='Valid')
plt.plot(y_hat_avg['Holt_linear'],label='Holt_linear')

9153734c-a3e7-11ec-952b-dac502259ad0.png

模型评价

用RMSE检验朴素法的的准确率

rms=sqrt(mean_squared_error(valid.Count,y_hat_avg.Holt_linear))
print(rms)
112.94278345314041

由于holts线性趋势,到目前为止具有最好的准确性,我们尝试使用它来预测测试数据集。

predict=fit1.forecast(len(test))
test['prediction']=predict
#计算每小时计数的比率
train_org['ratio']=train_org['Count']/train_org['Count'].sum()
#按小时计数分组
temp=train_org.groupby(['Hour'])['ratio'].sum()
#保存聚合后的数据
pd.DataFrame(temp,columns=['ratio']).to_csv('GROUPBY.csv')

temp2=pd.read_csv('GROUPBY.csv')
#按日、月、年合并test和test_org
merge=pd.merge(test,test_org,on=('day','Month','Year'),how='left')
merge['Hour']=merge['Hour_y']
merge['ID']=merge['ID_y']
merge.head()
9167badc-a3e7-11ec-952b-dac502259ad0.png

删除不需要的特征。

merge=merge.drop(['Year','Month','Datetime','Hour_x','Hour_y','ID_x','ID_y'],axis=1)
merge.head()
91784d7a-a3e7-11ec-952b-dac502259ad0.png

通过合并merge和temp2进行预测。

prediction=pd.merge(merge,temp2,on='Hour',how='left')
#将比率转换成原始比例
prediction['Count']=prediction['prediction']*prediction['ratio']*24
submission=prediction
pd.DataFrame(submission,columns=['ID','Count']).to_csv('Holt_Linear.csv')
918eeb84-a3e7-11ec-952b-dac502259ad0.png

霍尔特-温特法

霍尔特-温特(Holt-Winters)方法,在 Holt模型基础上引入了 Winters 周期项(也叫做季节项),可以用来处理月度数据(周期 12)、季度数据(周期 4)、星期数据(周期 7)等时间序列中的固定周期的波动行为。引入多个 Winters 项还可以处理多种周期并存的情况。

#HoltsWintermodel
y_hat_avg=valid.copy()
fit1=ExponentialSmoothing(np.asarray(Train['Count']),seasonal_periods=7,trend='add',seasonal='add',).fit()
y_hat_avg['Holts_Winter']=fit1.forecast(len(valid))
plt.plot(Train['Count'],label='Train')
plt.plot(valid['Count'],label='Valid')
plt.plot(y_hat_avg['Holts_Winter'],label='Holt_Winter')

91a226ae-a3e7-11ec-952b-dac502259ad0.png

模型评价

用RMSE检验朴素法的的准确率

rms=sqrt(mean_squared_error(valid.Count,y_hat_avg.Holts_Winter))
print(rms)
82.37292653831038

模型预测

predict=fit1.forecast(len(test))
test['prediction']=predict
#按日、月、年合并Test和test_original
merge=pd.merge(test,test_org,on=('day','Month','Year'),how='left')
merge['Hour']=merge['Hour_y']
merge=merge.drop(['Year','Month','Datetime','Hour_x','Hour_y'],axis=1)
#通过合并merge和temp2进行预测
prediction=pd.merge(merge,temp2,on='Hour',how='left')
#将比率转换成原始比例
prediction['Count']=prediction['prediction']*prediction['ratio']*24
prediction['ID']=prediction['ID_y']
submission=prediction.drop(['day','Hour','ratio','prediction','ID_x','ID_y'],axis=1)
#转换最终提交的csv格式
pd.DataFrame(submission,columns=['ID','Count']).to_csv('Holtwinters.csv')
91b6f2b4-a3e7-11ec-952b-dac502259ad0.png

迪基-福勒检验

函数执行迪基-福勒检验以确定数据是否为平稳时间序列。

在统计学里,迪基-福勒检验(Dickey-Fuller test)可以测试一个自回归模型是否存在单位根(unit root)。回归模型可以写为,是一阶差分。测试是否存在单位根等同于测试是否。

因为迪基-福勒检验测试的是残差项,并非原始数据,所以不能用标准t统计量。我们需要用迪基-福勒统计量。

fromstatsmodels.tsa.stattoolsimportadfuller
deftest_stationary(timeseries):
#确定滚动数据
rolmean=timeseries.rolling(24).mean()
rolstd=timeseries.rolling(24).std()

#会议滚动数据
orig=plt.plot(timeseries,color='blue',label='Original')
mean=plt.plot(rolmean,color='red',label='RollingMean')
std=plt.plot(rolstd,color='black',label='RollingStd')
plt.legend(loc='best')
plt.title('RollingMean&StandardDeviation')
plt.show(block=False)

#执行迪基-福勒检验
print('ResultsofDickey-FullerTest:')
dftest=adfuller(timeseries,autolag='AIC')
dfoutput=pd.Series(dftest[0:4],index=['TestStatistic','P-value','#lagsused','NoofObservationsused'])
forkey,valueindftest[4].items():
dfoutput['CriticalValue(%s)'%key]=value
print(dfoutput)

绘制检验图

test_stationary(train_org['Count'])
91ca66be-a3e7-11ec-952b-dac502259ad0.png
Results of Dickey-Fuller Test: 
Test Statistic                -4.456561
P-value                        0.000235
#lags used                    45.000000
No of Observations used    18242.000000
Critical Value (1%)           -3.430709
Critical Value (5%)           -2.861698
Critical Value (10%)          -2.566854
dtype: float64

检验统计数据表明,由于p值小于0.05,数据是平稳的。

移动平均值

在统计学中,移动平均(moving average),又称滑动平均是一种通过创建整个数据集中不同子集的一系列平均数来分析数据点的计算方法。移动平均通常与时间序列数据一起使用,以消除短期波动,突出长期趋势或周期。

对原始数据求对数。

Train_log=np.log(Train['Count'])
valid_log=np.log(Train['Count'])
Train_log.head()
Datetime
2012-08-25    1.152680
2012-08-26    1.299283
2012-08-27    0.949081
2012-08-28    0.882389
2012-08-29    0.916291
Freq: D, Name: Count, dtype: float64

绘制移动平均值曲线

moving_avg=Train_log.rolling(24).mean()
plt.plot(Train_log)
plt.plot(moving_avg,color='red')
plt.show()
91e12674-a3e7-11ec-952b-dac502259ad0.png

去除移动平均值后再进行迪基-福勒检验

train_log_moving_avg_diff=Train_log-moving_avg
train_log_moving_avg_diff.dropna(inplace=True)
test_stationary(train_log_moving_avg_diff)
92011a2e-a3e7-11ec-952b-dac502259ad0.png
Results of Dickey-Fuller Test: 
Test Statistic            -5.861646e+00
P-value                    3.399422e-07
#lags used                 2.000000e+01
No of Observations used    6.250000e+02
Critical Value (1%)       -3.440856e+00
Critical Value (5%)       -2.866175e+00
Critical Value (10%)      -2.569239e+00
dtype: float64

对数时序数据求二阶差分后再迪基-福勒检验

train_log_diff=Train_log-Train_log.shift(1)
test_stationary(train_log_diff.dropna())
921bee62-a3e7-11ec-952b-dac502259ad0.png
Results of Dickey-Fuller Test: 
Test Statistic            -8.237568e+00
P-value                    5.834049e-13
#lags used                 1.900000e+01
No of Observations used    6.480000e+02
Critical Value (1%)       -3.440482e+00
Critical Value (5%)       -2.866011e+00
Critical Value (10%)      -2.569151e+00
dtype: float64

季节性分解

对进行对数转换后的原始数据进行季节性分解。

decomposition=seasonal_decompose(
pd.DataFrame(Train_log).Count.values,freq=24)
trend=decomposition.trend
seasonal=decomposition.seasonal
residual=decomposition.resid

plt.plot(Train_log,label='Original')
plt.plot(trend,label='Trend')
plt.plot(seasonal,label='Seasonality')
plt.plot(residual,label='Residuals')
922acd1a-a3e7-11ec-952b-dac502259ad0.png

对季节性分解后的残差数据进行迪基-福勒检验

train_log_decompose=pd.DataFrame(residual)
train_log_decompose['date']=Train_log.index
train_log_decompose.set_index('date',inplace=True)
train_log_decompose.dropna(inplace=True)
test_stationary(train_log_decompose[0])
924ec422-a3e7-11ec-952b-dac502259ad0.png
Results of Dickey-Fuller Test: 
Test Statistic            -7.822096e+00
P-value                    6.628321e-12
#lags used                 2.000000e+01
No of Observations used    6.240000e+02
Critical Value (1%)       -3.440873e+00
Critical Value (5%)       -2.866183e+00
Critical Value (10%)      -2.569243e+00
dtype: float64

自相关和偏自相关图

fromstatsmodels.tsa.stattoolsimportacf,pacf
lag_acf=acf(train_log_diff.dropna(),nlags=25)
lag_pacf=pacf(train_log_diff.dropna(),nlags=25,method='ols')
plt.plot(lag_acf)
plt.axhline(y=0,linestyle='--',color='gray')
plt.axhline(y=-1.96/np.sqrt(len(train_log_diff.dropna())),linestyle='--',color='gray')
plt.axhline(y=1.96/np.sqrt(len(train_log_diff.dropna())),linestyle='--',color='gray')
plt.plot(lag_pacf)
plt.axhline(y=0,linestyle='--',color='gray')
plt.axhline(y=-1.96/np.sqrt(len(train_log_diff.dropna())),linestyle='--',color='gray')
plt.axhline(y=1.96/np.sqrt(len(train_log_diff.dropna())),linestyle='--',color='gray')
926d5388-a3e7-11ec-952b-dac502259ad0.png929041ea-a3e7-11ec-952b-dac502259ad0.png

AR模型

AR模型训练及预测

model=ARIMA(Train_log,order=(2,1,0))
#这里q值是零,因为它只是AR模型
results_AR=model.fit(disp=-1)
plt.plot(train_log_diff.dropna(),label='original')
plt.plot(results_AR.fittedvalues,color='red',label='predictions')
92a341b4-a3e7-11ec-952b-dac502259ad0.png
AR_predict=results_AR.predict(start="2014-06-25",end="2014-09-25")
AR_predict=AR_predict.cumsum().shift().fillna(0)
AR_predict1=pd.Series(np.ones(valid.shape[0])*np.log(valid['Count'])[0],index=valid.index)
AR_predict1=AR_predict1.add(AR_predict,fill_value=0)
AR_predict=np.exp(AR_predict1)
plt.plot(valid['Count'],label="Valid")
plt.plot(AR_predict,color='red',label="Predict")
plt.legend(loc='best')
plt.title('RMSE:%.4f'%(np.sqrt(np.dot(AR_predict,valid['Count']))/valid.shape[0]))
92ca9c78-a3e7-11ec-952b-dac502259ad0.png

MA模型

model=ARIMA(Train_log,order=(0,1,2))
#这里的p值是零,因为它只是MA模型
results_MA=model.fit(disp=-1)
plt.plot(train_log_diff.dropna(),label='original')
plt.plot(results_MA.fittedvalues,color='red',label='prediction')
92f605d4-a3e7-11ec-952b-dac502259ad0.png
MA_predict=results_MA.predict(start="2014-06-25",end="2014-09-25")
MA_predict=MA_predict.cumsum().shift().fillna(0)
MA_predict1=pd.Series(np.ones(valid.shape[0])*np.log(valid['Count'])[0],index=valid.index)
MA_predict1=MA_predict1.add(MA_predict,fill_value=0)
MA_predict=np.exp(MA_predict1)

plt.plot(valid['Count'],label="Valid")
plt.plot(MA_predict,color='red',label="Predict")
plt.legend(loc='best')
plt.title('RMSE:%.4f'%(np.sqrt(np.dot(MA_predict,valid['Count']))/valid.shape[0]))
9310a146-a3e7-11ec-952b-dac502259ad0.png

ARMA模型

model=ARIMA(Train_log,order=(2,1,2))
results_ARIMA=model.fit(disp=-1)
plt.plot(train_log_diff.dropna(),label='original')
plt.plot(results_ARIMA.fittedvalues,color='red',label='predicted')
9325d228-a3e7-11ec-952b-dac502259ad0.png
defcheck_prediction_diff(predict_diff,given_set):
predict_diff=predict_diff.cumsum().shift().fillna(0)
predict_base=pd.Series(np.ones(given_set.shape[0])*np.log(given_set['Count'])[0],index=given_set.index)
predict_log=predict_base.add(predict_diff,fill_value=0)
predict=np.exp(predict_log)

plt.plot(given_set['Count'],label="Givenset")
plt.plot(predict,color='red',label="Predict")
plt.legend(loc='best')
plt.title('RMSE:%.4f'%(np.sqrt(np.dot(predict,given_set['Count']))/given_set.shape[0]))

defcheck_prediction_log(predict_log,given_set):
predict=np.exp(predict_log)

plt.plot(given_set['Count'],label="Givenset")
plt.plot(predict,color='red',label="Predict")
plt.legend(loc='best')
plt.title('RMSE:%.4f'%(np.sqrt(np.dot(predict,given_set['Count']))/given_set.shape[0]))
plt.show()

ARIMA_predict_diff=results_ARIMA.predict(start="2014-06-25",
end="2014-09-25")
check_prediction_diff(ARIMA_predict_diff,valid)
93542b46-a3e7-11ec-952b-dac502259ad0.png

SARIMAX模型

y_hat_avg=valid.copy()
fit1=sm.tsa.statespace.SARIMAX(Train.Count,order=(2,1,4),seasonal_order=(0,1,1,7)).fit()
y_hat_avg['SARIMA']=fit1.predict(start="2014-6-25",end="2014-9-25",dynamic=True)
plt.plot(Train['Count'],label='Train')
plt.plot(valid['Count'],label='Valid')
plt.plot(y_hat_avg['SARIMA'],label='SARIMA')
8e6ef822-a3e7-11ec-952b-dac502259ad0.png

模型评价

rms=sqrt(mean_squared_error(valid.Count,y_hat_avg.SARIMA))
print(rms)
70.26240839723575

预测

predict=fit1.predict(start="2014-9-26",end="2015-4-26",dynamic=True)
test['prediction']=predict
#按日、月、年合并Test和test_original
merge=pd.merge(test,test_org,on=('day','Month','Year'),how='left')
merge['Hour']=merge['Hour_y']
merge=merge.drop(['Year','Month','Datetime','Hour_x','Hour_y'],axis=1
#通过合并merge和temp2进行预测
prediction=pd.merge(merge,temp2,on='Hour',how='left')
#将比率转换成原始比例
prediction['Count']=prediction['prediction']*prediction['ratio']*24
prediction['ID']=prediction['ID_y']
submission=prediction.drop(['day','Hour','ratio','prediction','ID_x','ID_y'],axis=1)
#转换最终提交的csv格式
pd.DataFrame(submission,columns=['ID','Count']).to_csv('SARIMAX.csv')

939ad33e-a3e7-11ec-952b-dac502259ad0.png

原文标题:时间序列分析和预测实战

文章出处:【微信公众号:数据分析与开发】欢迎添加关注!文章转载请注明出处。

审核编辑:彭菁

声明:本文内容及配图由入驻作者撰写或者入驻合作网站授权转载。文章观点仅代表作者本人,不代表电子发烧友网立场。文章及其配图仅供工程师学习之用,如有内容侵权或者其他违规问题,请联系本站处理。 举报投诉
  • 数据
    +关注

    关注

    8

    文章

    6504

    浏览量

    87444
  • 函数
    +关注

    关注

    3

    文章

    3845

    浏览量

    61223
  • 时间序列
    +关注

    关注

    0

    文章

    30

    浏览量

    10360

原文标题:时间序列分析和预测实战

文章出处:【微信号:DBDevs,微信公众号:数据分析与开发】欢迎添加关注!文章转载请注明出处。

收藏 人收藏

    评论

    相关推荐

    波导理论知识

    波导理论知识 学习资料! 
    发表于 11-03 08:36

    测试理论知识-美国国家仪器内部AE培训资料

    附件是美国国家仪器应用工程师内部培训资料,主要针对测试测量基础理论知识
    发表于 10-21 11:01

    基础理论之反射详解

    ,走线从细变宽,会增加一次反射,那是不是全程按照breakout区域走线会比较好?源端匹配电阻是不是也增加了一次反射? ......小编在此给大家分享下信号反射的基础理论知识,希望对大家有用。
    发表于 06-15 17:07

    检测技术基础理论知识

    检测技术基础理论知识
    发表于 11-27 14:51

    你都知道单片机的基础理论知识学习包括哪些吗

    你都知道单片机的基础理论知识学习包括哪些吗?
    发表于 01-21 07:20

    串行和并行通讯的基础理论知识分析

    基础理论知识串行和并行通讯比喻为单车道和多车道全双工、半双工及单工通讯通讯方式说明全双工在同一时刻,两个设备之间可以同时收发数据半双工两个设备之间可以收发数据,但不能在同一时刻进行单工在任何时刻都
    发表于 02-22 06:32

    科学数据时间序列预测方法

    针对传统的时间序列分析方法预测科学数据效果较差的特点,提出了一种结合自组织神经网络和灰色理论时间
    发表于 12-31 23:56 10次下载

    改进GP分形理论的最近邻序列预测算方法

    改进GP分形理论的最近邻序列预测算方法:针对现有的时间序列分析
    发表于 01-03 17:00 12次下载

    USB接口的基础理论知识

    USB接口的基础理论知识 USB的重要关键字:
    发表于 10-15 13:08 819次阅读
    USB接口的<b class='flag-5'>基础理论知识</b>

    开关电源的理论知识

    开关电源的理论知识 - 张兴柱博士,下来看看。
    发表于 03-31 10:26 0次下载

    维修电工中级理论知识试卷分享

    本文主要内容是维修电工的中级理论知识
    发表于 10-11 15:52 24次下载
    维修电工中级<b class='flag-5'>理论知识</b>试卷分享

    电信传输理论知识详解

    本文详细介绍了电信传输理论知识(新编)。
    发表于 11-16 15:57 3次下载
    电信传输<b class='flag-5'>理论知识</b>详解

    变频器的故障分析和解决 实践检验、理论知识及维修水平

    变频器的故障分析和解决不仅需要实践检验,还需要掌握相关变频器相关的理论知识理论分析与维修实践相结合以达到较高的维修水平,维修人员需要学习的相关理论
    的头像 发表于 07-17 08:54 3894次阅读
    变频器的故障<b class='flag-5'>分析</b>和解决 实践检验、<b class='flag-5'>理论知识</b>及维修水平

    单片机学习笔记:基础理论知识学习

    基础理论知识包括模拟电路、数字电路和C语言知识。模拟电路和数字电路属于抽象学科,要把它学好还得费点精神。在你学习单片机之前,觉得模拟电路和数字电路基础不好的话,不要急着学习单片机,应该先回顾所学过的模拟电路和数字电路知识,为学习
    的头像 发表于 05-30 10:45 3815次阅读

    详解单片机基础理论知识

    单片机基础理论知识包括模拟电路、数字电路和C语言知识,模拟电路与数字电路属于抽象学科,要把这些学好还得费点功夫,但是在学习单片机之前,觉得自己模拟电路和数字电路基础不好的话,先不要急着入门单片机
    发表于 03-30 17:19 4840次阅读
    详解单片机<b class='flag-5'>基础理论知识</b>