实验
 
名称
 
 
疏系数模型 和季节模型
 
 
实验
 
内容
 
 
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。
 
模型预测
 
 
 
 
 
> 到这里就结束了，如果对你有帮助你，欢迎点赞关注，你的点赞对我很重要
 
 
 
 

摘要
 
时间序列模型就是将预测对象按照时间顺序排列起来,用这一组时间序列过去的变化规律,推断今后变化的可能性及变化趋势、变化规律。
 
时间序列模型也是一种回归模型,其一方面承认事物发展的延续性,运用过去的数据来推测事物的发展趋势;另一方面又考虑到偶然因素产生的随机性,为了消除随机波动的影响,利用历史数据,进行统计分析,并对数据进行适当的处理,进行趋势预测。
 
优点是简单易行,便于掌握，能够充分运用原时间序列的各项数据,计算速度快,对模型参数有动态确定的能力,精度较好,采用组合的时间序列或者把时间序列和其他模型组合效果更好。
缺点是不能反映事物的内在联系,不能分析两个因素的相关关系,只适用于短期预测。
 
确定性时间序列分析
 
确定性因素分解
 
长期趋势变动（T）：它是指时间序列朝着一 定的方向持续_上升或下降，或停留在某一水平上的倾向,它反映了客观事物的主要变化趋势。
季节变动(S)：受季节变动影响所形成的一种长度和幅度固定的周期性波动。
循环变动©：通常是指周期为一年以上，由非季节因素引起的涨落起伏波形相似的波动。
不规则变动®：通常它分为突然变动和随机变动。
 
常用的确定性时间序列模型有：
 
 
随机时间序列的两个重要性质
 
通常根平稳性和随机性来描述一个完全没有信息的时间序列，也称为白噪声。
 
严平稳
 
严平稳是一种条件比较苛刻的平稳性定义，它认为只有当序列所有的统计性质都不会随着时间的推移而发生变化时，该序列才能被认为平稳。
 
宽平稳
 
宽平稳是使用序列的特征统计量（均值方差等）来定义的一种平稳性。它认为序列的统计性质主要由它的低阶矩决定，所以只要保证序列低阶矩平稳(二阶)，就能保证序列的主要性质近似稳定。
 
平稳性检验
 
图检验方法(时序图检验、自相关图检验)
 
特点:操作简单，应用广泛，但结论带有一定主观
 性。
 
统计检验方法(单位根检验)
 
特点:根据特征根是否在单位圆内来判断平稳性,有多种类型，结果具有客观性。
 
移动平均法
 



 
 
 

 

 
 


 移动就是移动，平均就是平均，比如我们取4个数来移动平均。就是先用前四个数求平均，来代表第5个值；2： 5求平均代表第6个值；8：11求平均代表第12个值。（个人认为这没啥逻辑性。。。也可能是我理解错了，懂的同学还望指出，这是司守奎168页的内容）。
 这里我们选了4和5个数来实验，其评判标准为标准差，越小精度越好、
 
clc,clear
y=[533.8  574.6  606.9  649.8   705.1  772.0  816.4  892.7  963.9  1015.1  1102.7];
m=length(y);   
n=[4,5];             %n为移动平均的项数
for i=1:length(n)    %由于n的取值不同，下面使用了细胞数组
    for j=1:m-n(i)+1
        yhat{i}(j)=sum(y(j:j+n(i)-1))/n(i);
    end
    y12(i)=yhat{i}(end)  %提出第12月份的预测值
    s(i)=sqrt(mean((y(n(i)+1:end)-yhat{i}(1:end-1)).^2)) %求预测的标准误差
end
% y12, s       %分别显示两种方法的预测值和预测的标准误差
plot([1:11],y,[1:12],[y(1:4) yhat{1}])
legend('原始数据图','预测图')

y12 =
  993.6000  958.1600

s =
  150.5121  182.3851
 
993的标准差更小，所以99205更好，但作图看看。。。这是书上的实例可以自行去看
 
 
指数平滑法
 
