山猫数量预测
某地区山猫的数量在连续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阶截尾
注意置信上下界对拖尾和截尾判断的影响
自相关函数和偏相关函数曲线。
通过自相关函数和偏相关函数的曲线判断使用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
注意差分还原情况
最后预测数据:
参考资料:
b站讲解资料
ARMA模型参数估计——增广最小二乘