python时间序列分析(ARIMA模型)

python时间序列分析(ARIMA模型)
转载请注明出处。
什么是时间序列
时间序列简单的说就是各时间点上形成的数值序列,时间序列分析就是通过观察历史数据预测未来的值。在这⾥需要强调⼀点的是,时间序列分析并不是关于时间的回归,它主要是研究⾃⾝的变化规律的(这⾥不考虑含外⽣变量的时间序列)。
为什么⽤python
  ⽤两个字总结“情怀”,爱屋及乌,个⼈⽐较喜欢python,就⽤python撸了。能做时间序列的软件很多,SAS、R、SPSS、Eviews甚⾄matlab等等,实际⼯作中应⽤得⽐较多的应该还是SAS和R,前者推荐王燕写的《应⽤时间序列分析》,后者推荐“”这篇博⽂()。python作为科学计算的利器,当然也有相关分析的包:statsmodels中tsa模块,当然这个包和SAS、R是⽐不了,但是python有另⼀个神器:pandas!pandas在时间序列上的应⽤,能简化我们很多的⼯作。
环境配置
  python推荐直接装Anaconda,它集成了许多科学计算包,有⼀些包⾃⼰⼿动去装还是挺费劲的。statsmodels需要⾃⼰去安装,这⾥我推荐使⽤0.6的稳定版,0.7及其以上的版本能在github上到,该版本在安装时会⽤C编译好,所以修改底层的⼀些代码将不会起作⽤。
时间序列分析
1.基本模型
  ⾃回归移动平均模型(ARMA(p,q))是时间序列中最为重要的模型之⼀,它主要由两部分组成: AR代表p阶⾃回归过程,MA代表q阶移动平均过程,其公式如下:
  依据模型的形式、特性及⾃相关和偏⾃相关函数的特征,总结如下:
