实验

			
名称
			
			

			
疏系数模型 和季节模型
			
		

			
实验

			
内容
			
			

			
1、简单季节模型
			
		

			
实验

			
目的
			
			

			
1、掌握疏系数模型

			
2、熟练建立季节模型
			
		
 

目录

简单季节模型结构

模型建立

时序图

差分平稳化

白噪声检验

模型定阶

参数估计和模型检验

模型预测

推荐阅读

使用Python完成时间序列分析基础
	
SPSS建立时间序列乘法季节模型实战案例
	
Python建立时间序列ARIMA模型实战案例
简单季节模型结构



模型建立

时序图



时序图显示，该序列既包含长期趋势又包含以年为周期的季节效应

 

差分平稳化

对原序列做1阶差分消去趋势，再做4步差分消去季节效应的影响，差分后序列时的时序图：



单位根检验：  



白噪声检验



检验结果显示，差分后序列时平稳非白噪声序列，需要对差分后的序列进行进一步拟合ARMA模型。

模型定阶



      自相关图显示出明显的下滑轨迹，这是典型的拖尾属性。偏自相关图除了1阶和4阶偏自相关系数显著大于2倍标准差。所以尝试拟合ARIMA(4,1,0)*(0,1,0)4

参数估计和模型检验



 

x2,x3，P>α，不通过显著性检验





模型的显著性检验：



检验结果显示，残差序列为白噪声序列，参数显著性检验显示两个参数均显著非0。

模型预测



 



> 到这里就结束了，如果对你有帮助你，欢迎点赞关注，你的点赞对我很重要

 

 

0 SARIMAX模型时间序列分析步骤

1. 用pandas处理时序数据

2. 检验时序数据的平稳性

3. 将时序数据平稳化

4. 确定order 的 p.d.q值

5. 确定season_order的四个值

6. 应用SARIMAX模型对时序数据进行预测

其实SARIMAX比ARIMA模型就多了个season_order参数的确定，但也是这里最费时间的一个步骤

1 将数据转化成为时序数据

先一股脑导入一下工具包

import pandas as pd
import datetime
import matplotlib.pyplot as plt
from pylab import mpl
mpl.rcParams['font.sans-serif']=['SimHei']
import seaborn as sns
import statsmodels.tsa.stattools as ts
import statsmodels.api as sm
from statsmodels.tsa.arima_model import ARIMA
from statsmodels.stats.diagnostic import unitroot_adf
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import itertools
import warnings
import numpy as np
from statsmodels.tsa.seasonal import seasonal_decompose

#读取数据
data = pd.read_csv('factor.csv')
data.index = pd.to_datetime(data['date'])
data.drop(['date'], axis=1, inplace=True)
data = data.result
data.head()

#数据大致情况展示
data.plot(figsize=(12,8))
plt.legend(bbox_to_anchor=(1.25, 0.5))
plt.title('result')
sns.despine()
plt.show()

 

2 序列平稳性检测

#数据平稳性检测 因为只有平稳数据才能做时间序列分析
def judge_stationarity(data_sanya_one):
    dftest = ts.adfuller(data_sanya_one)
    print(dftest)
    dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
    stationarity = 1
    for key, value in dftest[4].items():
        dfoutput['Critical Value (%s)'%key] = value 
        if dftest[0] > value:
                stationarity = 0
    print(dfoutput)
    print("是否平稳(1/0): %d" %(stationarity))
    return stationarity
stationarity = judge_stationarity(data)

3 序列平稳化 

#若不平稳进行一阶差分
if stationarity == 0:
    data_diff = data.diff()
    data_diff = data_diff.dropna()
    plt.figure()
    plt.plot(data_diff)
    plt.title('一阶差分')
    plt.show()

#再次进行平稳性检测
stationarity = judge_stationarity(data_diff)

 4 做一下季节性分解看看有没有季节性

