?? arma-model-for-oil-price.m
字號:
clear;
%--------------------------------油價序列零均值化后的數據如下----------------------------------------%：
P=[ 19.5900 14.9100 15.7400 15.4000 13.0600 19.0700 15.2800 15.8200 12.7700 12.0500...
11.6900 13.8500 13.8500 10.0700 9.1700 10.7900 13.4400 21.1700 18.6400 13.2100...
15.5400 21.9400 23.1100 18.6400 14.9400 16.9000 15.4600 11.1500 13.1300 12.4800...
12.9500 12.5900 10.5800 10.5800 12.3900 15.5300 13.0600 10.2200 16.3300 19.7200...
21.3100 18.8400 24.8400 15.6700 15.5700 12.7300 13.5600 15.5400 17.2200 12.1400...
11.0700 12.0200 11.5500 6.9200 10.3300 8.3800 12.1100 11.4600 12.7500 13.3200...
13.0000 11.9000 11.7900 12.5500 11.8400 11.2500 11.1500 10.9900 11.7000 14.0100...
17.5100 17.2700 16.9000 15.7900 15.4500 6.2400 16.7100 16.7700 16.6400 17.8000...
16.8700 16.1300 15.7600 15.6600 15.5400 15.3000 15.0500 14.6900 14.3900 14.1800...
13.70 13.66 13.27 13.56 13.14 14.19 ];
F=[ 19.5900 14.9100 15.7400 15.4000 13.0600 19.0700 15.2800 15.8200 12.7700 12.0500...
11.6900 13.8500 13.8500 10.0700 9.1700 10.7900 13.4400 21.1700 18.6400 13.2100...
15.5400 21.9400 23.1100 18.6400 14.9400 16.9000 15.4600 11.1500 13.1300 12.4800...
12.9500 12.5900 10.5800 10.5800 12.3900 15.5300 13.0600 10.2200 16.3300 19.7200...
21.3100 18.8400 24.8400 15.6700 15.5700 12.7300 13.5600 15.5400 17.2200 12.1400...
11.0700 12.0200 11.5500 6.9200 10.3300 8.3800 12.1100 11.4600 12.7500 13.3200...
13.0000 11.9000 11.7900 12.5500 11.8400 11.2500 11.1500 10.9900 11.7000 14.0100...
17.5100 17.2700 16.9000 15.7900 15.4500 6.2400 16.7100 16.7700 16.6400 17.8000...
16.8700 16.1300 15.7600 15.6600 15.5400 15.3000 15.0500 14.6900 14.3900 14.180];
%----------------------由于時間序列有不平穩(wěn)趨勢，進行兩次差分運算，消除趨勢性----------------------%
for i=2:96
Yt(i)=P(i)-P(i-1);
end
for i=3:96
L(i)=Yt(i)-Yt(i-1);
end
figure;
L=L(3:96);
Y=L(1:88);
plot(P);
title('原數據序列圖');
hold on;
pause
plot(Y,'r');
title('兩次差分后的序列圖和原數對比圖');
pause
%--------------------------------------對數據標準化處理----------------------------------------------%
Ux=sum(Y)/88 % 求序列均值
yt=Y-Ux;
b=0;
for i=1:88
b=yt(i)^2/88+b;
end
v=sqrt(b) % 求序列方差
Y=(Y-Ux)/v; % 標準化處理公式
f=F(1:88);
t=1:88;
figure;
plot(t,f,t,Y,'r')
title('原始數據和標準化處理后對比圖');
xlabel('時間t'),ylabel('油價y');
legend('原始數據 F ','標準化后數據Y ');
pause
%--------------------------------------對數據標準化處理----------------------------------------------%
%------------------------檢驗預處理后的數據是否符合AR建模要求，計算自相關和偏相關系數---------------%
%---------------------------------------計算自相關系數-----------------------------------%
R0=0;for i=1:88
R0=Y(i)^2/88+R0;
end
R0
for k=1:20
R(k)=0;
for i=k+1:88 R(k)=Y(i)*Y(i-k)/88+R(k); end
R %自協方差函數R
end
x=R/R0 %自相關系數x
figure;
plot(x)
title('自相關系數分析圖');
pause
%-----------------------------------計算自相關系數-------------------------------------%
%-----------------------解Y-W方程，其系數矩陣是Toeplit
矩陣。求得偏相關函數X-----------------------%
X1=x(1);
X11=x(1);
B=[x(1) x(2)]';
x2=[1 x(1)];
A=toeplitz(x2);
X2=A\B
X22=X2(2)
B=[x(1) x(2) x(3)]';
x3=[1 x(1) x(2)];
A=toeplitz(x3);
X3=A\B
X33=X3(3)
B=[x(1) x(2) x(3) x(4)]';
x4=[1 x(1) x(2) x(3)];
A=toeplitz(x4);
X4=A\B
X44=X4(4)
B=[x(1) x(2) x(3) x(4) x(5)]';
x5=[1 x(1) x(2) x(3) x(4)];
A=toeplitz(x5);
X5=A\B
X55=X5(5)
B=[x(1) x(2) x(3) x(4) x(5) x(6)]';
x6=[1 x(1) x(2) x(3) x(4) x(5)];
A=toeplitz(x6);
X6=A\B
X66=X6(6)
B=[x(1) x(2) x(3) x(4) x(5) x(6) x(7)]';
x7=[1 x(1) x(2) x(3) x(4) x(5) x(6)];
A=toeplitz(x7);
X7=A\B
X77=X7(7)
B=[x(1) x(2) x(3) x(4) x(5) x(6) x(7) x(8)]';
x8=[1 x(1) x(2) x(3) x(4) x(5) x(6) x(7)];
A=toeplitz(x8);
X8=A\B
X88=X8(8)
B=[x(1) x(2) x(3) x(4) x(5) x(6) x(7) x(8) x(9)]';
x9=[1 x(1) x(2) x(3) x(4) x(5) x(6) x(7) x(8)];
A=toeplitz(x9);
X9=A\B
X99=X9(9)
B=[x(1) x(2) x(3) x(4) x(5) x(6) x(7) x(8) x(9) x(10)]';
x10=[1 x(1) x(2) x(3) x(4) x(5) x(6) x(7) x(8) x(9)];
A=toeplitz(x10);
X10=A\B
X1010=X10(10)
B=[x(1) x(2) x(3) x(4) x(5) x(6) x(7) x(8) x(9) x(10) x(11)]';
x11=[1 x(1) x(2) x(3) x(4) x(5) x(6) x(7) x(8) x(9) x(10)];
A=toeplitz(x11);
X101=A\B
X1111=X101(11)
B=[x(1) x(2) x(3) x(4) x(5) x(6) x(7) x(8) x(9) x(10) x(11) x(12)]';
x12=[1 x(1) x(2) x(3) x(4) x(5) x(6) x(7) x(8) x(9) x(10) x(11)];
A=toeplitz(x12);
X12=A\B
X1212=X12(12)
X=[X11 X22 X33 X44 X55 X66 X77 X88 X99 X1010 X1111 X1212]
%-----------------------------------解Y-W方程，得偏相關函數X-------------------------------------%
figure;
plot(X);
title('偏相關函數圖');
pause
%-----根據偏相關函數截尾性，初判模型階次為5。用最小二乘法估計參數,計算10階以內的模型殘差方差和AIC值，應用AIC準則為模型定階------%
S=[R0 R(1) R(2) R(3) R(4)];
G=toeplitz(S);
W=inv(G)*[R(1:5)]' % 參數W(i) 與X5相同
K=0;
for t=6:88
r=0;
for i=1:5
r=W(i)*Y(t-i)+r;
end
at= Y(t)-r;
K=(at)^2+K;
end
U(5)=K/(88-5) % 5階模型殘差方差 0.4420
K=0;T=X1;
for t=2:88
at=Y(t)-T(1)*Y(t-1);
K=(at)^2+K;
end
U(1)=K/(89-1) % 1階模型殘差方差0.6954
K=0;T=X2;
for t=3:88
r=0;
for i=1:2
r=T(i)*Y(t-i)+r;
end
at= Y(t)-r;
K=(at)^2+K;
end
U(2)=K/(88-2) % 2階模型殘差方差 0.6264
K=0;T=X3;
for t=4:88
r=0;
for i=1:3
r=T(i)*Y(t-i)+r;
end
at= Y(t)-r;
K=(at)^2+K;
end
U(3)=K/(88-3) % 3階模型殘差方差 0.5327
K=0;T=X4;
for t=5:88
r=0;
for i=1:4
r=T(i)*Y(t-i)+r;
end
at= Y(t)-r;
K=(at)^2+K;
end
U(4)=K/(88-4) % 4階模型殘差方差 0.4751
K=0;T=X6;
for t=7:88
r=0;
for i=1:6
r=T(i)*Y(t-i)+r;
end
at= Y(t)-r;
K=(at)^2+K;
end
U(6)=K/(88-6) % 6階模型殘差方差 0.4365
K=0;T=X7;
for t=8:88
r=0;
for i=1:7
r=T(i)*Y(t-i)+r;
end
at= Y(t)-r;
K=(at)^2+K;
end
U(7)=K/(88-7) % 7階模型殘差方差 0.4331
K=0;T=X8;
for t=9:88
r=0;
for i=1:8
r=T(i)*Y(t-i)+r;
end
at= Y(t)-r;
K=(at)^2+K;
end
U(8)=K/(88-8) % 8階模型殘差方差0.4310
K=0;T=X9;
for t=10:88
r=0;
for i=1:9
r=T(i)*Y(t-i)+r;
end
at= Y(t)-r;
K=(at)^2+K;
end
U(9)=K/(88-9) %9階模型殘差方差 0.4297
K=0;T=X10;
for t=11:88
r=0;
for i=1:10
r=T(i)*Y(t-i)+r;
end
at= Y(t)-r;
K=(at)^2+K;
end
U(10)=K/(88-10) % 10階模型殘差方差 0.4317
U=10*U
for i=1:10
AIC2(i)=88*log(U(i))+2*(i) % AIC值分別為：172.6632 165.4660 153.2087 145.1442 140.7898 141.6824 142.9944 144.5601 146.3067 148.7036
end
%-----------------取使AIC值為最小值的階次，判斷模型階次為5。用最小二乘法估計參數--------------------%
%------------------檢驗{at}是否為白噪聲。求{at}的自相關系數，看其是否趨近于零-----------------------%
C=0;K=0;
for t=7:88
at=Y(t)-W(1)*Y(t-1)-W(2)*Y(t-2)-W(3)*Y(t-3)-W(4)*Y(t-4)-W(5)*Y(t-5)+Y(6)-W(1)*Y(5)-W(2)*Y(4)-W(3)*Y(3)-W(4)*Y(2)-W(5)*Y(1);
at1=Y(t-1)-W(1)*Y(t-2)-W(2)*Y(t-3)-W(3)*Y(t-4)-W(4)*Y(t-5)-W(5)*Y(t-6);
C=at*at1+C;
K=(at)^2+K;
end
p=C/K %若p接近于零，則{at}可看作是白噪聲
%--------------------------------{at}的自相關系數,趨近于零，模型適用--------------------------------%
%------------AR(5)模型方程為------------------------------------------------------------------------%
% X(t)=W(1)*X(t-1)-W(2)*X(t-2)-W(3)*X(t-3)-W(4)*X(t-4)-W(5)*X(t-5)+at (at=0.4420)
%------------------------------------------后六年的數據 進行預測和效果檢驗----------------------------------------------%
%-----------------------------單步預測 預測當前時刻后的六個數據----------------------------------%
XT=[L(84:94)];
for t=6:11
m(t)=0;
for i=1:5
m(t)=W(i)*XT(t-i)+m(t);
end
end
m=m(6:11);
%-------------預測值進行反處理---------------%
m(1)=Yt(90)+m(1); %一次反差分
z1(1)=P(90)+m(1); %二次反差分
m(2)=Yt(91)+m(2);
z1(2)=P(91)+m(2);
m(3)=Yt(92)+m(3);
z1(3)=P(92)+m(3);
m(4)=Yt(93)+m(4);
z1(4)=P(93)+m(4);
m(5)=Yt(94)+m(5);
z1(5)=P(94)+m(5);
m(6)=Yt(95)+m(6);
z1(6)=P(95)+m(6);
z1 % 單步預測的向后6個預測值:z1= 13.9423 13.4101 13.3588 12.9856 13.2594 12.9552
%---------------------------繪制數據模型逼近曲線-----------------------------------%
for t=6:88
r=0;
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -