测得某一地区一口井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:m1y(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 =logical1p1 =1.0000e-03adf =-7.4926ljz =-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.06583112.000000 38.603461 0.00012218.000000 46.794167 0.000227p值小于0.1，说明数据不是白噪声序列，可以建立ARIMA模型。
3.模型识别（AIC\BIC准则）
x1=x1';
LOGL=zeros(6,6);
PQ=zeros(6,6);
for p=0:5for q=0:5model=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.3506211.8987 203.2963 205.0547 208.4185 209.4344 205.0287205.7454 206.9559 208.8862 200.6190 207.4144 203.4940206.5786 208.2715 206.8524 196.0960 198.6505 204.3130208.4827 210.1279 200.9564 197.4272 197.1011 199.0800207.5084 209.4810 197.2431 199.0064 196.3561 186.3381BIC为：ans =212.8636 215.1666 212.9237 216.7755 219.3476 219.9266216.4240 210.0844 214.1055 219.7319 223.0105 220.8675212.5335 216.0066 220.1996 214.1951 223.2531 221.5954215.6293 219.5849 220.4284 211.9347 216.7519 224.6771219.7961 223.7040 216.7951 215.5286 217.4652 221.7068221.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=1ARIMA(1,0,1) Model (Gaussian Distribution):Value StandardError TStatistic PValue ________ _____________ __________ __________Constant -0.12332 0.21876 -0.56371 0.57295AR{1} -0.72233 0.10787 -6.6964 2.1366e-11MA{1} 0.99639 0.060558 16.454 7.9052e-61Variance 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 =logical0pValue =0.3979Qstat =6.2303CriticalValue =12.5916说明残差是白噪声序列。
6.运用模型进行预测
[yf,ymse]=forecast(m,n,'Y0',x1);
yhat=y(end)+cumsum(yf); %求y的预报值
for j=1:nx(m1+j)=yhat(j)+x(m1+j-s); %求x的预测值
end
xhat=x(m1+1:end)xhat =10.32189.773310.411711.42569.85849.981412.064311.393310.27029.98809.58139.5490此即为预测的第8年全年的数据。
