function?main
clear;clf;clc;
load?data_wind?wind;y=chafeng（wind2）;%x做一次差分得y
[cqcs]=fangcha（y）;%求y的自，偏自相關系數存于cqcs中并畫他們的序列圖
%選ARMA（4，0）模型；
%[p1p2]=ARMA（62y）;%求已知階數ARMA模型的參數
[p1p2]=ARMA（60y）;%求已知階數AR模型的參數
%[p1p2]=ARMA（01y）;%求已知階數MA模型的參數
result=[];N=size（y2）;%p1=[-0.0172?0.2932?0.2161?0.1603];
n=size（p11）;
for?k=1:n??????????????????????????????%%%%驗證模型對y的符合度
????result（k）=y（k）;?
end;
for?k=n+1:N
????sum1=0;
????for?i=1:n
????????sum1=sum1+p1（i）*y（k-i）;
????end;
????result（k）=sum1;
end;
for?k=1:10??????????????????????????????%%%%對y的預測
????sum1=0;
????for?i=1:n
????????sum1=sum1+p1（i）*result（k+N-i）;
????end;
????result（k+N）=sum1;
?
end;
?%?subplot（223）;plot（y）;hold?on;
?
??plot（result‘r‘）;title（‘一次差分‘）;%%%畫圖
result1（1）=wind（1）;
for?k=2:N???????????????%%%%回代得到x的預測，這是一次差分的逆過程
????result1（k）=wind（k-1）+result（k-1）;
end;
%for?k=N+1:N+10
?for?k=2:N
???result1（k）=result1（k-1）+result（k-1）???
end;
%subplot（224）;
plot（result1‘*-‘）;
hold?on;
plot（wind‘r+-‘）;
grid?on
xlabel（‘Time（h）‘）
ylabel（‘Vw（m/s）‘）
legend（‘預測‘‘實際‘）;
%title（‘擬合后‘）;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function?[faitheta]=ARMA（pqwind）
%%函數[faitheta]=ARMA（pqwind）
%%%pq:ARMA模型的階數，wind:原始序列值，
%%fai:對應p的參數，theta:對應q的參數

N=size（wind2）;
theta=[];
gama=[];G=[];g=[];
for?k=0:20%%%%%%%%%%%%%%%%%%%%%%%%%求fai
????sum=0;
????for?i=1:N-k
????????sum=sum+wind（i）*wind（i+k）;
????end;?
????gama（k+1）=sum/N;
end;
for?k=1:p
????g（k）=gama（q+k+1）;%%生成g向量
????for?i=1:p
????????G（ki）=gama（abs（q+k-i）+1）;%%%生成G矩陣
????end;
end;
fai=inv（G）*g‘;%%%%公式?5-83
%由于求MA模型參數時涉及到非線形方程組的求解，沒有統一的公式，本程序只就q=123時給出程序
y=[];?????????????????%%%求theta
k=0;i=0;
for?k=p+1:N
????sum=0;
????for?i=1:p
????????sum=sum+fai（i）*wind（k-i）;
????end;
????y（k-p）=wind（k）-sum;
end;
k=0;i=0;
gama1=[];
for?k=0:20
????sum=0;
????for?i=1:N-k-p
????????sum=sum+y（i）*y（i+k）;
????end;
????gama1（k+1）=sum/N;
end;
syms?a1?a2?a3?g1?g2?g3?g4?delta;
switch?q%%%%%分別根據q=123求非線性方程組???公式5-77
????case?1
????????s=solve（‘g1=delta*（1+a1）‘‘g2=delta*（-a1）‘deltaa1）;
????????theta（1）=subs（s.a1{g1gama1（1）}{g2gama1（2）}）;
????case?2
????????s=solve（‘g1=delta*（1+a1+a2）‘‘g2=delta*（-a1+a1*a2）‘‘g3=delta*（-a2）‘deltaa1a2）;
????????kk=subs（s.a1（1）{g1gama1（1）}{g2gama1（2）}{g3gama1（3）}）;
????????kk=subs（s.a2（1）{g1gama1（1）}{g2gama1（2）}{g3gama1（3）}）;
????case?3
????????s=solve（‘g1=delta*（1+a1+a2+a3）‘‘g2=delta*（-a1+a1*a2+a2*a3）‘‘g3=delta*（-a2+a2*a3）‘‘g4=delta*（-a3）‘deltaa1a2a3）;
????????theta（1）=subs（s.a1{g1gama1（1）}{g2gama1（2）}{g3gama1（3）}{g4gama1（4）}）;
????????theta（2）=subs（s.a2{g1gama1（1）}{g2gama1（2）}{g3gama1（3）}{g4gama1（4）}）;
????????theta（3）=subs（s.a3{g1gama1（1）}{g2gama1（2）}{g3gama1（3）}{g4gama1（4）}）;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function?y=chafeng（windN）%%函數差分
%%%%x:原始序列，N：差分的次數
%%%y:差分后的序列
p=[];
i=size（wind2）;
for?k=N+1:i
????p（k-N）=wind（k）-wind（k-N）;
end;
y=p;
%%%%%%%%

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----

?????文件???????4311??2009-05-05?03:20??time_series.m

-----------?---------??----------?-----??----

?????????????????4311????????????????????1