时间序列分析2|ARIMA模型|SARIMA模型

ARMA模型的定阶

自相关和偏自相关系数法

通过观察样本的自相关系数(ACF)和偏自相关系数(PACF)，进行大体的判断
![[Pasted image 20240820155402.png]]

模型定阶的经验方法

截尾：

最初的d阶样本(偏)自相关系数明显在2倍标准差范围外
95%的(偏)自相关系数都落在2倍标准差的范围以内
非零自相关系数衰减为小值波动的过程非常突然
这时，通常视为d阶截尾

拖尾：

有超过5%的相关系数落在2倍标准差范围之外
显著非零的相关系数衰减为小值波动的过程比较缓慢或者非常连续
通常视为相关系数拖尾

缺点

由于样本随机性，§ACF是截尾还是拖尾有时仍然会难以准确判断
ARMA(p.,q)模型阶数p和q一般难以准确判断，主观性强

信息准则函数定阶法

基于相关系数图形的定阶法具有很强的主观性，是一种较为粗略的方法
信息准则定阶法则可以帮助我们在一定标准下从待定模型中选择相对最优的模型

AIC信息准则
AIC准测函数(最小化)
$AIC(K)=-2\ln(L)+2K$
L是似然函数，K是模型参数个数
简化为
$AIC(K)=n\ln(\hat{\sigma}^{2})+2K$
$\hat{siamg}^{2}$ 是模型误差方差的ML估计(每一个样本的对应的拟合的误差的平方和再除以自由度)，n是样本个数
说明：

既考虑模型预测准确性(第一项)，也考虑模型本身的复杂程度(第二项)
在一定程度上避免过拟合
但是，当样本数量n较大时，更倾向于选择更高阶模型。

BIC信息准则
BIC准测函数(最小化)
$BIC(K)=-2\ln(L)+\ln(n)K$
L是似然函数，K是模型参数个数
说明：

既考虑模型预测准确性(第一项)，也考虑模型本身的复杂程度(第二项)
与AIC相比，BIC进一步加大了对模型复杂性的惩罚
更加倾向于选择参数少的简单模型

ARIMA模型的概念

非平稳时间序列

在实际应用中，经常会遇见不满足平稳性的时间序列，尤其在经济领域和商业领域中，多数时间序列都是平稳的

时间序列模型

平稳时间序列

特征：均值、方差，自协方差都是稳定的，与起止时间无关。
建模：ARMA模型
非平稳时间序列
特征：均值非平稳，方差和自协方差非平稳
处理方法：差分运算，平稳化变换等
建模：ARIMA模型，SARIMA模型
均值非平稳
均值非平稳性
特征：均值时变
常用模型
确定性趋势模型：把时间t作为自变量，序列观测值作为因变量，建立
序列值随时间变化的回归模型
随机趋势模型：差分运算、ARIMA模型，SARIMA模型

差值运算

一阶差分
$\nabla X_{t}=X_{t}-X_{t-1}=(1-B)X_{t}$
d阶差分
$\nabla^{d}X_{t}=\nabla^{d-1}X_{t}-\nabla^{d-1}X_{t-1}=(1-B)^{d}X_{t}$
k步差分
$\nabla_{k}X_{t}=X_{t}-X_{t-k}=(1-B^{k})X_{t}$

差分方法是一种非常简便、有效的确定性信息提取方法
差分运算的实质是使用自回归的方式提取确定性信息
$\nabla^{d}X_{t}=(1-B)^{d}X_{t}=\sum_{i=0}^{d}(-1)^{i}\begin{pmatrix} d \\ i \end{pmatrix}X_{t-i}$
实际中，很多非平稳时间序列都可以通过适当的差分运算，转化为平稳序列
蕴含着显著线性趋势的序列，一阶差分就可以实现趋势平稳；
蕴含着曲线趋势的序列，通常低阶（二阶或三阶）差分就可以提取出曲
线趋势的影响，实现趋势平稳；
对于蕴含着固定周期的序列，通过步长为周期长度的差分，提取出周
期的影响，实现趋势平稳。

ARIMA模型的定义