#季节性分解
decomposition = seasonal_decompose(data,freq=28)
trend = decomposition.trend
seasonal = decomposition.seasonal
residual = decomposition.resid

plt.figure(figsize=[15, 7])
decomposition.plot()
print("test: p={}".format(ts.adfuller(seasonal)[1]))

#季节平稳性检测
stationarity = judge_stationarity(residual)

 5 对order参数p、q定阶，下面两种都可以，但图我还是看不来，下一种傻瓜式（稍微费时间）

#画ACF图和PACF图来确定p、q值
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

def draw_acf_pacf(ts,lags):
    f = plt.figure(facecolor='white')
    ax1 = f.add_subplot(211)
    plot_acf(ts,ax=ax1,lags=lags)  #lags 表示滞后的阶数，值为30，显示30阶的图像
    ax2 = f.add_subplot(212)
    plot_pacf(ts,ax=ax2,lags=lags)  
    plt.subplots_adjust(hspace=0.5)
    plt.show()
draw_acf_pacf(myts_diff,30)

#对模型p,q进行定阶
warnings.filterwarnings("ignore") # specify to ignore warning messages
from statsmodels.tsa.arima_model import ARIMA 

pmax = int(5)    #一般阶数不超过 length /10
qmax = int(5)
bic_matrix = []
for p in range(pmax +1):
    temp= []
    for q in range(qmax+1):
        try:
            temp.append(ARIMA(data, (p, 1, q)).fit().bic)
        except:
            temp.append(None)
        bic_matrix.append(temp)
 
bic_matrix = pd.DataFrame(bic_matrix)   #将其转换成Dataframe 数据结构
p,q = bic_matrix.stack().idxmin()   #先使用stack 展平， 然后使用 idxmin 找出最小值的位置
print(u'BIC 最小的p值 和 q 值：%s,%s' %(p,q))  #  BIC 最小的p值 和 q 值：0,1

6  通过网格搜索对seasonal_order进行定阶

#通过网格搜索对seasonal_order进行定阶,目前就是pdq=011，seasonal_order=2, 2, 1, 52效果比较好，RMSE=202.4582
def get_ARIMA_params(data, pdq, m=12):
    p = d = q = range(0, 3)
    seasonal_pdq = [(x[0], x[1], x[2], m) for x in list(itertools.product(p, d, q))]
    score_aic = 1000000.0
    warnings.filterwarnings("ignore") # specify to ignore warning messages
    for param_seasonal in seasonal_pdq:
        mod = sm.tsa.statespace.SARIMAX(data,
                                        order=pdq,
                                        seasonal_order=param_seasonal,
                                        enforce_stationarity=False,
                                        enforce_invertibility=False)
        results = mod.fit()
        print('x{}12 - AIC:{}'.format(param_seasonal, results.aic))
        if results.aic < score_aic:
            score_aic = results.aic
            params = param_seasonal, results.aic
    param_seasonal, results.aic = params
    print('x{}12 - AIC:{}'.format(param_seasonal, results.aic))
pdq = [0, 1, 1]
get_ARIMA_params(data, pdq, m=52)

这上面最关键的是这个m值怎么设定，我设置默认为12。m是季节周期，参考别人代码的时候，月度数据的m为12，那我这里以周为单位，季节周期应该是52吧。一年12个月，52个星期，这个逻辑应该没有问题。

 7  根据定阶参数进行模型拟合

mod = sm.tsa.statespace.SARIMAX(data,
                                order=(0, 1, 1),
                                seasonal_order=(2, 1, 2, 52),
                                enforce_stationarity=False,
                                enforce_invertibility=False)
results = mod.fit()
print(results.summary().tables[1])
results.plot_diagnostics(figsize=(15, 12))
plt.show()

这里模型拟合的时候用的数据的原始数据，而不是差分后的数据，因为order参数中已经设置了d为1，在拟合的时候会自动进行一阶差分，并在预测的时候对预测结果进行差分还原。

8  对预测值和真实值作图，并计算RMSE值作为评估参数

