%% 根据后面网址程序进行修改后,可以预测。https://www.ilovematlab.cn/thread-494019-1-1.html
clear;
close all;
plength=1;
flength=10;
y= xlsread('oxygen.xlsx','I3:I103'); %%43 自己可以随便建一个数据
Dedata=detrend(y); %去趋势
Tredata=y-Dedata;
figure;
plot(y,'r');
hold on;
plot(Dedata,'g');
hold on;
plot(Tredata,'y');
grid
legend('orignal data','detrend Data ','trend data')
[r,c]=size(Dedata);
figure;
subplot(2,1,1)
autocorr(Dedata); %原序列的自相关函数图MA(q),观察系数是否在区间(-2T^(1/2),-2T^(1/2))内
subplot(2,1,2)
parcorr(Dedata); %原序列的偏相关函数图AR(p),观察系数是否在区间(-2T^(1/2),-2T^(1/2))内
%% 如果该序列不是平稳的做差分图,否则跳过该步
DX=Dedata;
d=0;
H=adftest(DX);
SaveDiffData=[];
while ~H %% H=1 表示平稳
SaveDiffData=[SaveDiffData;DX(1,1)];
DX=diff(DX);
d=d+1;
H=adftest(DX);
end
%% 计算差分序列MA和AR
figure;
subplot(2,1,1)
autocorr(DX); %差分序列DX自相关函数图MA(q),观察系数是否在区间(-2T^(1/2),-2T^(1/2))内
subplot(2,1,2)
parcorr(DX); %差分序列DX偏相关函数图AR(p),观察系数是否在区间(-2T^(1/2),-2T^(1/2))内
figure;
plot(DX); %差分序列DX自相关函数图MA(q),观察系数是否在区间(-2T^(1/2),-2T^(1/2))内
title('Dx');
%% 对差分后的序列做拟合和预测,求出最好的阶数
DX=reshape(DX,r-d,1);
DX_data=iddata(DX); %将DX转化为matlab接受的格式
test = [];
%p = [0 1 2 3];%自回归对应PACF,给定滞后长度上限p和q,一般取为T/10、ln(T)或T^(1/2),这里取
%q = [0 1 2 3];%移动平均对应ACF
for p = 1:10 %自回归对应PACF, p=3 原始的 给定滞后长度上限p和q,一般取为T/10、ln(T)或T^(1/2),这里取T/10=12
for q = 1:10 %移动平均对应ACF q=3是原始给的
m = armax(DX_data,[p,q]);
AIC = aic(m); %armax(p,q),选择对应FPE最小,AIC值最小的模型
test = [test;p q AIC];
end
end
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
%% 拟合过程
sys = armax(DX_data,[p_test q_test]);%armax(p,q),[p_test q_test]对应AIC值最小
figure;
e = resid(sys,DX_data); %拟合做残差分析
plot(e);
grid;
title('Residual analysis');
%% 检验残差的自相关和偏相关函数 是否是平稳序列
figure;
subplot(2,1,1);
autocorr(e.OutputData); %d阶差分序列z自相关函数图MA(q),置信水平0.95
subplot(2,1,2);
parcorr(e.OutputData); %d阶差分序列z偏相关函数图AR(p),置信水平0.95
%% 预测过程
p=predict(sys,DX_data,plength);%sys是模型 data是真实数据 plength预测长度
f=forecast(sys,DX_data,flength);
figure;
compare(DX,p,'g');
figure;
plot(DX,'-');
hold on;
plot(p,'-<');
hold on;
plot(f,':*');
grid
legend('DX data','predict Data ','forecast data')
%% 差分还原
PreP=p.OutputData;
PreF=f.OutputData;
PrePF=[PreP;PreF];
row=size(SaveDiffData,1);
if row~=0 % size(SaveDiffData,2)返回列数 1 返回行数
for index=row:-1:1
PreP=cumsum([SaveDiffData(index,1);PreP]);
PrePF=cumsum([SaveDiffData(index,1);PrePF]);
end
end
PreP=PreP+Tredata;
Taver=mean(Tredata);
Taver_matr=Taver*ones(flength,1);%%%趋势没有预测,只是对其求了个平均值,后期都是用平均值代替forecast趋势
Tredata=[Tredata;Taver_matr];
PrePF=PrePF+Tredata; %
figure;
plot(y,'-');
hold on;
% plot(PreP,'-<');
% hold on;
plot(PrePF,'-*');
grid
legend('original data','forecast data')
- 1
- 2
- 3
前往页