如果时间序列 $\left\{ X_{t}\right\}$ 的d阶差分 $Y_{t}=(1-B)^{d}X_{t}$ 是一个平稳的ARMA(p,q)序列，其中 $d\ge 1$ 是整数，则称 $X_{t}$ 为具有阶p，d和q的ARIMA模型
记为
$X_{t}\sim ARIMA(p,d,q)$
ARIMA模型的结构
$\left\{\begin{matrix} \Phi(B)\nabla^{d}X_{t}=\Theta(B)\varepsilon_{t} \\ \varepsilon_{t}\sim N(0,\sigma^{2}) \\ E(X_{s}\varepsilon_{t})=0, \forall s<t \end{matrix}\right.$
其中
$\nabla^{d}=(1-B)^{d}$
ARMA(p,q)的自回归系数多项式
$\Phi(B)=1-\phi_{1}B-\phi_{2}B^{2}-\dots-\phi_{p}B^{p}$
ARMA(p,q)的移动平均系数多项式
$\Theta(B)=1-\theta_{1}B-\theta_{2}B^{2}-\dots-\theta_{q}B^{q}$
ARIMA模型的实质就是差分运算与ARMA模型的组合
ARIMA模型的使用场合：差分平稳序列的建模问题

ARIMA模型的B-J建模方法

ARIMA模型的Box-Jenkins建模思想

Box-Jenkins建模思想可分为如下4个步骤：

对原序列进行平稳性检验，如果序列不平稳，通过差分变换或者其他变换，使序列平稳；
通过相关系数法或信息准则法，确定ARMA模型的阶数p和q；
估计模型的未知参数，并检验参数的显著性，以及模型本身合理性
进行诊断分析，残差检验。

建模过程中的检验

建模过程中需要一些必要的检验工作：

时间序列数据的平稳性检验
残差序列的白噪声检验
模型参数的显著性检验（参数对应的t统计量及相伴概率）

时间序列数据的平稳性检验

采用单位根检验，可基于两种方法实现：

Augmented Dickey-Fuller(ADF)检验
Phillips-Perron（PP）检验

[h, p] = adftest(X, 'Alpha', alpha)
[h, p] = pptest(X, 'Alpha', alpha)

在置信度alpha下，检验时序数据X是否趋势平稳
（原假设Ho：有单位根，不平稳）
h=0表示接受原假设，即X不平稳；
h=1表示拒绝原假设，即可以认为×平稳；
p为统计量的相伴概率。
$p<\alpha,\qquad h=1;$

残差序列的白噪声检验

主要检验残差序列是否为白噪声序列，是否存在自相关性。这里介绍Ljung-Box Q方法：

[h, p] = lbqtest(res, 'Alpha', alpha)

在置信度alpha下，检验残差序列res是否存在自相关性
（原假设Ho：不存在自相关性，是白噪声序列）
h=0表示接受原假设，即认为res是白噪声序列；
h=1表示拒绝原假设，即认为res不是白噪声序列；
p为统计量的相伴概率。
$p<\alpha,\qquad h=1;$

ARIMA模型建模实例

![[Pasted image 20240821111648.png]]

平稳性检验

Matlab代码(ARIMASample.m)

load data_CPI
y = dataCPI;
T = length(y);  %序列长度%%（1）平稳性检验
% 可视化分析法
figure(1), plot(1:T, y)

![[Pasted image 20240821113036.png]]

%单位根检验法
[h,p] = adftest(y);

h = 0
p = 0.999
接受原假设，y 序列非平稳

差分运算

%%(2）差分运算
diff_y = diff(y);%一阶差分运算 
figure(2), plot(diff_y)

![[Pasted image 20240821113719.png]]

[h, p] = adftest（diff_y); %平稳性adf检验

h = 1
p = 0.0241
拒绝原假设，可以认为，diff_y序列平稳

模型定阶

%% (3) 模型定阶
figure(3), autocorr（diff_y); %自相关系数图
figure(4), parcorr  (diff_y); %偏自相关系数图

![[Pasted image 20240821114044.png]]

自相关图中1，2，3，4阶显著不为0
偏自相系数中1，2阶显著不为0
difff_y序列可以尝试建立：AR(2)、MA(4)、ARMA(2.4)等模型
原序列可以尝试建立：ARIMA(0,1,4)、ARIMA(2,1.0)、ARIMA(2,1,4)等模型

ARIMA(2,1,0)

拟合ARIMA模型

%%（4）拟合ARIMA(2.1,0)模型 
mdl = arima(2, 1, 0);
EstMdl210 = estimate(mdl, y);

ARIMA(2,1,0) Model:
Value TStat PValue
Constant 0.0101 3.071 0.002
AR{1} 0.212 2.222 0.026
AR{2} 0.337 3.250 0.001
Variance 9.230e-05 8.307 9.850e-17

模型参数显著性检验

%%（5）模型参数显著性检验
%观察t统计量以及相伴概率

模型参数显著性检验：
模型参数AR{1}，AR{2}检验的t统计量相伴概率p值都小于0.05，模型参数显著。