predict_ts = results.predict(tpy='levels')  #tpy='levels'直接预测值，没有的话预测的是差值
myts = data[predict_ts.index]  # 过滤没有预测的记录

predict_ts.plot(color='blue', label='Predict',figsize=(12,8))

myts.plot(color='red', label='Original',figsize=(12,8))

plt.legend(loc='best')
plt.title('RMSE: %.4f'% np.sqrt(sum((predict_ts-myts)**2)/myts.size))
plt.show()

 



9 向后对forecast值作图

steps = 20
start_time = myts.index[-1]
forecast_ts = results.forecast(steps)
 
fore = pd.DataFrame()
fore['date'] = pd.date_range(start=start_time ,periods=steps, freq='7D')
fore['result'] = pd.DataFrame(forecast_ts.values)
fore.index = pd.to_datetime(fore['date'])
 
predict_ts['2019/1/18':].plot(color='blue', label='Predict',figsize=(12,8))
myts['2019/1/18':].plot(color='red', label='Original',figsize=(12,8))
fore.result.plot(color='black', label='forecast',figsize=(12,8))

plt.legend(loc='best')
plt.show()

这里有个函数pd.date_range是专门用于产生时间序列索引的，start = 开始时间，end = 结束时间，periods=时间索引的个数，freq=‘7D’表示7天为一个时间索引间隔，也可以是'7W'七周，'M'一个月等等。

由于预测的数据没有时间索引，只有序号所以我要在这给他生成时间索引，并合并到dataframe，这样就可以和其他值一起在图像上展示了。



 

最后forecast的效果还是可以的嘛，保存forecast文件

fore.to_csv('forecast_20steps.csv')

 

1.季节性时间序列

包含：长期趋势Trend，季节趋势Seasonal，周期循环Circle，随机项Random

这里分解为相加模型X=T+S+C+R

 

在对时间序列进行分解之前，应该对序列进行检验：（下次写）

2.decompose()函数

将时间序列进行上述分解

3.R分解操作过程

3.1数据读入与可视化


>#以纽约市月出生数量（1946.1-1959.12）的数据集为例

> births <-scan("http://robjhyndman.com/tsdldata/data/nybirths.dat")

Read 168 items

> birthstimeseries <- ts(births, frequency=12, start=c(1946,1))

> plot(birthstimeseries)




出生数量

 

从图上可以看出，出生数量具有一定的季节性（夏峰冬谷）和周期性，同时趋势性明显；但是每个周期内的波动幅度变化较小，且不随时间趋势而变化，随便波动项随时间变化页不明显。

3.时间序列分解

分解为加法模型


>birthcomponents <- decompose(birthstimeseries)

> plot(birthcomponents)




分解图

4.剔除季节因素

可以对季节性等进行剔除，现剔除季节因素


>birthstimeseriesseasonallyadjusted<-birthstimeseries-birthcomponents$seasonal

>plot(birthstimeseriesseasonallyadjusted)




出生数量（剔除季节因素）

 


			
实验

			
名称
			
			

			
疏系数模型 和季节模型
			
		

			
实验

			
内容
			
			

			
1、简单季节模型
			
		

			
实验

			
目的
			
			

			
1、掌握疏系数模型

			
2、熟练建立季节模型
			
		
推荐阅读

使用Python完成时间序列分析基础
	
SPSS建立时间序列乘法季节模型实战案例
	
Python建立时间序列ARIMA模型实战案例
疏系数模型的定义



 

  



 

时序图


 

白噪声检验



一阶差分后的序列为白噪声，所以我们对二阶序列进行后续分析



 

检验结果显示，差分后序列时平稳非白噪声序列，需要对差分后的序列进行进一步模型拟合





ACF拖尾，PACF二阶截尾，建立ARIMA(2,2,0)

参数估计和模型检验



 



模型的显著性检验和参数的显著性检验均通过

建立模型口径：



 

 

 

 模型预测