??Arma模型預(yù)測(cè)算法在兩年之前有看過，當(dāng)時(shí)沒有太仔細(xì)看沒能理解，最近結(jié)合網(wǎng)上幾篇比較Nice 的關(guān)于ARMA && ARIMA算法的博客，對(duì)該算法有了進(jìn)一步了解，將自己的理解進(jìn)行整理。
1 概述
??Arma模型（自回歸移動(dòng)平均模型）是時(shí)間序列分析中常用的模型之一，它可以用于預(yù)測(cè)未來的時(shí)間序列值。Arma模型的核心思想是將時(shí)間序列看作是自回歸和移動(dòng)平均過程的組合。其中，自回歸過程指的是時(shí)間序列值與其前一時(shí)刻值之間的關(guān)系；移動(dòng)平均過程指的是時(shí)間序列值與其前一時(shí)刻的噪聲誤差之間的關(guān)系。
Arma模型可以表示為ARMA(p, q)，其中p表示自回歸項(xiàng)的階數(shù)，q表示移動(dòng)平均項(xiàng)的階數(shù)。具體地，Arma模型可以寫成如下形式：
Y t = α 1 Y t ? 1 + α 2 Y t ? 2 + . . . + α p Y t ? p + ? t + β 1 ? t ? 1 + β 2 ? t ? 2 + . . . + β q ? t ? q Y_t=\alpha_1Y_{t-1}+\alpha_2Y_{t-2}+...+\alpha_pY_{t-p}+\epsilon_t+\beta_1\epsilon_{t-1}+\beta_2\epsilon_{t-2}+...+\beta_q\epsilon_{t-q} Yt?=α1?Yt?1?+α2?Yt?2?+...+αp?Yt?p?+?t?+β1??t?1?+β2??t?2?+...+βq??t?q?
??其中， Y t Y_t Yt?表示時(shí)間序列在時(shí)刻t的值， ? t \epsilon_t ?t?表示在時(shí)刻t的噪聲誤差， α 1 \alpha_1 α1?到 α p \alpha_p αp?和 β 1 \beta_1 β1?到 β q \beta_q βq?分別表示自回歸系數(shù)和移動(dòng)平均系數(shù)。
??Arma模型可以用于對(duì)時(shí)間序列進(jìn)行預(yù)測(cè)和模擬，通常使用最小二乘法或極大似然法來估計(jì)模型參數(shù)。同時(shí)，Arma模型也可以通過對(duì)模型殘差的分析來檢驗(yàn)?zāi)Ｐ偷臄M合程度和合理性。
??更加詳細(xì)的理論可以參考這篇博文
2 算法主要流程 及 測(cè)試方案
??練習(xí)時(shí)沒有使用真實(shí)時(shí)間序列，因?yàn)檎鎸?shí)序列的周期項(xiàng)、趨勢(shì)項(xiàng)中包含隨機(jī)誤差影響不可控，使不能夠直觀準(zhǔn)確反應(yīng)預(yù)測(cè)效果，這里使用模擬的時(shí)間序列進(jìn)行測(cè)試，考慮到一個(gè)完整的序列是由周期項(xiàng)、半周期項(xiàng)、趨勢(shì)項(xiàng)和隨機(jī)誤差組成，所以在模擬的時(shí)候?qū)@些因素都進(jìn)行了考慮；下面簡(jiǎn)單介紹下Arma的流程和測(cè)試的幾種方案下得到的的結(jié)果。
主要流程
- 加載時(shí)間序列，檢驗(yàn)序列是否平穩(wěn)，常見的平穩(wěn)性檢驗(yàn)方法是adf和kpss檢驗(yàn)，兩種檢驗(yàn)方法可能對(duì)同一序列有不同檢驗(yàn)結(jié)果;
%% matlab自帶
ad1 = adftest(a4); % 1平穩(wěn) 0非平穩(wěn)
ad2 = kpsstest(a4); % 0 平穩(wěn) 1 非平穩(wěn)
- 如果檢驗(yàn)顯示序列非平穩(wěn)，對(duì)序列進(jìn)行差分直至序列達(dá)到平穩(wěn)為止，一般差分1次就能滿足平穩(wěn)性要求；
%% 對(duì)時(shí)間序列進(jìn)行差分，并記錄差分了幾次
count = 0;
arr1 = a4;
while true
arr21 = diff(arr1);
count = count + 1;
ad1 = adftest(arr21); % 1平穩(wěn) 0非平穩(wěn)
ad2 = kpsstest(arr21); % 0 平穩(wěn) 1 非平穩(wěn)
if ad1==1 && ad2 == 0
break
end
arr1 = arr21;
end
- 模型定階：定階對(duì)預(yù)測(cè)結(jié)果還是比較重要的。參考博文，一般都是先看自相關(guān)(ACF) 和 偏相關(guān)(PACF)圖，根據(jù)拖尾和截尾可以人工選擇一個(gè)合適的階數(shù)；如果人工選擇比較困難，可以根據(jù)AIC、BIC準(zhǔn)則定階，具體理論可參考其他博文；
% 自相關(guān)和偏相關(guān)圖
figure(2)
subplot(2,1,1)
coll1 = autocorr(arr22);
stem(coll1)%繪制經(jīng)線圖
title('自相關(guān)')
subplot(2,1,2)
coll2 = parcorr(arr22);
stem(coll2)%繪制經(jīng)線圖
title('偏相關(guān)')
% Aic準(zhǔn)則
value = aic(AmMode);
- 序列預(yù)測(cè)：在預(yù)測(cè)的時(shí)候要考慮序列是否進(jìn)行過差分(diff)，因?yàn)榉瞧椒€(wěn)序列差分后進(jìn)行預(yù)測(cè)的，預(yù)測(cè)結(jié)果也是差分后的結(jié)果，此時(shí)要對(duì)序列還原 才能得到原序列的預(yù)測(cè)；
% 預(yù)測(cè)
aaa = predict(model, arr2, ntest); %利用arr2去向后預(yù)測(cè)
% 序列還原
hhh = cumsum([a4(1);difvalue]);%還原差分
幾種方案
??由于時(shí)間序列是模擬的，所以可以控制生成周期、周期+趨勢(shì)、周期+趨勢(shì)+噪聲、周期+噪聲這幾個(gè)類型的時(shí)間序列，根據(jù)對(duì)不同序列的預(yù)測(cè)效果，可以看到ARMA對(duì)不同序列的適用性；
a. 純周期序列
??對(duì)比發(fā)現(xiàn)，周期序列的預(yù)測(cè)效果較好；
b. 純周期+噪聲
??這種情況下，預(yù)測(cè)效果就沒有特別好了，所以在一些情況下可以對(duì)序列進(jìn)行濾波后再進(jìn)行預(yù)測(cè)；
c. 純周期+趨勢(shì)+噪聲
??這種情況預(yù)測(cè)效果不太好；
d. 周期+趨勢(shì)
??沒有噪聲時(shí)這個(gè)預(yù)測(cè)結(jié)果也是可以的。
所以，經(jīng)過對(duì)比能看出，噪聲對(duì)預(yù)測(cè)結(jié)果影響較大，沒有噪聲的情況下，ARMA/ARIMA的預(yù)測(cè)效果大概有一定可信性。文章來源：http://www.zghlxwxcb.cn/news/detail-763804.html
3 Matlab代碼運(yùn)行程序
%% 主程序
%% 1 生成模擬信號(hào)
clear;
n = 132; % 模擬11年的數(shù)據(jù),用10年的數(shù)據(jù)作為訓(xùn)練值
ntest = 12; %用最后1年數(shù)據(jù)作為預(yù)測(cè)比對(duì)值
x = 1:n;
% 周期信號(hào)
a1 = cos(2*pi*x/12) + 0.6 * sin(2*pi*x/24);
figure(1)
subplot(2,2,1)
plot(x,a1);
title('a1 周期信號(hào)')
% 趨勢(shì)信號(hào)
subplot(2,2,2);
x2=0.14;
a2 = mapminmax(x2 * x,0,1);
plot(x,a2);
title('a2 趨勢(shì)信號(hào)')
% 隨機(jī)信號(hào)(0-1的隨機(jī)值)
subplot(2,2,3)
a3 = rand(1,n);
plot(x,a3);
title('a3 隨機(jī)信號(hào)')
% 三種信號(hào)疊加 生成模擬信號(hào)
subplot(2,2,4) %
%a41 = 2 .* a1; % 1 只有周期項(xiàng)
%a41 = 2 .* a1 + 3.2 * a3; % 2 周期 + 隨機(jī)
%a41 = 2 .* a1 + 5 * a2 + 3.2 * a3; % 3 周期 + 趨勢(shì) + 隨機(jī)
a41= 2 .* a1 + 5 * a2;
plot(x,a41);
title('a4 疊加序列')
%% 平穩(wěn)檢驗(yàn)
a4 = a41(1:n-ntest);
ad1 = adftest(a4); % 1平穩(wěn) 0非平穩(wěn)
ad2 = kpsstest(a4); % 0 平穩(wěn) 1 非平穩(wěn)
if (ad1 == 1 && ad2 == 0)
% 平穩(wěn)
type = 3;
else
% 非平穩(wěn)
type = 4;
end
%% type 可以自行輸入
switch type
case 1
datas = 0;
case 2
datas = 0;
case 3 % 平穩(wěn)序列 arma
datas = function_arma(a4,ntest);
%% 繪圖
figure(3)
plot(1:n,a41,'--','LineWidth',1)
hold on
plot((n-ntest)+1:n,datas((n-ntest)+1:n),'--','LineWidth',1.5)
hold on
xlabel('time')
ylabel('value')
legend('真實(shí)值','預(yù)測(cè)值')
case 4 % 非平穩(wěn)序列 arima
[Y,lower,upper] = function_arima(a4,ntest);
figure(3);
plot(1:n,a41,'--','LineWidth',1)
hold on;
plot( (n - ntest+1 : n), Y, 'g', 'LineWidth', 1);
end
%% 平穩(wěn)平穩(wěn) arma預(yù)測(cè)
%% a4訓(xùn)練數(shù)據(jù) ntest預(yù)測(cè)長度
function data = function_arma(a4,ntest)
%自相關(guān)、偏自相關(guān)
figure(2)
subplot(2,1,1)
coll1 = autocorr(a4);
stem(coll1)%繪制經(jīng)線圖
title('自相關(guān)')
subplot(2,1,2)
coll2 = parcorr(a4);
stem(coll2)%繪制經(jīng)線圖
title('偏相關(guān)')
%
count2 = length(a4);
limvalue = round(count2 / 10);
if(limvalue > 10)
limvalue = 10;
end
% 計(jì)算最佳的p 和 q
[bestp,bestq] = Select_Order_arma(a4, limvalue, limvalue);
% 得到最佳模型
xa4 = iddata(a4');
model = armax(xa4,[bestp,bestq]);
% 開始預(yù)測(cè),獲取預(yù)測(cè)數(shù)據(jù)
arr1 = [a4';zeros(ntest,1)];
arr2 = iddata(arr1);
arr3=predict(model, arr2, ntest); %利用arr2去向后預(yù)測(cè)
arr4 = get(arr3);
dataPre = arr4.OutputData{1,1}(length(a4)+1:length(a4)+ntest);
data=[a4';dataPre]; % 獲取訓(xùn)練結(jié)果加預(yù)測(cè)結(jié)果
end
%% Arma找最優(yōu)階數(shù)
function [p1,p2] = Select_Order_arma(x2, limp, limq)
array = zeros(limp, limq);
x = iddata(x2');
for i =1:limp
for j = 1:limq
AmMode = armax(x,[i,j]);
value = aic(AmMode);
array(i,j) = value;
end
end
[p1,p2]=find(array==min(min(array)));
end
%% ARIMA
%% 非平穩(wěn)序列
function [preData,LLL,UUU] = function_arima(a4,ntest)
% 差分
count = 0;
arr1 = a4;
while true
arr21 = diff(arr1);
count = count + 1;
ad1 = adftest(arr21); % 1平穩(wěn) 0非平穩(wěn)
ad2 = kpsstest(arr21); % 0 平穩(wěn) 1 非平穩(wěn)
if ad1==1 && ad2 == 0
break
end
arr1 = arr21;
end
arr22 = arr21';
% 自相關(guān)
figure(2)
subplot(2,1,1)
coll1 = autocorr(arr22);
stem(coll1)%繪制經(jīng)線圖
title('自相關(guān)')
subplot(2,1,2)
coll2 = parcorr(arr22);
stem(coll2)%繪制經(jīng)線圖
title('偏相關(guān)')
% 此時(shí)是平穩(wěn)的
count2 = length(arr22);
limvalue = round(count2 / 10);
if(limvalue > 3)
limvalue = 3;
end
% 計(jì)算最佳的p 和 q
[bestp,bestq] = Select_Order_arima(arr22, limvalue, limvalue, count);
modelo = arima(bestp,count,bestq);
md1 = estimate(modelo, arr22);
[Y, YMSE] = forecast(md1, ntest, 'Y0', arr22);
lower = Y - 1.96*sqrt(YMSE); %95置信區(qū)間下限
upper = Y + 1.96*sqrt(YMSE);
difvalue = [arr22;Y]; %差分序列
difvalueL = [arr22;lower];
difvalueU = [arr22;upper];
count2 = length(a4);
if(count >=1)
hhh = cumsum([a4(1);difvalue]);%還原差分值
LLL1 = cumsum([a4(count2);difvalueL]);
UUU1 = cumsum([a4(count2);difvalueU]);
end
preData = hhh(length(a4)+1 : (length(a4) + ntest));
LLL = LLL1(length(a4)+1 : (length(a4) + ntest));
UUU = UUU1(length(a4)+1 : (length(a4) + ntest));
end
%% Arima找最優(yōu)階數(shù)
function [best_p,best_q] = Select_Order_arima(data, pmax, qmax, d)
min_aic = Inf;
min_bic = Inf;
best_p = 0;
best_q = 0;
for p=0:pmax
for q=0:qmax
model = arima(p,d,q);
try
[fit,~,logL]=estimate(model,data);
[aic, bic] = aicbic(logL, p + q + 1, length(data));
catch
continue
end
if aic < min_aic
min_aic = aic;
min_bic = bic;
best_p = p;
best_q = q;
end
end
end
end
4 總結(jié)
??通過此次練習(xí)，對(duì)ARMA的預(yù)測(cè)過程基本了解，寫Matlab程序的時(shí)候參考了一些博文的邏輯，如有錯(cuò)誤可以進(jìn)行探討。文章來源地址http://www.zghlxwxcb.cn/news/detail-763804.html
到了這里，關(guān)于Arma模型預(yù)測(cè)時(shí)間序列的Matlab實(shí)現(xiàn)的文章就介紹完了。如果您還想了解更多內(nèi)容，請(qǐng)?jiān)谟疑辖撬阉鱐OY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！