残差的白噪声检验

%%（6）残差的白噪声检验
% 检验残差是否为白噪声序列，是否存在相关性
res = infer(EstMdl, y); %残差序列
%可视化分析 
figure(5),
subplot(2,2,1), plot(res./sqrt(EstMdl.Variance))
subplot(2,2,2), qqplot(res) 
subplot(2,2,3), autocorr(res) 
subplot(2,2,4), parcorr(res)

![[Pasted image 20240821135157.png]]

res序列是白噪声序列，不存在自相关性

模型白噪声检验

%%（6）模型白噪声检验
% 检验残差是否为白噪声序列，是否存在自相关性
% Ljung-Box Q (lbq) 检验
[h,p] = lbqtest(res)

h = 0
p = 0.28
接受原假设，即可以认为res序列是白噪声序列，不存在自相关性
所建立的ARIMA(2,1,0)模型，较好地实现了原数据序列规律的学习，可以用于进一步预测分析

ARIMA(2,1,4)

%%（4）拟合ARIMA（2.1,4）模型
mdl = arima(2, 1, 4);
EstMdl214 = estimate(mdl, y);%%（5）模型参数显著性检验
%观察t统计量以及相伴概率

![[Pasted image 20240821140040.png]]

模型AR(1)，MA(3)不显著

剔除AR(1)

mdl = arima('ARLags', [2], ...
'D', 1, ...
'MALags',[1:4]);
EstMdl214n1 = estimate(mdl, y);

![[Pasted image 20240821140238.png]]

模型参数显著性检验：模型MA{2} MA{3}的系数不显著

剔除MA(3)

mdl = arima('ARLags', [2], ...
'D', 1, ...
'MALags',[1,2,4]);
EstMdl214n2 = estimate(mdl,y);

![[Pasted image 20240821140522.png]]

模型参数都是显著的

模型合理性检验

%%模型合理性检验
%检验残差是否为白噪声序列，是否存在自相关性
res214n2 = infer（EstMdl214n2, y); %残差序列%Ljung-BoxQ(lbq)检验
[h, p] = lbqtest(res214n2)

h =0
p = 0.8720
接受原假设，即可以认为res214n2序列是白噪声序列，不存在自相关性
修正后的ARIMA(2,1,4)模型，也较好地实现了原数据序列规律的学习，可以用于进一步预测分析

基于AIC，BIC信息准则分析模型优劣

%%AIC/BIC的计算
mdl = arima(2, 1, 0);
[EstMdl, EstParamCov, LogL, info] = estimate(mdl, y);
[AICValue, BICValue] = aicbic(LogL, length(info.X), length(y))

EstMdl，训练好的模型
EstParamCov，带估计的所有参数的协方差矩阵
LogL，模型的最大自然函数
info，模型的结构变量
info.X，表示所有的待定参数组成的向量
length(info.X)，待定参数的个数
length(y)，样本个数n
结果：
AICValue = -488.8444
BICValue = -479.4691

结果：

AIC准则下
ARIMA（2,1,4）优于其它模型
BIC准则下：
ARIMA(2,1,0) ARIMA([2],1,4)
ARIMA([2],1,[1,2,4]) 优于其它模型
以上模型都是可接受的

模型预测

%%（7）模型预测
[yF, yMSE] = forecast(EstMdl210, 20, 'Y0', y); 
UB = yF + 1.96*sqrt(yMSE); %95%置信区间上限 
LB = yF - 1.96*sqrt(yMSE; %95%置信区间下限 
figure(6), hold on, plot(y, 'Color', [.75,.75,.75]); plot(T+1 : T+20, yF, 'r');
plot(T+1 : T+20, UB, 'k--', T+1 : T+20, LB, 'k--')
legend（'CPI'，预测值'，'95%置信区间） 
title（'CPI预测值'）

forecast，模型预测
EstMdl210，训练好的模型
20，向后预测的步数
Y0, y，指定时间序列模型最初的初始化序列数据
yF，预测值
yMSE，与每一个预测值对应的MSE的值，误差的平方的平均值

SARIMA模型及建模实例

季节波动时间序列

在某些时间序列中，由于季节性变化（包括季度、月度、周度等变化）或其他一些固有因素的变化，会存在一些明显的周期性，这类序列称为季节性序列
描述这类序列的模型之一是季节时间序列模型（seasonal ARIMA model)，用SARIMA表示

季节时间序列的重要特征

处理季节性时间序列的一个重要工具：

