【MATLAB第82期】基于MATLAB的季节性差分自回归滑动平均模型SARIMA时间序列预测模型含预测未来

一、模型介绍

1、模型简介

季节性差分自回归移动平均模型（Seasonal Autoregressive Integrated Moving Average Model, SARIMA），又称为周期性差分自回归移动平均模型，是时间序列预测常用的分析方法之一，常应用于包含趋势和季节性的单变量数据的预测。SARIMA对于时间序列数据的季节性变动等周期性属性具有较高的敏感性。

2、模型参数

SARIMA结构参数有七个：(p,d,q) (P,D,Q,s)

1、季节性与非季节性差分数
d：代表时序数据需要进行几阶差分化，才是稳定的，也叫Integrated项。使用SARIMA模型要求数据平稳。SARIMA的差分项有两个，分别是季节性差分D与非季节性差分d。通常季节性差分经过一次即可，非季节性差分通常在0~3之间。确定非季节性差分数d从0至3循环，平稳后停止，当检验模型参数时d=1时数据已经平稳。

D: 周期性差分阶数。季节性差分通常在0~3之间。确定季节性差分数D从0至3循环，平稳后停止，当检验模型参数时D=1时数据已经平稳。

2、确定SARIMA模型阶数
这个步骤中需要确定的阶数有四个：AR阶数p，MA阶数q，SAR阶数P，SMA阶数Q。用基于AICBIC准则的方法定阶。
p：代表预测模型中采用的时序数据本身的滞后数(lags) ,也叫做AR/Auto-Regressive项。p通常在1~3之间，通过循环可得p=2时，AICBIC值最小。
q：代表预测模型中采用的预测误差的滞后数(lags)，也叫做MA/Moving Average项。q通常在1~3之间，通过循环可得q=1时，AICBIC值最小。
P: 周期性自回归阶数。P通常在1~2之间，通过循环可得P=2时，AICBIC值最小。
Q: 周期性移动平均阶数。Q通常在1~2之间，通过循环可得Q=1时，AICBIC值最小。

3、确定周期时间间隔
s: 周期时间间隔。过对现有数据集的观察，案例数据在3年尺度上存在一定的周期性，可以认为既定的年份的数据与其前一年的数据（也即前36个月的数据）之间存在较强的相关性，故选取36为周期时间间隔。

4、残差检测
为了确保确定的阶数合适，还需要进行残差检验。残差即原始信号减掉模型拟合出的信号后的残余信号。如果残差是随机正态分布的、不自相关的，这说明残差是一段白噪声信号，也就说明有用的信号已经都被提取到模型中了
在这里插入图片描述
上图为残差检验的结果图。Standardized Residuals是查看残差是否接近正态分布，理想的残差要接近正态分布；ACF和PACF检验残差的自相关和偏自相关，理想的结果应该在图中不存在超出蓝线的点；最后一张QQ图是检验残差是否接近正太分布的，理想的结果中蓝点应该靠近红线。
除了上述图像检验方法，还可以通过Durbin-Watson对相关性进行检验：
Durbin-Watson 统计是计量经济学分析中最常用的自相关度量，该值接近2，则可以认为序列不存在一阶相关性。
运算结果为1.96，这个值越接近2越说明残差不存在一阶相关性。
上述检验可以证明，残差接近正态分布，且相互独立，可以认为SARIMA建模符合要求。

二、预测效果

在这里插入图片描述

表5 SARIMA训练集和测试集预测结果评价指标
样本 RMSE MSE MAE R2
训练集 2.101 4.4143 0.88828 0.93821
测试集 1.6282 2.6511 1.1039 0.97196

从评价指标数值表的角度看，SARIMA在销量预测上有较高的精度，测试所得的误差值较小，能深入挖掘长时间序列数据的深层规律。从对比LSTM图像的角度看，SARIMA对周期性波动性的挖掘较为深入，但对整体态势的预测上表现较LSTM差，对个别点的预测不够精确，缺乏对离散数值的关注。

三、部分代码展示

close all
clear all
%% 1.加载数据
xall= importdata('经营数据2.xlsx');%导入数据
time=xall.textdata;%时间数据
xnum = datenum(time(2:end,1)); % 将日期转为数值
data= xall.data(:,2);% 时间序列数据 

data1=data;
addpath('funs');
S = 36; %季节性序列变化周期
step = 12;
% 通常P和Q不大于3
%% 2.确定季节性与非季节性差分数，D取默认值1，d从0至3循环，平稳后停止
for d = 0:3
          
    dY = filter(D,data)%对原数据进行差分运算  
    if(getStatAdfKpss(dY)) %数据平稳
        disp(['非季节性差分数为',num2str(d),'，季节性差分数为1']);
        break;
    end
end
%% 3.确定阶数ARlags,MALags,SARLags,SMALags
max_ar = 3;    %ARlags上限
max_ma = 3;    %MALags上限
max_sar = 2;   %SARLags上限
max_sma = 2;   %SMALags上限