时间序列的ARIMA预测分析

测得某一地区一口井7年的地下水埋深数据，试预报第8年全年的地下水埋深。

 

1,平稳性检验

画原始序列图和自相关图

x=[9.40,8.81,8.65,10.01,11.07,11.54,12.73,12.43,11.64,11.39,11.10,10.85,...
    10.71,10.24,8.48,9.88,10.31,10.53,9.55,6.51,7.75,7.80,5.96,5.21,...
    6.39,6.38,6.51,7.14,7.26,8.49,9.39,9.71,9.65,9.26,8.84,8.29,...
    7.21,6.93,7.21,7.82,8.57,9.59,8.77,8.61,8.94,8.40,8.35,7.95,...
    7.66,7.68,7.85,8.53,9.38,10.09,10.59,10.83,10.49,9.21,8.66,8.39,...
    8.27,8.14,8.71,10.43,11.47,11.73,11.61,11.93,11.55,11.35,11.11,10.49,...
    10.16,9.96,10.47,11.70,10.10,10.37,12.47,11.91,10.83,10.64,10.29,10.34];
x=reshape(x',numel(x),1);
subplot(1,2,1);
plot(x)
title('原始数据图像');
subplot(1,2,2);
autocorr(x)
title('自相关函数图像');

画1阶12步差分序列图和自相关图

s=12;%周期
n=12;%预报数据个数
m1=length(x);
for i=s+1:m1
    y(i-s)=x(i)-x(i-s);
end
x1=diff(y);
m2=length(x1);
figure;
subplot(1,2,1);
plot(x1)
title('查分后序列图像');
subplot(1,2,2);
autocorr(x1)
title('自相关函数图像');

对差分后的序列进行平稳性检验

>> [h1,p1,adf,ljz]=adftest(x1)

h1 =

  logical

   1


p1 =

   1.0000e-03


adf =

   -7.4926


ljz =

   -1.9450

表示1阶12步差分后的数据是平稳的。

 

2.白噪声检验

yanchi=[6,12,18];
[H,pValue,Qstat,CriticalValue]=lbqtest(x1,'lags',yanchi);
fprintf('%15s%15s%15s','延迟阶数','卡方统计量','p值');
fprintf('\n');
for i=1:length(yanchi)
    fprintf('%18f%19f%19f',yanchi(i),Qstat(i),pValue(i));
    fprintf('\n');
end

           延迟阶数          卡方统计量             p值
          6.000000          11.831682           0.065831
         12.000000          38.603461           0.000122
         18.000000          46.794167           0.000227

p值小于0.1，说明数据不是白噪声序列，可以建立ARIMA模型。

 

3.模型识别（AIC\BIC准则）

x1=x1';
LOGL=zeros(6,6);
PQ=zeros(6,6);
for p=0:5
    for q=0:5
        model=arima(p,0,q);
        [fit,~,logL]=estimate(model,x1);  %指定模型的结构
        LOGL(p+1,q+1)=logL;
        PQ(p+1,q+1)=p+q;  %计算拟合参数的个数
    end
end
LOGL=reshape(LOGL,36,1);
PQ=reshape(PQ,36,1);
[aic,bic]=aicbic(LOGL,PQ+1,m2);
fprintf('AIC为：\n');
reshape(aic,6,6)
fprintf('BIC为：\n');
reshape(bic,6,6)

AIC为：

ans =

  210.6009  210.6412  206.1356  207.7248  208.0342  206.3506
  211.8987  203.2963  205.0547  208.4185  209.4344  205.0287
  205.7454  206.9559  208.8862  200.6190  207.4144  203.4940
  206.5786  208.2715  206.8524  196.0960  198.6505  204.3130
  208.4827  210.1279  200.9564  197.4272  197.1011  199.0800
  207.5084  209.4810  197.2431  199.0064  196.3561  186.3381

BIC为：

ans =

  212.8636  215.1666  212.9237  216.7755  219.3476  219.9266
  216.4240  210.0844  214.1055  219.7319  223.0105  220.8675
  212.5335  216.0066  220.1996  214.1951  223.2531  221.5954
  215.6293  219.5849  220.4284  211.9347  216.7519  224.6771
  219.7961  223.7040  216.7951  215.5286  217.4652  221.7068
  221.0844  225.3197  215.3446  219.3705  218.9829  211.2275

寻找使AIC\BIC值最小的p值和q值（p、q值越小越好）

得p=q=1

 

4.模型估计

p=input('输入阶数p=');
q=input('输入阶数q=');
model=arima(p,0,q);  %指定模型的结构
m=estimate(model,x1);  %拟合参数

输入阶数p=1
输入阶数q=1
 
    ARIMA(1,0,1) Model (Gaussian Distribution):
 
                 Value      StandardError    TStatistic      PValue  
                ________    _____________    __________    __________

    Constant    -0.12332       0.21876        -0.56371        0.57295
    AR{1}       -0.72233       0.10787         -6.6964     2.1366e-11
    MA{1}        0.99639      0.060558          16.454     7.9052e-61
    Variance     0.94265       0.12861          7.3298     2.3048e-13

可得模型结构为

 

5.模型检验

对残差序列进行白噪声检验

z=iddata(x1);
ml1=armax(z,[p,q]);
e=resid(ml1,z);  %拟合做残差分析
[H,pValue,Qstat,CriticalValue]=lbqtest(e.OutputData,'lags',6)

H =

  logical

   0


pValue =

    0.3979


Qstat =

    6.2303


CriticalValue =

   12.5916

说明残差是白噪声序列。

 

6.运用模型进行预测

[yf,ymse]=forecast(m,n,'Y0',x1);
yhat=y(end)+cumsum(yf);  %求y的预报值
for j=1:n
    x(m1+j)=yhat(j)+x(m1+j-s);  %求x的预测值
end
xhat=x(m1+1:end)

xhat =

   10.3218
    9.7733
   10.4117
   11.4256
    9.8584
    9.9814
   12.0643
   11.3933
   10.2702
    9.9880
    9.5813
    9.5490

此即为预测的第8年全年的数据。