一次移动平均实际上认为最近 N 期数据对未来值影响相同，都加权
 1/N ；而 N 期以前的数据对未来值没有影响，加权为 0。但是，二次及更高次移动平均数的权数却不是 1/N ，且次数越高，权数的结构越复杂，但永远保持对称的权数，即两端项权数小，中间项权数大，不符合一般系统的动态性。一般说来历史数据对未来值的影响是随时间间隔的增长而递减的。所以，更切合实际的方法应是对各期观测值依时间顺序进行加权平均作为预测值。指数平滑法可满足这一要求，而且具有简单的递推形式。指数平滑法根据平滑次数的不同，又分为一次指数平滑法、二次指数平滑法和三次指数平滑法等，分别介绍如下。
 
一次指数平滑法
 
和移动平均法差不多，只不过是用了前两次的平均值代表后一个数，然后加了个权重。
 还是看例子吧，数学公式写了也看不很明白。
 

 例题、某市 1976～1987 年某种电器销售额如表 4 所示。试预测 1988 年该电器销售额。
 
解 采用指数平滑法，并分别取α = 0.2,0.5和0.8进行计算，初始值
 
 

  
 
 
 
 


 

  表 4 某种电器销售额及指数平滑预测值计算表 （单位：万元） 
 
 

clc,clear
yt=[50 52 47 51 49 48 51 40 48 52 51 59]';   %实际销售额数据以列向量的方式存放在纯文本文件中
n=length(yt);
alpha=[0.2 0.5 0.8];
m=length(alpha);
yhat(1,[1:m])=(yt(1)+yt(2))/2;
for i=2:n
    yhat(i,:)=alpha*yt(i-1)+(1-alpha).*yhat(i-1,:);
end
yhat
err=sqrt(mean((repmat(yt,1,m)-yhat).^2)) 
xlswrite('dianqi.xls',yhat) %把预测数据写到Excel文件，准备在word表格中使用
yhat1988=alpha*yt(n)+(1-alpha).*yhat(n,:)

yhat =
   51.0000   51.0000   51.0000
   50.8000   50.5000   50.2000
   51.0400   51.2500   51.6400
   50.2320   49.1250   47.9280
   50.3856   50.0625   50.3856
   50.1085   49.5313   49.2771
   49.6868   48.7656   48.2554
   49.9494   49.8828   50.4511
   47.9595   44.9414   42.0902
   47.9676   46.4707   46.8180
   48.7741   49.2354   50.9636
   49.2193   50.1177   50.9927

err =
    4.5029    4.5908    4.8426

yhat1988 =
   51.1754   54.5588   57.3985
 
这里三个权值对应的标准差相差并不大，但肯定要取最小的，所以是4.5也就是权值为0.2的那个，计算出值为51.1754，这个数据本身就不是很好，感觉好没有太大的变化规律或者趋势，预测的肯定也不怎么样。
 
 这个模型用于第一个例题试试
 
 相对来说好多了，最终预测值为1082
 
二次指数平滑法
 
仍以例 3 我国 1965～1985 年的发电总量资料为例，试用二次指数平滑法预
 测 1986 年和 1987 年的发电总量。
 


 表 6 我国发电总量及一、二次指数平滑值计算表
 

 
 
 

 
clc,clear
yt=[676,825,774,716,940,1159,1384,1524,1668,1688,1958,2031,2234,2566,2820,3006,3093,3277,3514,3770,4107];
n=length(yt), alpha=0.3; st1(1)=yt(1); st2(1)=yt(1);

for i=2:n
    st1(i)=alpha*yt(i)+(1-alpha)*st1(i-1);%一次平滑值
    st2(i)=alpha*st1(i)+(1-alpha)*st2(i-1);%二次平滑值，和一次一个样，相当于把一次平滑值的结果作为原式数据进行计算
