大家好，又见面了，我是你们的朋友全栈君。
我于2019年发布此篇文章至今收获了许多人的指点，当时的代码的确晦涩难懂，近期有空，将代码重新整理了一遍，重新发送至此。希望能够帮助大家更好地理解。
建模步骤：
目录
数据包和版本申明
步骤一：数据准备与数据预处理
步骤二：数据重采样
步骤三：平滑处理
步骤四：平稳性检验
步骤五： 时间序列定阶
（2）信息准则定阶
步骤六：模型构建
步骤七：模型评价
总结
数据包和版本申明
申明：本实验环境为python 3.7.4 statsmodels版文为：0.10.1
代码语言：javascript复制import pandas as pd
import numpy as np
import seaborn as sns #热力图
import itertools
import datetime
import matplotlib.pyplot as plt #画图
import statsmodels.api as sm
from sklearn.metrics import mean_squared_error
from statsmodels.tsa.stattools import adfuller #ADF检验
from statsmodels.stats.diagnostic import acorr_ljungbox #白噪声检验
from statsmodels.graphics.tsaplots import plot_acf,plot_pacf #画图定阶
from statsmodels.tsa.arima_model import ARIMA #ARIMA模型
from statsmodels.tsa.arima_model import ARMA #ARMA模型
from statsmodels.stats.stattools import durbin_watson #DW检验
from statsmodels.graphics.api import qqplot #qq图步骤一：数据准备与数据预处理
自动生成2018年1月1日至2018年9月1日数据，数据服从标准正态分布，存入old_data.csv中。
代码语言：javascript复制#### Part:generate raw data and save in the old_data.csv
#### 创建一个时间列表，从20180101到20180901数据，存入 old_data.csv
def genertate_data():
index = pd.date_range(start='2018-1-1',end = '2018-9-1',freq='10T') # 10分钟采样一次
index = list(index)
data_list = []
for i in range(len(index)):
data_list.append(np.random.randn()) # 数据是符合标准正态分布的样本
dataframe = pd.DataFrame({'time':index,'values':data_list})
dataframe.to_csv('G:\WX\2\old_data.csv',index=0)
print('the data is existting')故意去将文件夹中的某些值，改成了-10000，弄成了异常值，（因为老师说尽可能显得步骤完整，最后分数才会高-，-所以我自己手动添加异常）。这块的主要工作就是利用pandas里面的函数，去查看一下刚特殊操作后的数据。
代码语言：javascript复制#### Step 1 数据预处理
#### delete or revise some values in data and make data preprocessing
#### 删掉或者修改创建的数据后，进行简单数据预处理
def data_preprocessing():
data = pd.read_csv('G:\WX\2\old_data.csv')
#print(data.describe()) #查看统计信息,发现最小值有-10000的异常数据
#print((data.isnull()).sum()) #查看是否存在缺失值
#print((data.duplicated()).sum()) #重复值
def change_zero(x):
if x == -10000:
return 0
else :
return x
data['values'] = data['values'].apply(lambda x: change_zero(x))
#利用均值填充缺失值
mean = data['values'].mean()
def change_mean(x):
if x == 0:
return mean
else:
return x
data['values'] = data['values'].apply(lambda x: change_mean(x))
#保存处理过的数据
data.to_csv('G:\WX\2\new_data.csv',index=0)
print('new data is existing')步骤二：数据重采样
为了得高分（-，-），做了很多个数据，然后一共有34992个数据，然后进行了一下重采样，数据以天进行重采样。
代码语言：javascript复制#### Step 2 重采样
#### Resample Data and Sampling frequency is days
#### 重采样，将采样频率换成以天为单位
def Resampling(): #重采样
df = pd.read_csv('G:\WX\2\new_data.csv')
#将默认索引方式转换成时间索引
df['time'] = pd.to_datetime(df['time'])
df.set_index("time", inplace=True)
train_data = df['2018-1-1':'2018-8-1'] ## 取到20180101 至 20180801 做训练
test = df['2018-8-1':'2018-9-1'] ## 取到20180801 至 20180901 做预测
train_data = train_data.resample('D').mean() ## 以天为时间间隔取均值,重采样
test_data = test.resample('D').mean()
return train_data,test_data步骤三：平滑处理
由于ARMA和ARIMA需要时间序列满足平稳性和非白噪声的要求，所以要用差分法和平滑法（滚动平均和滚动标准差）来实现序列的平稳性操作。一般情况下，对时间序列进行一阶差分法就可以实现序列的平稳性，有时需要二阶差分。
代码语言：javascript复制#### Step 3 差分转平稳
def stationarity(timeseries): #平稳性处理（timeseries 时间序列）
## 差分法,保存成新的列
diff1 = timeseries.diff(1).dropna() # 1阶差分 dropna() 删除缺失值
diff2 = diff1.diff(1) #在一阶差分基础上再做一次一阶差分，即二阶查分
## 画图
#diff1.plot(color = 'red',title='diff 1',figsize=(10,4))
#diff2.plot(color = 'black',title='diff 2',figsize=(10,4))
## 平滑法
rollmean = timeseries.rolling(window=4,center = False).mean() ## 滚动平均
rollstd = timeseries.rolling(window=4,center = False).std() ## 滚动标准差
## 画图
#rollmean.plot(color = 'yellow',title='Rolling Mean',figsize=(10,4))
#rollstd.plot(color = 'blue',title='Rolling Std',figsize=(10,4))
return diff1,diff2,rollmean,rollstd差分法处理结果图：由图可以看出 一阶差分基本就满足了平稳性需要。
平滑法处理结果如图所示。
可以看出，平滑法不太适合我造出来的数据。一般情况下，平滑法更适合带有周期性稳步上升的数据类型。
步骤四：平稳性检验
利用ADF检验判断序列是否平稳，利用白噪声检验判断序列是否为随机性序列。
代码语言：javascript复制#### Step 4 平稳性检验
def ADF_test(timeseries): ## 用于检测序列是否平稳
x = np.array(timeseries['values'])
adftest = adfuller(x, autolag='AIC')
#print (adftest)
if adftest[0] < adftest[4]["1%"] and adftest[1] < 10**(-8):
# 对比Adf结果和10%的时的假设检验 以及 P-value是否非常接近0(越小越好)
print("序列平稳")
return True
else:
print("非平稳序列")
return False
def random_test(timeseries) : #随机性检验（白噪声检验）
p_value = acorr_ljungbox(timeseries, lags=1) # p_value 返回二维数组，第二维为P值
if p_value[1] < 0.05:
print("非随机性序列")
return True
else:
print("随机性序列,即白噪声序列")
return False（1）ADF检验结果如下：
如何确定该序列能否平稳呢？主要看：
（1）1%、%5、不同程度拒绝原假设的统计值和ADF Test result的比较，ADF Test result同时小于1%、5%、10%即说明非常好地拒绝该假设，本数据中，adf结果为-6.9， 小于三个level的统计值。
（2）P-value是否非常接近0.本数据中，P-value 为 7.9e-10，接近0。
ADF结果如何查看参考了这篇博客：
Python时间序列中ADF检验详解_学渣渣-CSDN博客_python进行adf检验
（2）白噪声结果如图：
统计量的P值小于显著性水平0.05，则可以以95%的置信水平拒绝原假设，认为序列为非白噪声序列（否则，接受原假设，认为序列为纯随机序列。）
由于P值为0.315远大于0.05所以接受原假设，认为时间序列是白噪声的，即是随机产生的序列，不具有时间上的相关性。（解释一下，由于老师没有给数据，所以只能硬着头皮，假设它是非白噪声的做）
步骤五： 时间序列定阶
定阶方法主要为两种：
（1）ACF和PACF 利用拖尾和截尾来确定
（2）信息准则定阶（AIC、BIC、HQIC）
热力图为辅助方法使得p与q的取值更加明确,可视化,实际上依旧为方法二-信息准则定阶。
代码语言：javascript复制def determinate_order_acf(timeseries): #利用ACF和PACF判断模型阶数
plot_acf(timeseries,lags=40) #延迟数
plot_pacf(timeseries,lags=40)
plt.show()
def detetminante_order_AIC(timeseries): #信息准则定阶：AIC、BIC、HQIC
#AIC
AIC = sm.tsa.arma_order_select_ic(timeseries,
max_ar=6,max_ma=4,ic='aic')['aic_min_order']
#BIC
BIC = sm.tsa.arma_order_select_ic(timeseries,max_ar=6,
max_ma=4,ic='bic')['bic_min_order']
#HQIC
HQIC = sm.tsa.arma_order_select_ic(timeseries,max_ar=6,
max_ma=4,ic='hqic')['hqic_min_order']
print('the AIC is{},nthe BIC is{}n the HQIC is{}'.format(AIC,BIC,HQIC))
def heatmap_AIC(timeseries):
#设置遍历循环的初始条件，以热力图的形式展示，原理同AIC，BIC，HQIC定阶
p_min = 0
q_min = 0
p_max = 5
q_max = 5
d_min = 0
d_max = 5
# 创建Dataframe,以BIC准则
results_aic = pd.DataFrame(index=['AR{}'.format(i)
for i in range(p_min,p_max 1)],
columns=['MA{}'.format(i) for i in range(q_min,q_max 1)])
# itertools.product 返回p,q中的元素的笛卡尔积的元组
for p,d,q in itertools.product(range(p_min,p_max 1),
range(d_min,d_max 1),range(q_min,q_max 1)):
if p==0 and q==0:
results_aic.loc['AR{}'.format(p), 'MA{}'.format(q)] = np.nan
continue
try:
model = sm.tsa.ARIMA(timeseries, order=(p, d, q))
results = model.fit()
#返回不同pq下的model的BIC值
results_aic.loc['AR{}'.format(p), 'MA{}'.format(q)] = results.aic
except:
continue
results_aic = results_aic[results_aic.columns].astype(float)
#print(results_bic)
fig, ax = plt.subplots(figsize=(10, 8))
ax = sns.heatmap(results_aic,
#mask=results_aic.isnull(),
ax=ax,
annot=True, #将数字显示在热力图上
fmt='.2f',
)
ax.set_title('AIC')
plt.show() 直接利用步骤3的一阶差分来进行定阶，结果如图所示：
上面分别是ACF和PACF的图，至于如何定阶不详细叙述了。一般是通过截尾和拖尾来确定阶数。
热力图定阶结果如下所示：
黑色的位置最好，可以看出p，q取（1,1）（3,1）（1,4）都可以。一般情况下是越小越好。
热力图实现过程参考了了一篇文章，但是博客链接丢失，如果有侵权，请告知，会删除相关部分。
步骤六：模型构建
代码语言：javascript复制def ARMA_model(train_data,order): # 训练数据,测试数据,定阶
arma_model = ARMA(train_data,order) #ARMA模型
arma = arma_model.fit()#激活模型
#print(result.summary()) #给出一份模型报告
############ in-sample ############ 样本内预测
in_sample_pred = arma.predict()
############ out-sample ########## 样本外预测
#### 样本外预测需要从train_data 样本内的某一个时间节点开始
#### 利用start和end 控制样外预测 起止时间
out_sample_pred = arma.predict(start=len(train_data)-2,end = len(train_data) 30,
dynamic=True)
#in_sample_pred.plot()
#train_data.plot()
return arma,in_sample_pred,out_sample_pred
def ARIMA_model(train_data,order):
arima_model = ARIMA(train_data,order) #ARIMA模型
arima = arima_model.fit()
#print(result.summary()) #给出一份模型报告
########样本内预测#########
in_sample_pred = arima.predict()
####### 样本外预测##########
out_sample_pred = arima.predict(start=len(train_data)-2,end = len(train_data) 30,
dynamic=True)
return arima,in_sample_pred,out_sample_pred预测过程有两种预测方式，一种是样本内的预测（in_sample_pred），一种是样本外的预测(out_sample_pred)。样本内预测就是的是2018-1-1到2018-8-1的。但是要预测的是8-1到9-1的情况，是out-sample预测，一般情况下，out-sample是我们想要的，而不是样本内的预测。
样本外预测是由dynamic参数决定的，特别注意：样本外的预测也要从样本内的某一个时间点开始才能进行预测。因此样本外的预测开始时间要从train_data长度内的某一个时间节点开始。
步骤七：模型评价
主要分为四种方法：（1）QQ图检验残差是否满足正态分布（2）利用D-W检验,检验残差的自相关性（3）计算预测值和真实值的标准差，误差相关等 （4）还原预测序列和测试序列，用图来直观评价模型
代码语言：javascript复制def evaluate_model(model,train_data,predict_data):
###（1）利用QQ图检验残差是否满足正态分布
resid = model.resid # 求解模型残差
plt.figure(figsize=(12,8))
qqplot(resid,line='q',fit=True)
###（2）利用D-W检验,检验残差的自相关性
print('D-W检验值为{}'.format(durbin_watson(resid.values)))
###（3）利用预测值和真实值的误差检测，这里用的是标准差
#row_train_data 是从 2018-1-1开始的,经过差分后train_data发生变化
print('标准差为{}'.format(mean_squared_error(train_data,predict_data,sample_weight=None,
multioutput='uniform_average'))) #标准差（均方差）
def string_toDatetime(string): # 截取时间
return datetime.datetime.strptime(string, "%Y-%m-%d %H:%M:%S")
#### 绘制图像,查看预测效果
def draw_picture(row_train_data,out_sample_pred,test_data):
#print(out_sample_pred)
# 样本外预测传入 test_data,out_sample_pred
# 由于预测都是由差分后的平稳序列得出,因此需要对差分后的数据进行还原
# 还原后绘制同一起点的曲线
#######还原 out_sample_pred #########
#### out_sample 2018-07-31
#### test_data 2018-8-1
##2018-8-1 00:00 到 2018-9-1 00:00 ###
#将差分后的序列还原,re_out_sample_pred为还原之后
re_out_sample_pred = pd.Series(np.array(row_train_data)[-2][0],
index=[row_train_data.index[-2]]).append(out_sample_pred[1:]).cumsum()
#### 横坐标
x = []
for i in range(32):
x.append(i 1)
x = np.array(x)
#### 纵坐标
y1 = np.array(test_data)
y2 = np.array(re_out_sample_pred[1:])
#### 画图
plt.plot(x,y1,color='blue')
plt.plot(x,y2,color='red')
plt.show()（1）qq图如下所示：
通过qq图可以看出，残差基本满足了正态分布。
（2）D-W检验结果为：
当D-W检验值接近于2时，不存在自相关性，说明模型较好。
D-W检验如何数学说明，可以参考下面链接。
DW值判断准则 – 百度文库
（3）利用标准差来评价模型时，尤其为样本外预测时，注意时间序列的时间对齐。
在利用图来还原预测数据的过程中，主要利用cumsum（）函数，主要作用是累加操作。
代码语言：javascript复制 re_out_sample_pred = pd.Series(np.array(row_train_data)[-2][0],
index=[row_train_data.index[-2]]).append(out_sample_pred[1:]).cumsum()调用以上步骤函数代码如下：
代码语言：javascript复制if __name__ == "__main__":
## Step 1 and Step 2 都只运行一遍
genertate_data() # 生成数据
data_preprocessing() # 1：数据预处理
train_data,test_data = Resampling() #:2：数据重采样,返回训练数据和测试数据
row_train_data = train_data # 保存差分前的序列,为了后面做评估
### 差分或者滚动平均,利用步骤4函数确定
Smooth_data = stationarity(train_data) # 4 差分
for data in zip(Smooth_data,range(4)):# range(4) 用于判断哪种方法 满足平稳性和白噪声
if ADF_test(data[0]) and random_test(data[0]) : # 平稳性和白噪声检测
train_data = data[0] # 先用差分，再用平滑,分别对应4个序列
method = data[1]
print(method) #### 如果是差分做的,那么后面ARIMA模型中要使用这个参数
break
## 三种选择一种即可
determinate_order_acf(train_data) # ACF定阶
detetminante_order_AIC(train_data) # BIC 定阶
heatmap_AIC(train_data) # 热力图 显示
#### 模型建议和模型评价
#### order 由差分和定阶给出
order = (1,1) ## ARMA p,q
order = (2,1,0) ## ARIMA p,d,q
#### 调用模型
arma,in_sample_pred,out_sample_pred = ARMA_model(train_data,order)
arima,in_sample_pred,out_sample_pred = ARIMA_model(train_data,order)
#### 模型评价（样本内外均可,此处只用于样本内）
evaluate_model(arma,train_data,in_sample_pred) # 样本内预测
evaluate_model(arima,train_data,in_sample_pred) # 样本内预测
#### 画图-差异比较（样本外预测）
draw_picture(row_train_data,out_sample_pred,test_data)总结
关于ARMA和ARIMA模型，从数据处理到最后建模实现，就完成了。
但是，里面其实有一个很大的问题，就是当数据不是平稳性的数据的时候，用到了差分法进行处理，用到了dropna()这个函数，这个函数的意思是去掉序列中nan（在这个了里面是0）。因此当序列中两列相邻值相等时，就会去掉前面那一列，因此处理后的数据可能不是按照每一天的数据分布的，但是预测出来的是每一天都存在的。
如果不加dropna（）这个函数的话，定阶那些都会报错（错误信息是存在nan值），但是模型不会报错。因此这块是一个存在的问题，还亟待处理，但是我看了很多的文章，对于这里好像没有过多深入的研究。
整篇博客都是代码实现的，具体的数学公式，自行百度吧~~
2022年1月5日更新：
此版本代码基于statsmodels 0.10.1实现，目前最新statsmodels版本，已经不能像下面这样调用ARMA 和ARIMA模型了。
代码语言：javascript复制from statsmodels.tsa.arima_model import ARIMA #ARIMA模型
from statsmodels.tsa.arima_model import ARMA #ARMA模型新版本为：
代码语言：javascript复制from statsmodels.tsa.arima.model import ARIMA #ARIMA模型关于ARMA好像封装到了 其他模块中，如果想要利用ARMA模型去讨论，还需去官网自行查看。
由于文章过长，完整版代码，会同步到下面公众号，支持提问并第一时间回复。欢迎关注。
发布者：全栈程序员栈长，转载请注明出处：https://javaforall.cn/142302.html原文链接：https://javaforall.cn