在时间序列中,ARIMA模型是在ARMA模型的基础上多了差分的操作。
2.pandas时间序列操作
⼤熊猫真的很可爱,这⾥简单介绍⼀下它在时间序列上的可爱之处。和许多时间序列分析⼀样,本⽂同样使⽤航空乘客数据(AirPassengers.csv)作为样例。
数据读取:
# -*- coding:utf-8 -*-
import numpy as np
import pandas as pd
from datetime import datetime
import matplotlib.pylab as plt
# 读取数据,pd.read_csv默认⽣成DataFrame对象,需将其转换成Series对象
虹膜定位df = pd.read_csv('AirPassengers.csv', encoding='utf-8', index_col='date')
df.index = pd.to_datetime(df.index)  # 将字符串索引转换成时间索引
ts = df['x']  # ⽣成pd.Series对象
# 查看数据格式
ts.head()
ts.head().index
查看某⽇的值既可以使⽤字符串作为索引,⼜可以直接使⽤时间对象作为索引
ts['1949-01-01']
ts[datetime(1949,1,1)]
两者的返回值都是第⼀个序列值:112
如果要查看某⼀年的数据,pandas也能⾮常⽅便的实现
ts['1949']
切⽚操作:
ts['1949-1' : '1949-6']
注意时间索引的切⽚操作起点和尾部都是包含的,这点与数值索引有所不同
pandas还有很多⽅便的时间序列函数,在后⾯的实际应⽤中在进⾏说明。
3. 平稳性检验
我们知道序列平稳性是进⾏时间序列分析的前提条件,很多⼈都会有疑问,为什么要满⾜平稳性的要求呢?在⼤数定理和中⼼定理中要求样本同分布(这⾥同分布等价于时间序列中的平稳性),⽽我们的建模过程中有很多都是建⽴在⼤数定理和中⼼极限定理的前提条件下的,如果它不满⾜,得到的许多结论都是不可靠的。以虚假回归为例,当响应变量和输⼊变量都平稳时,我们⽤t统计量检验标准化系数的显著性。⽽当响应变量和输⼊变量不平稳时,其标准化系数不在满⾜t分布,这时再⽤t检验来进⾏显著性分析,导致拒绝原假设的概率增加,即容易犯第⼀类错误,从⽽得出错误的结论。
平稳时间序列有两种定义:严平稳和宽平稳
严平稳顾名思义,是⼀种条件⾮常苛刻的平稳性,它要求序列随着时间的推移,其统计性质保持不变。对于任意的τ,其联合概率密度函数满⾜:
严平稳的条件只是理论上的存在,现实中⽤得⽐较多的是宽平稳的条件。
宽平稳也叫弱平稳或者⼆阶平稳(均值和⽅差平稳),它应满⾜:
常数均值
常数⽅差
跳线帽
常数⾃协⽅差
平稳性检验:观察法和单位根检验法
基于此,我写了⼀个名为test_stationarity的统计性检验模块,以便将某些统计检验结果更加直观的展现出来。
# -*- coding:utf-8 -*-
from statsmodels.tsa.stattools import adfuller
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
aphics.tsaplots import plot_acf, plot_pacf
# 移动平均图
def draw_trend(timeSeries, size):
f = plt.figure(facecolor='white')
# 对size个数据进⾏移动平均
rol_mean = lling(window=size).mean()
# 对size个数据进⾏加权移动平均
rol_weighted_mean = pd.ewma(timeSeries, span=size)
timeSeries.plot(color='blue', label='Original')
rolmean.plot(color='red', label='Rolling Mean')
rol_weighted_mean.plot(color='black', label='Weighted Rolling Mean')
plt.legend(loc='best')
plt.title('Rolling Mean')
plt.show()
def draw_ts(timeSeries):
f = plt.figure(facecolor='white')
timeSeries.plot(color='blue')
plt.show()
'''
  Unit Root Test
The null hypothesis of the Augmented Dickey-Fuller is that there is a unit
root, with the alternative that there is no unit root. That is to say the
bigger the p-value the more reason we assert that there is a unit root
'''
方形磁铁def testStationarity(ts):
dftest = adfuller(ts)
# 对上述函数求得的值进⾏语义描述
dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
for key,value in dftest[4].items():
ttx2基板dfoutput['Critical Value (%s)'%key] = value
return dfoutput
# ⾃相关和偏相关图,默认阶数为31阶
def draw_acf_pacf(ts, lags=31):
楔形塞尺
f = plt.figure(facecolor='white')
ax1 = f.add_subplot(211)
plot_acf(ts, lags=31, ax=ax1)
ax2 = f.add_subplot(212)
plot_pacf(ts, lags=31, ax=ax2)
plt.show()
观察法,通俗的说就是通过观察序列的趋势图与相关图是否随着时间的变化呈现出某种规律。所谓的规律就是时间序列经常提到的周期性因素,现实中遇到得⽐较多的是线性周期成分,这类周期成分可以采⽤差分或者移动平均来解决,⽽对于⾮线性周期成分的处理相对⽐较复杂,需要采⽤某些分解的⽅法。下图为航空数据的线性图,可以明显的看出它具有年周期成分和长期趋势成分。平稳序列的⾃相关系数会快速衰减,下⾯的⾃相关图并不能体现出该特征,所以我们有理由相信该序列是不平稳的。
单位根检验:ADF是⼀种常⽤的单位根检验⽅法,他的原假设为序列具有单位根,即⾮平稳,对于⼀个平稳的时序数据,就需要在给定的置信⽔平上显著,拒绝原假设。ADF只是单位根检验的⽅法之⼀,如果想采⽤其他检验⽅法,可以安装第三⽅包arch,⾥⾯提供了更加全⾯的单位根检验⽅法,个⼈还是⽐较钟情ADF检验。以下为检验结果,其p值⼤于0.99,说明并不能拒绝原假设。
3. 平稳性处理
由前⾯的分析可知,该序列是不平稳的,然⽽平稳性是时间序列分析的前提条件,故我们需要对不平稳的序列进⾏处理将其转换成平稳的序列。a. 对数变换
对数变换主要是为了减⼩数据的振动幅度,使其线性规律更加明显(我是这么理解的时间序列模型⼤部分都是线性的,为了尽量降低⾮线性的因素,需要对其进⾏预处理,也许我理解的不对)。对数变换相当于增加了⼀个惩罚机制,数据越⼤其惩罚越⼤,数据越⼩惩罚越⼩。这⾥强调⼀下,变换的序列需要满⾜⼤于0,⼩于0的数据不存在对数变换。
ts_log = np.log(ts)
test_stationarity.draw_ts(ts_log)
b. 平滑法
根据平滑技术的不同,平滑法具体分为移动平均法和指数平均法。
移动平均即利⽤⼀定时间间隔内的平均值作为某⼀期的估计值,⽽指数平均则是⽤变权的⽅法来计算均值
test_stationarity.draw_trend(ts_log, 12)
从上图可以发现窗⼝为12的移动平均能较好的剔除年周期性因素,⽽指数平均法是对周期内的数据进⾏了加权,能在⼀定程度上减⼩年周期因素,但并不能完全剔除,如要完全剔除可以进⼀步进⾏差分操作。
c.  差分
时间序列最常⽤来剔除周期性因素的⽅法当属差分了,它主要是对等周期间隔的数据进⾏线性求减。前⾯我们说过,ARIMA模型相对ARMA模型,仅多了差分操作,ARIMA模型⼏乎是所有时间序列软件都⽀持的,差分的实现与还原都⾮常⽅便。⽽statsmodel中,对差分的⽀持不是很好,它不⽀持⾼阶和多阶差分,为什么不⽀持,这⾥引⽤作者的说法:
作者⼤概的意思是说预测⽅法中并没有解决⾼于2阶的差分,有没有感觉很牵强,不过没关系,我们有pandas。我们可以先⽤pandas将序列差分好,然后在对差分好的序列进⾏ARIMA拟合,只不过这样后⾯会多了⼀步⼈⼯还原的⼯作。
快速插头
diff_12 = ts_log.diff(12)
diff_12.dropna(inplace=True)
diff_12_1 = diff_12.diff(1)
diff_12_1.dropna(inplace=True)
stStationarity(diff_12_1)
从上⾯的统计检验结果可以看出,经过12阶差分和1阶差分后,该序列满⾜平稳性的要求了。
d. 分解
所谓分解就是将时序数据分离成不同的成分。statsmodels使⽤的X-11分解过程,它主要将时序数据分离成长期趋势、季节趋势和随机成分。与其它统计软件⼀样,statsmodels也⽀持两类分解模型,加法模型和乘法模型,这⾥我只实现加法,乘法只需将model的参数设置为"multiplicative"即可。
from statsmodels.tsa.seasonal import seasonal_decompose
decomposition = seasonal_decompose(ts_log, model="additive")
trend = d
seasonal = decomposition.seasonal
residual = sid
得到不同的分解成分后,就可以使⽤时间序列模型对各个成分进⾏拟合,当然也可以选择其他预测⽅法。我曾经⽤过⼩波对时序数据进⾏过分解,然后分别采⽤时间序列拟合,效果还不错。由于我对⼩波的理解不是很好,只能简单的调⽤接⼝,如果有谁对⼩波、傅⾥叶、卡尔曼理解得⽐较透,可以将时序数据进⾏更加准确的分解,由于分解后的时序数据避免了他们在建模时的交叉影响,所以我相信它将有助于预测准确性的提⾼。
4. 模型识别
在前⾯的分析可知,该序列具有明显的年周期与长期成分。对于年周期成分我们使⽤窗⼝为12的移动平进⾏处理,对于长期趋势成分我们采⽤1阶差分来进⾏处理。
rol_mean = lling(window=12).mean()
rol_mean.dropna(inplace=True)
ts_diff_1 = rol_mean.diff(1)
ts_diff_1.dropna(inplace=True)
stStationarity(ts_diff_1)
观察其统计量发现该序列在置信⽔平为95%的区间下并不显著,我们对其进⾏再次⼀阶差分。再次差分后的序列其⾃相关具有快速衰减的特点,t统计量在99%的置信⽔平下是显著的,这⾥我不再做详细说明。
ts_diff_2 = ts_diff_1.diff(1)
ts_diff_2.dropna(inplace=True)
数据平稳后,需要对模型定阶,即确定p、q的阶数。观察上图,发现⾃相关和偏相系数都存在拖尾的特点,并且他们都具有明显的⼀阶相关性,所以我们设定p=1, q=1。下⾯就可以使⽤ARMA模型进⾏数据拟合了。这⾥我不使⽤ARIMA(ts_diff_1, order=(1, 1, 1))进⾏拟合,是因为含有差分操作时,预测结果还原⽼出问题,⾄今还没弄明⽩。
from statsmodels.tsa.arima_model import ARMA
model = ARMA(ts_diff_2, order=(1, 1))
result_arma = model.fit( disp=-1, method='css')
5. 样本拟合
模型拟合完后,我们就可以对其进⾏预测了。由于ARMA拟合的是经过相关预处理后的数据,故其预测值需要通过相关逆变换进⾏还原。
predict_ts = result_arma.predict()
# ⼀阶差分还原
diff_shift_ts = ts_diff_1.shift(1)
diff_recover_1 = predict_ts.add(diff_shift_ts)
# 再次⼀阶差分还原

本文发布于:2024-09-21 19:25:46,感谢您对本站的认可!

本文链接:https://www.17tex.com/tex/3/276664.html

版权声明:本站内容均来自互联网,仅供演示用,请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系,我们将在24小时内删除。

标签:序列   时间   数据
留言与评论(共有 0 条评论)
   
验证码:
Copyright ©2019-2024 Comsenz Inc.Powered by © 易纺专利技术学习网 豫ICP备2022007602号 豫公网安备41160202000603 站长QQ:729038198 关于我们 投诉建议