ARIMA模型

一、ARIMA模型介绍

ARIMA模型全称为自回归积分滑动平均模型(Autoregressive Integrated Moving Average Model,简记ARIMA)，是由博克思(Box)和詹金斯(Jenkins)于70年代初提出一著名时间序列预测方法^[1]  ，所以又称为box-jenkins模型、博克思-詹金斯法。其中ARIMA（p，d，q）称为差分自回归移动平均模型，AR是自回归， p为自回归项； MA为移动平均，q为移动平均项数，d为时间序列成为平稳时所做的差分次数。所谓ARIMA模型，是指将非平稳时间序列转化为平稳时间序列，然后将因变量仅对它的滞后值以及随机误差项的现值和滞后值进行回归所建立的模型。ARIMA模型根据原序列是否平稳以及回归中所含部分的不同，包括移动平均过程（MA）、自回归过程（AR）、自回归移动平均过程（ARMA）以及ARIMA过程。

ARIMA模型的基本思想是：将预测对象随时间推移而形成的数据序列视为一个随机序列，用一定的数学模型来近似描述这个序列。这个模型一旦被识别后就可以从时间序列的过去值及现在值来预测未来值。

二、ARIMA模型建模过程

1. 检查平稳性

平稳性就是围绕着一个常数上下波动且波动范围有限，即有常数均值和常数方差。如果有明显的趋势或周期性，那它通常不是平稳序列。 

不平稳序列可以通过差分转换为平稳序列。d阶差分就是相距d期的两个序列值之间相减。如果一个时间序列经过差分运算后具有平稳性，则该序列为差分平稳序列，可以使用ARIMA模型进行分析。 

2、确定模型阶数

AIC准则：即最小信息准则，同时给出ARMA模型阶数和参数的最佳估计，适用于样本数据较少的问题。目的是判断目标的发展过程与哪一个随机过程最为接近。因为只有样本量足够大时，样本的自相关函数才非常接近原时间序列的自相关函数。具体运用时，在规定范围内使模型阶数由低到高，分别计算AIC值，最后确定使其值最小的阶数，就是模型的合适阶数。

模型参数的最大似然估计时：AIC =(n-d)loge + 2(p+q+2)

模型参数的最小二乘估计时：AIC =(n-d)loge + (p+q+1)log n

clc;
x=csvread('E:\某事件.csv',1,4);
Data=x;

len = length(Data)
len_t = 100;
step=len - len_t;

SourceData=Data(1:len_t); 
TempData=SourceData;
TempData=detrend(TempData);%去趋势线
TrendData=SourceData-TempData;%趋势函数
%--------差分，平稳化时间序列---------

H=adftest(SourceData)
difftime=0;
SaveDiffData=[];
while ~H
SaveDiffData=[SaveDiffData,TempData(1,1)];
TempData=diff(TempData);%差分，平稳化时间序列
difftime=difftime+1;%差分次数
H=adftest(TempData);%adf检验，判断时间序列是否平稳化
end
%---------模型定阶或识别--------------

u = iddata(TempData);

test = [];
for p = 1:24                       %自回归对应PACF,给定滞后长度上限p和q，一般取为T/10、ln(T)或T^(1/2),这里取T/10=12
for q = 1:5                    %移动平均对应ACF
m = armax(u,[p q]);        
AIC = aic(m);              %armax(p,q),计算AIC
test = [test;p q AIC];
end
end
test;
size(test);

for k = 1:size(test,1)
if test(k,3) == min(test(:,3)) %选择AIC值最小的模型
p_test = test(k,1);
q_test = test(k,2);
break;
end
end

%------1阶预测-----------------
TempData=[TempData;zeros(step,1)];
m = armax(u,[p_test q_test]);        %m = armax(u(1:ls),[p_test q_test]);        %armax(p,q),[p_test q_test]对应AIC值最小，自动回归滑动平均模型
data1 = iddata(SourceData);
P1 = predict(m,data1,1);
PreR1 = P1.OutputData;
PreR1 = PreR1';

result = []
for k = len_t:len-1
    TempData = Data(1:k);
    n = iddata(TempData);
    P2 = forecast(m,n,1);
    PreR2 = P2.OutputData;
    PreR2 = PreR2';
    result = [result, PreR2];
end
%----------还原差分-----------------

if size(SaveDiffData,2)~=0
for index=size(SaveDiffData,2):-1:1
PreR=cumsum([SaveDiffData(index),PreR]);
end
end
 
  %-------------------预测趋势并返回结果----------------

mp1=polyfit([1:size(TrendData',2)],TrendData',1);
xt=[];
for j=1:step
xt=[xt,size(TrendData',2)+j];
end
TrendResult=polyval(mp1,xt);
PreData=result+TrendResult;
ttt = PreData - Data(len_t+1,:);
error = sum(abs(ttt))/step
result = [PreR1, result];
tempx =[TrendData',TrendResult] + result;% tempx为预测结果
tempx = tempx';
plot(Data,'b');
hold on
plot(tempx,'r');
print(gcf,'-dpng','E:/pre901/name.png')