你有沒有遇到過這樣的問題：我有一段資料，它是随着時間等間隔采樣的，現在想用某種方法預測出後續一段時間的趨勢。這就是所謂的時間序列的預測問題。
時間序列預測的應用主要是在經濟領域，比如預測股市行情，預測GDP走勢；也可以用來預測機械性能退化趨勢或者其他諸如某個量随時間變化這樣的場景當中。
AR/MA/ARMA模型是分析時間序列的重要方法，在本篇文章中将重點介紹如何使用ARMA(p,q)模型對時間序列信号進行一次完整的預測。而對于ARMA(p,q)模型的相關理論知識則不進行重點介紹。沒有涉及到的理論部分，建議閱讀相關專業書籍。
本文的目的是，對于不了解時間序列與ARMA(p,q)的同學，可以套用文中的方法得到想要的預測結果；對于了解理論知識的同學，通過代碼加深了解。
是以，這不是一篇ARMA模型入門的文章，而是一篇應用ARMA進行預測的傻瓜式示範教程，就算你對裡邊的理論完全不懂，相信也照樣可以依樣畫葫蘆做出預測（雖然不推薦這樣！）。
1.導入實驗資料
close all
clear all
load Data_EquityIdx %納斯達克綜合指數
len = 120;
Y = DataTable.NASDAQ(1:len);
plot(Y)
2.平穩性檢驗
使用ARMA模型要求時間序列必須是平穩的，是以第一步是對原始資料進行平穩性檢驗。檢驗方法有很多種，包括ADF、KPSS、P-P等。這裡用ADF檢驗和KPSS檢驗。
y_h_adf = adftest(Y)
y_h_kpss = kpsstest(Y)
得到y_h_adf =0，y_h_kpss =1，沒有通過檢驗。
此時對Y取差分，再進行檢驗：
% 一階差分，結果平穩。如果依舊不平穩的話，再次求差分，直至通過檢驗
Yd1 = diff(Y);
yd1_h_adf = adftest(Yd1)
yd1_h_kpss = kpsstest(Yd1)
得到yd1_h_adf =1，yd1_h_kpss =0，通過檢驗了。
如果此時并沒有通過檢驗的話，可以取二階差分，再不行就三階，直到信号變得平穩，此時差分的階數就是ARIMA(p,i,q)中的i的值。上述例子中一階差分就平穩了，是以i=1。
3.确定ARMA模型階數
模型階數的确定方法主要包括兩種：（1）ACF和PACF法。（2）基于準則的确定法。（都是我自己起的名）
（1）ACF和PACF法
所謂ACF即自相關函數，PACF為偏自相關函數。
% ACF和PACF法，确定階數
figure
autocorr(aimY)
figure
parcorr(aimY)
做出的結果如下兩圖：
自相關函數ACF
偏自相關函數PACF
怎樣根據這兩站圖确定p和q呢，有這樣一張圖：
通過ACF和PACF一般特征
對于我們關注的ARMA(p,q)，通俗地說，PACF最後一個在藍線外（即門檻值外）的Lag值就是p值；ACF最後一個在藍線外（即門檻值外）的Lag值就是q值。
如果按照這個标準，p似乎要定到13，q也要定到12，這個值有些偏大了。
（2）基于準則的确定法
如果方法（1）不理想，可以用方法（2），這個方法不用糾結，可以直接給出推薦結果。代碼如下：
% 通過AIC，BIC等準則暴力標明階數
max_ar = 3;
max_ma = 3;
[AR_Order,MA_Order] = ARMA_Order_Select(Y,max_ar,max_ma,1);
ARMA_Order_Select是我自己寫的一個函數，隻需要輸入信号、最大p值、最大q值和差分階數就可以得到推薦p和q。
運作結果為：
AR_Order = 3
MA_Order = 3
即p=1,q=3。
4.殘差檢驗
為了確定确定的階數合适，還需要進行殘差檢驗。殘差即原始信号減掉模型拟合出的信号後的殘餘信号。如果殘差是随機正态分布的、不自相關的，這說明殘差是一段白噪聲信号，也就說明有用的信号已經都被提取到ARMA模型中了。
Mdl = arima(AR_Order, 1, MA_Order); %第二個變量值為1，即一階差分
EstMdl = estimate(Mdl,data);
[res,~,logL] = infer(EstMdl,data); %res即殘差
stdr = res/sqrt(EstMdl.Variance);
figure('Name','殘差檢驗')
subplot(2,3,1)
plot(stdr)
title('Standardized Residuals')
subplot(2,3,2)
histogram(stdr,10)
title('Standardized Residuals')
subplot(2,3,3)
autocorr(stdr)
subplot(2,3,4)
parcorr(stdr)
subplot(2,3,5)
qqplot(stdr)
殘差檢驗
上圖為殘差檢驗的結果圖。Standardized Residuals是檢視殘差是否接近正态分布，理想的殘差要接近正态分布；ACF和PACF檢驗殘差的自相關和偏自相關，理想的結果應該在圖中不存在超出藍線的點；最後一張QQ圖是檢驗殘差是否接近正太分布的，理想的結果中藍點應該靠近紅線。
除了上述圖像檢驗方法，還可以通過Durbin-Watson對相關性進行檢驗：
% Durbin-Watson 統計是計量經濟學分析中最常用的自相關度量
diffRes0 = diff(res);
SSE0 = res'*res;
DW0 = (diffRes0'*diffRes0)/SSE0 % Durbin-Watson statistic，該值接近2，則可以認為序列不存在一階相關性。
運算結果為1.9997，這個值越接近2越說明殘差不存在一階相關性。
上述檢驗可以證明，殘差接近正态分布，且互相獨立，可以認為ARMA模組化符合要求。
5.預測
萬事俱備，終于來到最後一步。
%% 5.預測
step = 300; %預測步數為300
[forData,YMSE] = forecast(EstMdl,step,'Y0',data); %matlab2018及以下版本寫為Predict_Y(i+1) = forecast(EstMdl,1,'Y0',Y(1:i));
lower = forData - 1.96*sqrt(YMSE); %95置信區間下限
upper = forData + 1.96*sqrt(YMSE); %95置信區間上限
figure()
plot(data,'Color',[.7,.7,.7]);
hold on
h1 = plot(length(data):length(data)+step,[data(end);lower],'r:','LineWidth',2);
plot(length(data):length(data)+step,[data(end);upper],'r:','LineWidth',2)
h2 = plot(length(data):length(data)+step,[data(end);forData],'k','LineWidth',2);
legend([h1 h2],'95% 置信區間','預測值',...
'Location','NorthWest')
title('Forecast')
hold off
多步預測結果
上圖中灰線為用來訓練的1200個資料，黑線為未來值的預測，紅線為95%置信區間上下限。也就是說未來真實值有95%的機率落在這個範圍内。
可以看到使用ARIMA方法進行長期預測的結果是趨勢性的。
6.封裝
上述整個過程可以封裝成一個函數，如下：
%% 進行使用ARIMA進行預測的函數
function [forData,lower,upper] = Fun_ARIMA_Forecast(data,step,max_ar,max_ma,figflag)
% 使用ARIMA進行預測的函數（使用n階差分、不使用對數），可以直接調用，差分階數自動确定。
% 輸入：
% data為待預測資料，一維資料，最小11個資料。但是資料長度處于11~15時依舊可能出現報錯的情況。
% step為拟預測步數
% max_ar 為最大p值
% max_ma 為最大q值
% figflag 為畫圖示志位，'on'為畫圖，'off'為不畫
% 輸出：
% forData為預測結果，其長度等于step
% lower為預測結果的95%置信下限值
% upper為預測結果的95%置信上限值
直接調用就可以出結果。
使用上邊封裝好的函數試一下循環疊代的單步預測效果：其思路就是比如從30開始，用前30個資料預測第31個；再用前31個預測第32個...用前100個預測第101個。然後再把所有單步預測結果畫出來和原資料做對比。
%% 使用Fun_ARIMA_Forecast函數實作預測
% 單步預測
close all
clear all
addpath ../funs %将funs檔案夾添加進路徑
load Data_EquityIdx %納斯達克綜合指數
len = 100;
data = DataTable.NASDAQ(1:len); %如果要替換資料，将此處data替換即可。
forData1 = zeros(1,len); %全部初始化為0
for i = 30:len
forData1(i+1) = Fun_ARIMA_Forecast(data(1:i),1,2,2,'off');
end
figure()
plot(31:len,data(31:len))
hold on
plot(31:len,forData1(31:len))
legend('原始資料','單步預測資料')
畫出的結果如下（對比從30~100的資料段）：
如果要擷取本文的完整代碼，可以關注我的公衆号“括号的城堡”，微信号為“khscience”，回複“時間序列”，公衆号裡可能還會有更多有趣的東西分享。
為了不再開坑，一篇講完了
參考：
MATLAB 文檔中心 - MathWorks 中國
【時間序列分析】ARMA預測GDP的python實作