季节差分(s步差分)：可消除周期性变化
$\nabla_{s}X_{t}=(1-B^{s})X_{t}=X_{t}-X_{t-s}$
对于非平稳季节性时间序列，有时需要进行D次季节差分之后才能转换为平稳的序列。
$\nabla_{s}^{D}X_{t}=(1-B^{s})^{D}X_{t}$

季节时间序列(SARIMA)模型

随机季节模型
乘积季节模型
简单季节模型(一种相对简单的情形)
简单季节模型是指序列中的季节效应和其它效应之间是加法关系
通过简单的趋势差分、季节差分之后序列即可转化为平稳
模型形式可表示为：
$\Phi(B)\nabla_{s}\nabla^{d}X_{t}=\Theta(B)\varepsilon_{t}$

例

![[Pasted image 20240821152640.png]]

平稳性检验

Matlab代码(SARIMASample) 
load data_GDP;
y = yGDP_1992_2007';
T = length(y); %序列长度%%（1）可视化分析（平稳性检验） 
figure(1)
plot(y, 'linewidth', 3)
title('1992-2017年中国季度GDP')

![[Pasted image 20240821152734.png]]

该数据有明显的增长趋势和季节性，季节周期长度为4
2. 差分运算

%%(2）差分运算
%增长趋势差分
diffD1_y = diff(y); 
figure(2),
plot(diffD1_y, '-o', 'linewidth', 1.5)
title('季度GDP一阶差分序列')

![[Pasted image 20240821161738.png]]

%季节差分
diffD1S1_y = diffD1_y(5 : end) - diffD1_y(1 : end-4); 
figure(3),
plot(diffD1S1_y, '-o', 'linewidth', 1.5) 
title('季度GDP趋势和季节差分序列')

![[Pasted image 20240821161803.png]]

平稳性检验

%% （3）平稳性检验
[h, p] = adftest(diffD1S1_y);  %adf单位根检验

h = 1
p = 1.0000e-03
拒绝原假设，可以认为diffD1S1_y序列平稳
4. 模型定阶

%%（4）模型定阶（准则函数法） 
maxLags = 4;
AICSet = zeros(maxLags, maxLags); 
BICSet = zeros(maxLags, maxLags);
for p = 1:maxlags 
for q = 1:maxLagsmdl = arima('ARLags', [1:p], 'MALags', [1:q]);[EstMdl, EstParamCov, LogL, info] = estimate(mdl, diffD1S1_y);[AICSet(p, q), BICSet(p, q)] = aicbic(LogL, length(info.X), length(diffD1S1_y));
end 
end
figure(4), heatmap(AlCSet/1000) 
figure(5), heatmap(BICSet/1000)

![[Pasted image 20240821163325.png]]
![[Pasted image 20240821163337.png]]

原序列x可用SARIMA(4,1,3)建模
5. 拟合SARIMA(4,1,3)模型

%% (5)拟合SARIMA(4,1.,)模型
mdl = arima('Seasonality', 4,...
'D', 1,...
'ARLags', 1:4,... 
'MALags', 1:3)
[EstMdl] = estimate(mdl, y);

![[Pasted image 20240821163529.png]]

SARIMA(4,1,3)模型具体形式
$\Phi(B)\nabla_{4}\nabla^{1}X_{t}=c+\Theta(B)\varepsilon_{t}$
其中
$c = 182.23$
$\Phi(B)=1-0.31B-0.30B^{2}-0.27B^{3}+0.69B^{4}$
$\Theta(B)=1-0.68B-1.02B^{2}-0.65B^{3}$
6. 残差白噪声检验

%%（6）残差白噪声检验
res = infer(EstMdl,y);%残差序列%Ljung-BoxQ(lbq）检验
[h,p] =lbqtest(res)

h=0
p = 0.4210
接受原假设，即可以认为res序列是白噪声序列，不存在自相关性
可以认为所建SARIMA模型实现了原数据序列总体规律的学习，可用于进一步的预测分析
7. 模型预测

[yF, yMSE] = forecast(EstMdl, 20, 'Y0', y');
UB = yF + 1.96*sqrt(yMSE); %95%置信区间上限 
LB = yF - 1.96*sqrt(yMSE); %95%置信区间下限 
figure(6),hold on
plot(y, 'Color', [.75, .75, .75]); 
plot(T+1 : T+20, yF, 'r');
plot(T+1 : T+20, UB, 'k:', T+1 : T+20, LB, 'k:')
legend('GDP', '预测值','95%置信区间'） 
title('GDP预测')

![[Pasted image 20240821164545.png]]