【数学建模入门】时间序列｜山猫数量预测

山猫数量预测

某地区山猫的数量在连续114年的统计数据如表下所示。分析该数据得出山猫的生长规律，并预测以后两个年度山猫的数量。

数据的文本格式

269	321	585	871	1475	2821	3928	5943	4950	2577	523	98
184	279	409	2285	2685	3409	1824	409	151	45	68	213
546	1033	2129	2536	957	361	377	225	360	731	1638	2725
2871	2119	684	299	236	245	552	1623	3311	6721	4254	687
255	473	358	784	1594	1676	2251	1426	756	299	201	229
469	736	2042	2811	4431	2511	389	73	39	49	59	188
377	1292	4031	3495	537	105	153	387	758	1307	3465	6991
6313	3794	1836	345	382	808	1388	2713	3800	309	2985	3790
674	71	80	108	229	399	1132	2432	3575	2935	1537	529
485	662	1000	1520	2657	3396						

ARIMA vs ARMAX

ARIMA(Autoregressive Integrated Moving Average Model):自回归积分滑动平均模型，用于非平稳序列模型。

模型的原理为：先将非平稳序列经过d 阶差分得到平稳序列，对得到的平稳序列进行ARMA的模型预测。

对于一个非平稳序列模型，若必须对时间序列做 d 阶差分才能得到一个平稳序列，则使用ARIMA(p,d,q)模型，其中 d 是差分的阶数，p为自回归项数，q为移动平均项数。

---

ARMAX(Autoregressive Moving Average with Extra Input):带额外输入（可理解为外部干扰项）的自回归移动平均模型，模型结构为

矩阵表示可简化为：

其中，y(t)——随时间的输出值

       na——系统极点个数

       nb——系统零点个数 + 1

       nc——系数C的个数

       nk——在输入影响输出之前（也称系统的死区时间）发生的输入样本数

       y(t−1)…y(t−na) ——当前输出所依赖的先前的一系列输出

       u(t−nk)…u(t−nk−nb+1)——当前输出所依赖的先前和延迟的一系列输入

       e(t−1)…e(t−nc)——白噪声干扰值

  参数na、nb和nc是ARMAX模型的阶数，nk是延迟值。q是延迟因子，q的几次方作用在yt上，yt就倒退几个时间，如q^3* yt= y( t-3)。其中，

流程图

数据处理

clc,clear
a=textread('data84.txt');%把原始数据保存到纯文本文件data84.txt
%进行数据处理
a=a';
a=nonzeros(a);
n=length(a);
figure,plot(a,'.-');

选择ARMA模型

figure,subplot(2,1,1),autocorr(a); %自相关函数图像
subplot(2,1,2),parcorr(a); %偏相关函数图像

拖尾：始终有非零取值，不会在k大于某个常数后就恒等于零(或在0附近随机波动)

截尾：在大于某个常数k后快速趋于0为k阶截尾

1.截尾和拖尾判断
2.截尾和拖尾判断
3.截尾和拖尾判断

注意置信上下界对拖尾和截尾判断的影响

自相关函数和偏相关函数曲线。

通过自相关函数和偏相关函数的曲线判断使用ARMA模型。

可以看到数量成周期性变化，故选择ARIMA模型。

差分运算

for i = 11:n
    b(i-10) = a(i)-a(i-10);
end
b=b';
figure,plot(b,'.-');

AIC BIC 准则确定p,q

%AIC BIC 准则进行定阶
k=0; %初始化试探模型的个数
for i = 0:3
    for j = 0:3
        if i==0 && j==0
            continue
        elseif i==0
            ToEstMd=arima('MALags',1:j,'Constant',0); %指定模型的结构
        elseif j==0
            ToEstMd=arima('ARLags',1:i,'Constant',0); %指定模型的结构
        else
            ToEstMd=arima('ARLags',1:i,'MALags',1:j,'Constant',0);%指定模型结构
        end
        k=k+1;R(k)=i;M(k)=j;
        [EstMd,EstParamCov,logL,info]=estimate(ToEstMd,w);
        numParams=sum(any(EstParamCov)); %计算拟合参数的个数
        %AIC BIC
        [aic(k),bic(k)]=aicbic(logL,numParams,n2);
    end
end
fprintf('R,M,AIC,BIC的对应值如下\n % f');
check=[R',M',aic',bic']

p=1,q=1

AIC、BIC越小越好

注意差分还原情况

最后预测数据：

---

参考资料：
b站讲解资料
ARMA模型参数估计——增广最小二乘