%AR時間序列模型估計

clear
tic
%%%%%%%%%%%%%%%第一步，模擬一個AR模型并繪制ACF，PACF圖%%%%%%%%%%%%%
%s首先設定AR模型的多項式系數。AR模型中只有多項式A（q）和C（q），
a1?=?0.6;
a2?=?0.8;
a3?=?0;
a4?=?0;
c1?=?0;
c2?=?0;
c3?=?0;
c4?=?0;?
obv?=?700;????????????????????????%obv是模擬的觀測數目。?
A?=?[1?a1?a2?a3?a4];
B?=?[];???????????????????????????%因為AR模型沒有輸入，因此多項式B是空的。
C?=?[1?c1?c2?c3?c4];
D?=?[];???????????????????????????%把D也設為空的。
F?=?[];???????????????????????????%AR模型里的F多項式也是空的。
m?=?idpoly（ABCDF11）????%這樣就生成了AR模型，把它存儲在m中。NoiseVariance被設定為1，1也是默認值。抽樣間隔Ts設為1。?
error?=?randn（obv1）;???????????%生成一個obv*1的正態隨機序列。準備用作模型的誤差項。
e?=?iddata（[]error1）;????????%用randn函數生成一個噪聲序列。存儲在e中。抽樣間隔是1秒。???????????????????????。
y?=?sim（me）;?
get（y）????????????????????????????%使用get函數來查看動態系統的所有性質。
Y=y.OutputData;?%把y.OutputData的全部值賦給變量r，r就是一個obv*1的向量。
dlmwrite（‘data.txt‘Y）?%將生成的700個數據保存在data.txt中
Y=Y‘
figure（1）
plot（Y）?????????????????%繪出y隨時間變化的曲線。?
title（‘時間序列圖像‘）
%%%%%%%%%%%%%%%%-檢驗預處理后的數據是否符合AR建模要求，計算自相關和偏相關系數%%%%%%%%%%%
R0=0;
for?i=1:500
??R0=Y（i）^2/500+R0;
end
R0
for?k=1:20
??R（k）=0;
??for?i=k+1:500
??R（k）=Y（i）*Y（i-k）/500+R（k）;
??end
??R?%自協方差函數R???
end
x=R/R0?%自相關系數x
figure;
plot（x）
title（‘自相關系數分析圖‘）;??
%%%%%%%%%%%%%%%%%%%%%%%%%%計算自相關系數%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%解Y-W方程，其系數矩陣是Toeplit矩陣。求得偏相關函數X%%%%%%%%%%%
X1=x（1）;
X11=x（1）/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）
X=[X11?X22?X33?X44?X55?]??
%%%%%%%%%%%%%%%%%%%%%%%解Y-W方程，得偏相關函數X%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure;??
plot（X）;
grid?on
title（‘偏相關函數圖‘）;

%%%%%%%%%%%%%%%根據偏相關函數截尾性，初判模型階次。%%%%%%%%%%%%%%%%%%%
%用最小二乘法估計參數計算5階以內的模型殘差方差和AIC值，應用AIC準則為模型定階%
??S=[R0?R（1）];
??G=toeplitz（S）;
??W=inv（G）*[R（1:2）]‘???????????
????
??K=0;
??for?t=3:498???
??r=0;??
??for?i=1:2
??r=W（i）*Y（t-i）+r;
??end
??at=?Y（t）-r;
??K=（at）^2+K;??
??end
??U（2）=K/（500-2）?%?2階模型殘差方差?

????
K=0;T=X1;
for?t=2:498
??at=Y（t）-T（1）*Y（t-1）;
??K=（at）^2+K;??
end???
??U（1）=K/（4