end
xlswrite('fadian.xls',[st1',st2']) %把数据写入表单Sheet1中的前两列
at=2*st1-st2;
bt=alpha/(1-alpha)*(st1-st2);
yhat=at+bt;  %最后的一个分量为1986年的预测值
xlswrite('fadian.xls',yhat','Sheet1','C2') %把预测值写入第3列
str=['C',int2str(n+2)] %准备写1987年预测值位置的字符串
xlswrite('fadian.xls',at(n)+2*bt(n),'Sheet1',str)%把1987年预测值写到相应位置
plot([1965:1985],yt,[1966:1986], yhat)
legend('原始数据图','预测图')
title("对比图")
 
得到86 87两年的数据为
 4.2239e+03
 4.4342e+03
 总体来看效果可以
 
 还是用于第一题的数据
 
 简直翻天覆地啊。这个很可以。
 
三次指数平滑法
 
某省 1978～1988 年全民所有制单位固定资产投资总额如表 7 所示，试预测
 1989 年和 1990 年固定资产投资总额
 



 表 7 某省全民所有制单位固定资产投资总额及一、二、三次指数平滑值计算表（单位：亿元)
 
  

 

 

 
 
clc,clear
yt=load('touzi.txt'); %原始投资总额数据以列向量的方式存放在纯文本文件中
n=length(yt); alpha=0.25; st0=mean(yt(1:3));
st1(1)=alpha*yt(1)+(1-alpha)*st0;
st2(1)=alpha*st1(1)+(1-alpha)*st0;
st3(1)=alpha*st2(1)+(1-alpha)*st0;
for i=2:n
    st1(i)=alpha*yt(i)+(1-alpha)*st1(i-1);
    st2(i)=alpha*st1(i)+(1-alpha)*st2(i-1);
    st3(i)=alpha*st2(i)+(1-alpha)*st3(i-1);
end
xlswrite('touzi.xls',[st1',st2',st3'])  %把数据写在前三列
at=3*st1-3*st2+st3;
bt=0.5*alpha/(1-alpha)^2*((6-5*alpha)*st1-2*(5-4*alpha)*st2+(4-3*alpha)*st3);
ct=0.5*alpha^2/(1-alpha)^2*(st1-2*st2+st3);
yhat=at+bt+ct;
xlswrite('touzi.xls',yhat','Sheet1','D2')  %把数据写在第4列第2行开始的位置
plot(1:n,yt,'D',2:n,yhat(1:end-1),'*')
legend('实际值','预测值','Location','northwest')  %图注显示在左上角
xishu=[ct(end),bt(end),at(end)]; %二次预测多项式的系数向量
yhat1990=polyval(xishu,2)  %求预测多项式m=2时的值
plot([1:11],yt,[1:11], yhat)
legend('原始数据图','预测图')
title("对比图")
 

 用于第一题，看起来也不错。
 
 后面的有点难，自己没怎么看懂，先不写了吧 😂

 
目录
 
0、前言
一、数据准备&探索
1、平稳性
1.1 平稳性检验
1.2 数据处理（平滑、变换、差分、分解）
1.2.1 对数变换
1.2.2 平滑法（移动平均&指数平均）
1.2.3 差分法
1.2.4 分解
   
   
2、非白噪声检验
  
二、模型（ARIMA）
1、 自回归移动平均模型（ARMA）
2、 模型参数确定
2.1 PACF和ACF 图 找到
2.2 信息准则确定
2.3 热力图确定
   
3、模型预测
3、 其他模型
  
参考文献

 
0、前言
 
什么是时间序列
 时间序列简单的说就是各时间点上形成的数值序列，时间序列分析就是通过观察历史数据预测未来的值。在这里需要强调一点的是，时间序列分析并不是关于时间的回归，它主要是研究自身的变化规律的（这里不考虑含外生变量的时间序列）。
 时间序列的分析的步骤是先对数据进行平稳性和非白噪声检验（如不满足需对数据进行平滑或差分等预处理），然后才是模型调参跟预测，因此本文分为2大部分介绍，答题思路见脑图。
 
 
一、数据准备&探索
 
1、平稳性
 
序列平稳性是进行时间序列分析的前提条件为什么要满足平稳性的要求呢？在大数定理和中心定理中要求样本同分布（这里同分布等价于时间序列中的平稳性），而我们的建模过程中有很多都是建立在大数定理和中心极限定理的前提条件下的，如果它不满足，得到的许多结论都是不可靠的。以虚假回归为例，当响应变量和输入变量都平稳时，我们用t统计量检验标准化系数的显著性。而当响应变量和输入变量不平稳时，其标准化系数不在满足t分布，这时再用t检验来进行显著性分析，导致拒绝原假设的概率增加，即容易犯第一类错误，从而得出错误的结论。
 判断一个序列是不是平稳序列有三个评判标准：
 
均值 ，是与时间t 无关的常数。
方差 ，是与时间t 无关的常数。这个特性叫做方差齐性。
协方差 ，只与时期间隔k有关，与时间t 无关的常数。
 
如下所示分别为稳定的数据和均值不稳定、方差不稳定、自协方差随时间变化
 
1.1 平稳性检验
 
时间序列平稳性检验主要是观察法和单位根检验法，观察法
 观察法，通俗的说就是通过观察序列的趋势图与相关图是否随着时间的变化呈现出某种规律。所谓的规律就是时间序列经常提到的周期性因素，现实中遇到得比较多的是线性周期成分，这类周期成分可以采用差分或者移动平均来解决，而对于非线性周期成分的处理相对比较复杂，需要采用某些分解的方法。
 
一般都是使用单位根检验法 Dickey-Fuller Test 进行判断，即在一定置信水平下，假设时序数据Null hypothesis: 非稳定、序列具有单位根。因此对于一个平稳的时序数据，就需要在给定的置信水平上显著，拒绝原假设。如下通过statsmodels包中的adfuller来进行检测：
 数据下载
 
import os
import sys
import pandas as pd
import pandas_datareader.data as web
import numpy as np
import statsmodels.formula.api as smf
import statsmodels.tsa.api as smt
import statsmodels.api as sm
import scipy.stats as scs
from arch import arch_model
import matplotlib.pyplot as plt
import matplotlib as mpl
from statsmodels.stats.diagnostic import acorr_ljungbox #白噪声检验
#######读取数据等
data_ts = pd.read_csv("test0923.csv")
data_ts.columns = ['week_no','send_date','day_week','send_date_pre','export_type','depart_thr_letter_code',
       'send_weight','real_weight_qty','meterage_weight_qty','send_weight_pre','real_weight_qty_pre','meterage_weight_qty_pre']

data_ts.columns
data_ts.head()
data_ts['send_date'] = pd.to_datetime(data_ts['send_date'])
data_ts.dtypes
data_ts.index=data_ts['send_date']
df_data_ts=data_ts[data_ts['export_type']=='HKG']
#ads=df_data_ts['send_weight'][-100:]
ads=df_data_ts['send_weight'][-65:]

###########ADF检验

adftest = adfuller(ads, autolag='AIC')
print (adftest)
adfuller(np.diff(ads), autolag='AIC')#检查不通过用一阶差分后的数据看看，即np.diff(ads)
 
结果如下：
 
 如何确定该序列能否平稳呢？主要看：
 （1）1%、%5、%10不同程度拒绝原假设的统计值和ADF Test result的比较，ADF Test result同时小于1%、5%、10%即说明非常好地拒绝该假设，本数据中，adf结果为-10.6， 小于三个level的统计值，可以拒绝原假设，即数据是平稳的。
 （2）P-value是否非常接近0.本数据中，P-value 为 4.8e-19，接近0。
 
1.2 数据处理（平滑、变换、差分、分解）
 
如果，序列不平稳，或者数据预测不是很好，就要考虑进行数据处理。处理方法主要有以下四种：
 
1.2.1 对数变换
 
对数变换主要是为了减小数据的振动幅度，使其线性规律更加明显。对数变换相当于增加了一个惩罚机制，数据越大其惩罚越大，数据越小惩罚越小。这里强调一下，变换的序列需要满足大于0，小于0的数据不存在对数变换。
 
ads_log = np.log(ads)
adfuller(ads_log, autolag='AIC')
 
1.2.2 平滑法（移动平均&指数平均）
 
移动平均即利用一定时间间隔内的平均值作为某一期的估计值，也有升级版的加权平均。参考代码，如下：
 
############ 移动平均图----------
def draw_trend(timeSeries, size):
    f = plt.figure(facecolor='white')
    # 对size个数据进行移动平均,min_periods=1 则前面size数不为空
    rol_mean = timeSeries.rolling(window=size,min_periods=1).mean()
    # 对size个数据进行加权移动平均
    # 加权移动平均是越近是越远的2倍
    #rol_weighted_mean = pd.ewma(timeSeries, span=size)
    rol_weighted_mean=timeSeries.ewm(span=size,min_periods=1).mean()

    timeSeries.plot(color='blue', label='Original')
    rol_mean.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()
    return rol_mean
    
ads_rol_mean=draw_trend(ads,3)    
 
用平均法的数据带入模型后预测结果需要还原
 
######移动平均还原
ads_rec=ads_rol_mean*3 - ads.rolling(window=2,min_periods=1).sum().shift(1)
##前size个数值需专门计算 

##or 另一种还原方法（直接计算肯定更不容易出错）
ads_rec=ads_rol_mean
for i in range(len(ads_rol_mean)) :
    print(i)
    if i<=1 :
        ads_rec[0]=ads_rol_mean[0]
        ads_rec[1]=ads_rol_mean[1]*2-ads_rol_mean[0]
    else :
        ads_rec[i]=ads_rol_mean[i]*3-ads_rec[i-1]-ads_rec[i-2]
 
指数平均则是用变权的方法来计算均值。具体怎么做的呢？一个简单的数学式为：
 $$y_t^s=\mu \ast y_t+(1-\mu )\ast y_{t-1}^s,$$
 式子中的 $\mu$表示平滑因子，它定义我们“遗忘”当前真实观测值的速度有多快。 $\mu$越小，表示当前真实观测值的影响力越小，而前一个模型预测值的影响力越大，最终得到的时间序列将会越平滑。
 那么指数体现在哪呢？指数就隐藏在递归函数之中，我们上面的函数，每次都要用( 1 − $\mu$)乘以模型的上一个预测值。
 代码可参考如下（自己项目未用的到，代码来自这篇博客）：
 
def exponential_smoothing(series, alpha):
    """
        series - dataset with timestamps
        alpha - float [0.0, 1.0], smoothing parameter
    """
    result = [series[0]] # first value is same as series
    for n in range(1, len(series)):
        result.append(alpha * series[n] + (1 - alpha) * result[n-1])
    return result
    
def plotExponentialSmoothing(series, alphas):
    with plt.style.context('seaborn-white'):    
        plt.figure(figsize=(15, 7))
        for alpha in alphas:
            plt.plot(exponential_smoothing(series, alpha), label="Alpha {}".format(alpha))
        plt.plot(series.values, "c", label = "Actual")
        plt.legend(loc="best")
        plt.axis('tight')
        plt.title("Exponential Smoothing")
        plt.grid(True);
        
plotExponentialSmoothing(ads.Ads, [0.3, 0.05])
plotExponentialSmoothing(currency.GEMS_GEMS_SPENT, [0.3, 0.05])
 
上述是但指数平滑法，还有双指数平滑法等，可参考上面的博客。
 
1.2.3 差分法
 
时间序列最常用来剔除周期性因素的方法当属差分了，它主要是对等周期间隔的数据进行线性求减。
 参考代码如下：
 
# 进行差分
ads=pd.Series(ads) ###也可不用pd.Series，用np.diff(ads)
ads.head()
ads_diff=ads.diff(1)
ads_diff.head()
ads_diff=ads_diff.dropna()
 
############## 差分还原
diff_shift=ads.shift(1)
ads_recov=ads_diff.add(diff_shift)
ads_recov[0]=ads[0]
# or第二种方法
ads_recov=ads_diff  ###看0 is nan?
for i in range(len(ads_diff)) :
    print(i)
    if i<=0 :
        ads_recov[0]=ads[0]
    else :
        ads_recov[i]=ads_recov[i]+ads_recov[i-1]
############
 
1.2.4 分解
 
分解就是将时序数据分离成不同的成分。statsmodels使用的X-11分解过程，它主要将时序数据分离成长期趋势、季节趋势和随机成分。与其它统计软件一样，statsmodels也支持两类分解模型，加法模型和乘法模型，这里我只实现加法，乘法只需将model的参数设置为"multiplicative"即可。本人未用该方法，代码参考于这篇文章
 
from statsmodels.tsa.seasonal import seasonal_decompose
decomposition = seasonal_decompose(ads, model="additive")

trend = decomposition.trend
seasonal = decomposition.seasonal
residual = decomposition.resid
 
2、非白噪声检验
 
非白噪声检验也可以用初步观察的方法，可用如下代码观察残差图等：
 
def tsplot(y, lags=None, figsize=(10, 8), style='bmh'):
    if not isinstance(y, pd.Series):
        y = pd.Series(y)
    with plt.style.context(style):    
        fig = plt.figure(figsize=figsize)
        #mpl.rcParams['font.family'] = 'Ubuntu Mono'
        layout = (3, 2)
        ts_ax = plt.subplot2grid(layout, (0, 0), colspan=2)
        acf_ax = plt.subplot2grid(layout, (1, 0))
        pacf_ax = plt.subplot2grid(layout, (1, 1))
        qq_ax = plt.subplot2grid(layout, (2, 0))
        pp_ax = plt.subplot2grid(layout, (2, 1))
        
        y.plot(ax=ts_ax)
        ts_ax.set_title('Time Series Analysis Plots')
        smt.graphics.plot_acf(y, lags=lags, ax=acf_ax, alpha=0.5)
        smt.graphics.plot_pacf(y, lags=lags, ax=pacf_ax, alpha=0.5)
        sm.qqplot(y, line='s', ax=qq_ax)
        qq_ax.set_title('QQ Plot')        
        scs.probplot(y, sparams=(y.mean(), y.std()), plot=pp_ax)

        plt.tight_layout()
    return 
 
非白噪声检验代码如下：
 
###########纯随机性检验（白噪声检验）
p_value = acorr_ljungbox(ads, lags=1) 
print (p_value)
acorr_ljungbox(np.diff(ads), lags=1)
 
结果如图：
 
 统计量的P值小于显著性水平0.05，则可以以95%的置信水平拒绝原假设，认为序列为非白噪声序列（否则，接受原假设，认为序列为纯随机序列。）
 
由于P值为3.9e-5远大于0.05所以拒绝原假设，认为时间序列是非白噪声的，即非是随机产生的序列，具有时间上的相关性。
 从严格意义上说，只有通过平稳性检验和非白噪声检验才能用时间序列模型。不过很多博客都没有非白噪声检验这一步，可能这跟数据分析挖掘最终只看效果有关吧。
 
二、模型（ARIMA）
 
1、 自回归移动平均模型（ARMA）
 
ARIMA比ARMA模型多了I这个差分，其他一样。自回归移动平均模型(ARMA(p，q))是时间序列中最为重要的模型之一，它主要由两部分组成： AR代表p阶自回归过程，MA代表q阶移动平均过程，其公式特征如下：
 
 
AR§：autoregression model(自回归模型)，即时间序列对自身的回归。基本假设是当前的序列值取决于它之前的值，并且存在一定的滞后。 模型最大的滞后值称为 p ，要确定初始的 p 值，我们需要在 PACF 图中找到最大的滞后点。
MA(q)：moving average model(移动平均值模型)，这里不讨论它的细节，它基于当前的误差依赖于先前的误差，有一定的滞后性（滞后值记为q） 这一假设，对时间序列进行建模。和上面一样，可以在 ACF 图找到滞后值 q qq 的初始值；
 
详细例子和说明可参考该文章
 
2、 模型参数确定
 
根据ARIMA模型理论很容易知道我们需要确定 p(自回归滞后)、d（差分）、q（误差滞后）三个参数。主要有如下三种方法
 
2.1 PACF和ACF 图 找到
 
p值可以通过 PACF 图找到，q值可以通过 ACF 图找到。d可以根据差分后的平稳性检验得到。参考代码如下：
 
def tsplot(y, lags=None, figsize=(10, 8), style='bmh'):
    if not isinstance(y, pd.Series):
        y = pd.Series(y)
    with plt.style.context(style):    
        fig = plt.figure(figsize=figsize)
        #mpl.rcParams['font.family'] = 'Ubuntu Mono'
        layout = (3, 2)
        ts_ax = plt.subplot2grid(layout, (0, 0), colspan=2)
        acf_ax = plt.subplot2grid(layout, (1, 0))
        pacf_ax = plt.subplot2grid(layout, (1, 1))
        qq_ax = plt.subplot2grid(layout, (2, 0))
        pp_ax = plt.subplot2grid(layout, (2, 1))
        
        y.plot(ax=ts_ax)
        ts_ax.set_title('Time Series Analysis Plots')
        smt.graphics.plot_acf(y, lags=lags, ax=acf_ax, alpha=0.5)
        smt.graphics.plot_pacf(y, lags=lags, ax=pacf_ax, alpha=0.5)
        sm.qqplot(y, line='s', ax=qq_ax)
        qq_ax.set_title('QQ Plot')        
        scs.probplot(y, sparams=(y.mean(), y.std()), plot=pp_ax)

        plt.tight_layout()
    return 

tsplot(df_data_ts['send_weight'], lags=30)
 
2.2 信息准则确定
 
就是遍历所有参数值，看信息准则最小的（AIC/BIC/HQIC）,一般选取AIC即可。
 信息准则的区别可参考文章
 
best_aic = np.inf 
best_order = None
best_mdl = None
#一般阶数不超过整体数据的十分之一，因此这里采用穷举法。分别从0~10取p、q。
pq_rng = range(10) # [0,1,2,3,4]
d_rng = range(3) # [0,1] ARMA
for i in pq_rng:
    for d in d_rng:
        for j in pq_rng:
            try:
                tmp_mdl = smt.ARIMA(ads, order=(i,d,j)).fit(method='mle', trend='nc')
                tmp_aic = tmp_mdl.aic
                if tmp_aic < best_aic:
                    best_aic = tmp_aic
                    best_order = (i, d, j)
                    best_mdl = tmp_mdl
            except: continue


print('aic: {:6.5f} | order: {}'.format(best_aic, best_order))
_ = tsplot(best_mdl.resid, lags=8)
print(best_mdl.summary())

## 德宾-沃森（Durbin-Watson）检验
print(sm.stats.durbin_watson(best_mdl.resid.values))
 
【注意】 代码中专门加了德宾-沃森（Durbin-Watson）检验。对于一阶差分，当D-W检验值接近于2时，不存在自相关性，说明模型较好。
 
2.3 热力图确定
 
其实热力图定阶的方式和信息准则定阶的方式类似，只是用热力图的方式呈现了。
 

#设置遍历循环的初始条件，以热力图的形式展示，跟AIC定阶作用一样
    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(ads, 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、模型预测
 
模型预测，有2个函数predict()、forecast()，2个函数的各自参数不做介绍，自己查看帮助文档。区别除了模型参与以外，有以下2个重要区别：
 
predict()可以预测样本内和样本外，forecast只能预测样本外，而且第一个数必须是样本最后一个。predict()函数进行预测，其可以传入start、end参数代表预测数据的开始和结束坐标，坐标值不仅可以是int、str，还可以是 datetime时间类型，如果start、end介于原有数据区域内，即为对原数据的预测拟合，如果end超过了原有数据长度，即代表对未来数据进行预测。forecast()函数接收一个steps参数，代表对未来多少个数据进行预测，也可以传入datetime时间，代表从样本数据结束一直预测到该时间。
如果模型参数d=0不需要差分，则2个函数预测结果一样。若d>0需要差分时,则predict()返回的是差分值，而forecast()则返回的是原始值，从这一点上讲建议使用forecast()。
 
参考代码分别如下：
 
best_aic = np.inf 
best_order = None
best_mdl = None

pq_rng = range(10) # [0,1,2,3,4]
d_rng = range(2) # [0,1] ARMA
for i in pq_rng:
    for d in d_rng:
        for j in pq_rng:
            try:
                tmp_mdl = smt.ARIMA(ads, order=(i,d,j)).fit(method='mle', trend='nc')
                tmp_aic = tmp_mdl.aic
                if tmp_aic < best_aic:
                    best_aic = tmp_aic
                    best_order = (i, d, j)
                    best_mdl = tmp_mdl
            except: continue


print('aic: {:6.5f} | order: {}'.format(best_aic, best_order))
_ = tsplot(best_mdl.resid, lags=8)
print(best_mdl.summary())
## 德宾-沃森（Durbin-Watson）检验
 print(sm.stats.durbin_watson(best_mdl.resid.values))
 
 ###最好用forecast，不用predict，predict在需要查分时数据结果需要还原,2种方法结果不同
n_steps = 14
f, err95, ci95 = best_mdl.forecast(steps=n_steps) # 95% CI
_, err99, ci99 = best_mdl.forecast(steps=n_steps, alpha=0.01) # 99% CI
##预测的第一个数是样本最后的日期
idx = pd.date_range(ads.index[-1], periods=n_steps, freq='D')
fc_95 = pd.DataFrame(np.column_stack([f, ci95]), 
                     index=idx, columns=['forecast', 'lower_ci_95', 'upper_ci_95'])
fc_99 = pd.DataFrame(np.column_stack([ci99]), 
                     index=idx, columns=['lower_ci_99', 'upper_ci_99'])
fc_all = fc_95.combine_first(fc_99)
fc_all.head()

############predict
pred = best_mdl.predict(1,66)#不加开始和结束数据则为原始时间索引
########差分还原------
pred= best_mdl.predict()
pred= pred.add(ads.shift(1))
#这种方法简单但不能还原第一个数
 
3、 其他模型
 
还有增加季节因素的SARIMA模型，以及ARCH、GARCH模型，可参考文章。
 
参考文献
 
1、http://www.blackarbs.com/blog/time-series-analysis-in-python-linear-models-to-garch/11/1/2016
 2、https://www.cnblogs.com/foley/p/5582358.html
 3、https://blog.csdn.net/theVicTory/article/details/104941483
 4、https://blog.csdn.net/Earl211/article/details/50957029
 5、https://blog.csdn.net/jh1137921986/article/details/90257764
 6、https://blog.csdn.net/foneone/article/details/90141213?utm_medium=distribute.pc_aggpage_search_result.none-task-blog-2allfirst_rank_v2~rank_v25-2-90141213.nonecase&utm_term=arima%E4%BA%8C%E9%98%B6%E5%B7%AE%E5%88%86%E7%9A%84%E9%A2%84%E6%B5%8B%E5%80%BC%E8%BF%98%E5%8E%9F#%EF%BC%883%EF%BC%89ADF%E6%A3%80%E9%AA%8C

 
 
实验数据：链接: https://pan.baidu.com/s/1SFy2Zc6A6KAT5rVUVUMWEA 提取码: 8191
 
 
本部分过于复杂，以例题形式讲解
 
 
例题1：某产品销售数据（2014年第1季度-2018年第4季度）如下表，根据此表预测之后到2020年第四季度的销量。
 
（数据文件为——“某产品的销售数据（2014-2018年季度数据）.sav”）
 
 
spss操作过程如下：
 
1、首先根据数据做出数据的时间序列图
 
但在此之间要构建虚变量，便于观察（即将带中文的或复杂表达变成一个或几个表达）
 
 
之后做出时间序列图（即折线图）