Python和R使用指數加權平均(EWMA)，ARIMA自迴歸移動平均模型預測時間序列

原文連接：http://tecdat.cn/?p=21773

概述

• 學習建立時間序列預測的步驟
• 關注Dickey-Fuller檢驗和ARIMA（自迴歸移動平均）模型
• 從理論上學習這些概念以及它們在python中的實現

介紹

時間序列（從如今起稱爲TS）被認爲是數據科學領域中不爲人知的技能之一。

使用python建立時間序列預測

咱們使用如下步驟：

1. 時間序列是什麼
2. 加載和處理時間序列
3. 如何檢驗時間序列的平穩性？
4. 如何使時間序列平穩？
5. 預測時間序列

1.什麼是時間序列?

顧名思義，TS是以固定時間間隔收集的數據點的集合。

對這些數據進行分析，以肯定長期趨勢，從而預測將來或進行其餘形式的分析。

可是什麼使TS不一樣於常規迴歸問題呢？

1. 它取決於時間。因此線性迴歸模型的基本假設，即觀測值是獨立的，在這種狀況下不成立。
2. 隨着一個增長或減小的趨勢，大多數TS具備某種形式的季節性趨勢，即特定於特定時間範圍的變化。例如，若是你看到一件羊毛夾克的銷售隨着時間的推移，你必定會發現冬季的銷售量更高。

因爲其固有的特性，對其進行分析涉及到各類步驟。下面將詳細討論這些問題。讓咱們從在Python中加載一個TS對象開始。咱們使用航空乘客數據。

2加載和處理時間序列

Pandas有專門的庫來處理TS對象，特別是datatime64[ns]類，它存儲時間信息，咱們能夠快速執行一些操做。讓咱們從啓動所需的庫開始：

%matplotlib inline

from matplotlib.pylab import rcParams 


data = pd.read_csv('AirPassengers.csv', parse_dates=['Month'], index_col='Month')

讓咱們逐一理解這些論點：

1. parse  dates：指定包含日期時間信息的列。如上所述，列名是'Month'。
2. index\_col：將Pandas用於TS數據背後的一個關鍵思想是，索引必須是描述日期時間信息的變量。因此這個參數告訴panda使用'Month'列做爲索引。
3. date parse：指定一個函數，將輸入字符串轉換爲datetime變量。默認狀況下，讀取格式爲「YYYY-MM-DD HH:MM:SS」的數據。若是數據不是此格式，則必須手動定義格式。相似於這裏定義的dataparse函數的東西能夠用於此目的。

如今咱們能夠看到數據以time對象做爲索引，#Passengers做爲列。咱們可使用如下命令交叉檢查索引的數據類型：

注意dtype='datetime[ns]'確認它是一個datetime對象。做爲我的偏好，我會將列轉換爲Series對象，以防止每次使用TS時都引用列名稱。

ts = data[‘#Passengers’] ts.head(10)

在進一步討論以前，我將討論一些TS數據的索引技術。讓咱們從選擇Series對象中的特定值開始。這能夠經過如下兩種方式實現：

#1.指定索引爲字符串常量:
ts['1949-01-01']

#2.導入datetime庫並使用'datetime'函數:
from datetime import datetime

二者都將返回值「112」，這也能夠從之前的輸出中確認。假設咱們想要1949年5月以前的全部數據。這能夠經過兩種方式實現：

#1. 指定整個範圍：
ts['1949-01-01':'1949-05-01']

#2. 若是其中一個索引位於末尾，請使用「：」：
ts[:'1949-05-01']

二者都將產生如下輸出：

這裏有兩點須要注意：

1. 與數字索引不一樣的是，結束索引包含在這裏。例如，若是咱們將列表索引爲[：5]，那麼它將返回索引處的值–[0,1,2,3,4]。但這裏的索引「1949-05-01」包含在輸出中。
2. 必須對索引進行排序，才能使範圍發揮做用。

考慮另外一個須要1949年全部值的例子。這能夠經過如下方式實現：

ts['1949']

月份部分被省略了。相似地，若是您選擇某個月的全部日期，則能夠省略日期部分。
如今，讓咱們繼續分析TS。

3.如何檢驗時間序列的平穩性?

若是TS的統計特性（如均值、方差）隨時間保持不變，則稱其爲平穩的。但爲何重要呢？大多數TS模型都假設TS是平穩的。直覺上，咱們能夠假設，若是一個TS在一段時間內有一個特定的行爲，那麼它頗有可能在未來也會有一樣的行爲。與非平穩序列相比，平穩序列的相關理論更爲成熟和易於實現。
平穩性是使用很是嚴格的標準來定義的。然而，出於實際目的，咱們能夠假設序列是平穩的，若是它隨時間具備恆定的統計特性，即：

1. 恆定平均值
2. 恆定方差
3. 不依賴於時間的自協方差。

我將跳過這些細節，由於在本文中對其進行了很是明確的定義。讓咱們開始測試平穩性的方法。首先也是最重要的是簡單地繪製數據並進行可視化分析。可使用如下命令繪製數據：

plt.plot(ts)

很明顯，隨着一些季節性變化，數據整體呈上升趨勢。然而，可能並不老是可以作出這樣的視覺直觀推斷（咱們稍後會看到這樣的狀況）。所以，更正式地說，咱們可使用如下方法檢查平穩性：

1. 繪製滾動統計：咱們能夠繪製移動平均值或移動方差，看它是否隨時間變化。移動平均值/方差，個人意思是，在任什麼時候刻‘t’，咱們將取去年的平均值/方差，即過去12個月的平均值/方差。但這更像是一種視覺技術。
2. Dickey-Fuller檢驗：這是檢驗平穩性的統計檢驗之一。這裏的零假設是TS是非平穩的。檢驗結果包括檢驗統計量和不一樣置信水平的一些臨界值。若是「檢驗統計量」小於「臨界值」，咱們能夠拒絕零假設並說序列是平穩的。

在這一點上，這些概念聽起來可能不是很直觀。
回到平穩性檢查，咱們將使用滾動統計圖以及Dickey-Fuller檢驗結果不少，因此我定義了一個函數，它以TS做爲輸入併爲咱們生成它們。請注意，我繪製了標準差而不是方差，以保持單位與平均值類似。

def test_stat(timeseries):

    orig = plt.plot(timeseries, color='blue',label='原始序列')
    mean = plt.plot(rolmean, color='red', label='滾動平均')
    std = plt.plot(rolstd, color='black', label = '滾動標準差')
    plt.title('滾動平均數和標準誤差',fontproperties=myfont)


    #執行Dickey-Fumer 檢驗：
    dftest = adfuller(timeseries, maxlag=13, regression='ctt', autolag='BIC')

讓咱們爲輸入序列運行它：

test_stat(ts)

雖然標準差的變化很小，但平均值明顯隨時間增長，這不是一個平穩序列。並且，檢驗統計量遠大於臨界值。請注意，應該比較有符號的值，而不是絕對值。
下一步，咱們將討論能夠用來將這個TS轉換平穩的方法。

4如何使時間序列平穩？

雖然許多TS模型都採用平穩性假設，但實際時間序列中幾乎沒有一個是平穩的。統計學家們已經找到了使序列平穩的方法，咱們如今就來討論。實際上，要使一個序列徹底平穩幾乎是不可能的，但咱們要儘可能接近它。
讓咱們瞭解是什麼使TS不穩定。TS的非平穩性有兩個主要緣由：

1. 趨勢-隨時間變化的平均值。例如，在這個例子中，咱們看到平均而言，乘客人數隨着時間的推移而增加。
2. 季節性-特定時間範圍內的變化。因爲加薪或節日的緣由，人們可能有在特定月份出現的傾向。

其基本原理是對序列中的趨勢和季節性進行建模或估計，並從序列中去除這些趨勢和季節性以獲得平穩序列。而後統計預測技術就能夠在這個序列上實現。最後一步是經過應用趨勢和季節性約束將預測值轉換爲原始規模。
注意：我將討論一些方法。在這種狀況下，有些可能工做得很好，而有些可能不行。但咱們的想法是掌握全部的方法，而不是隻關注手頭的問題。
讓咱們從趨勢部分開始。

估計和消除趨勢

減小趨勢的第一個技巧之一能夠轉換。例如，在這種狀況下，咱們能夠清楚地看到存在顯着的增加趨勢。所以，咱們能夠應用轉換，以懲罰更高的值，而不是較小的值。這些能夠取對數，平方根，立方根等。在此處取對數轉換簡便簡單：

np.log(ts)

在這種簡單的狀況下，很容易看到數據的將來趨勢。但在有噪音的狀況下，它不是很直觀。所以，咱們可使用一些技術來估計或模擬這一趨勢，而後將其從序列中刪除。有不少方法能夠作到這一點，其中最經常使用的有：

1. 聚合-取一段時間的平均值，如月/周平均值
2. 平滑-取滾動平均值
3. 多項式擬合-擬合迴歸模型

我將在這裏討論平滑，你應該嘗試其餘技術，以及可能解決其餘問題。平滑是指採起滾動估計，即考慮過去的幾個實例。有不少種方法，但我將在這裏討論其中的兩種。

移動（滾動）平均

在這種方法中，咱們根據時間序列的頻率取「k」連續值的平均值。這裏咱們能夠取過去一年的平均值，也就是最近12個值。pandas有肯定滾動統計的特定功能。

pd.rolling_mean(ts_log,12)

紅線表示滾動平均值。咱們從原來的數列中減去這個。注意，由於咱們取最後12個值的平均值，因此滾動平均值沒有定義前11個值。這能夠觀察到:

ts_log - moving_avg

注意前11個是Nan。讓咱們刪除這些NaN值並檢查圖以檢驗平穩性。

moving_avg_diff.dropna(inplace=True)

這看起來是一個更好的系列。滾動值彷佛略有變化，但沒有具體趨勢。並且，檢驗統計量小於5%的臨界值，因此咱們能夠用95%的置信度說這是一個平穩序列。

加權移動平均值（EWMA）

然而，這種方法的一個缺點是，時間段必須嚴格定義，在這種狀況下，咱們能夠取年平均數，但在預測股票價格等複雜狀況下，很可貴出一個數字。因此咱們取一個「加權移動平均值」，其中最近的值被賦予更高的權重。分配權重的方法有不少種，一種流行的方法是指數加權移動平均法，即用衰減因子將權重分配給全部先前的值。請在此處查找詳細信息。這能夠實現爲：

ewma(ts_log, half=12)

請注意，這裏的參數「半衰期」用於定義指數衰減的量。這只是一個假設，很大程度上取決於業務領域。其餘參數，如跨度和質心，也能夠用來定義衰變，這將在上面共享的連接中討論。如今，讓咱們從序列中刪除此項並檢查平穩性：

ts_log - expwighted_avg

這個TS在平均值和標準誤差上的變化甚至更小。此外，檢驗統計量小於1%的臨界值，這比前一種狀況要好。注意，在這種狀況下，不會出現缺失值，由於從開始算起的全部值都是給定的權重。因此即便沒有之前的值，它也能工做。

消除趨勢和季節性

之前討論過的簡單的趨勢減小技術並不適用於全部狀況，特別是那些具備高季節性的狀況。讓咱們討論兩種消除趨勢和季節性的方法：

1. 差分-用特定的時間差取差分
2. 分解–對趨勢和季節性進行建模，並將其從模型中移除。

差分

處理趨勢性和季節性最經常使用的方法之一是差分。在這項技術中，咱們取某一時刻的觀測值與前一時刻的觀測值之差。這在改善平穩性方面效果很好。一階差分可按以下方式進行：

ts_log - ts_log.shift()

這彷佛大大減小了這一趨勢。讓咱們使用繪圖進行驗證：

咱們能夠看到平均值和標準差隨時間的變化很小。同時，Dickey-Fuller檢驗統計量小於10%的臨界值，所以TS是平穩的，置信度爲90%。咱們也能夠取二階或三階差，在某些應用中可能會獲得更好的結果。

分解

在這種方法中，趨勢和季節性都被分別建模，序列的其他部分被返回。我將跳過統計數據得出結果：

decompose(ts_log)

 decomposition.trend
 decomposition.seasonal
 decomposition.resid

在這裏咱們能夠看到趨勢，季節性是從數據中分離出來的，咱們能夠對殘差進行建模。讓咱們檢查殘差的平穩性：

ts_log_decompose = residual

Dickey-Fuller檢驗統計量顯著低於1%的臨界值。因此這個t很是接近平穩。另外，您應該注意到，在本例中，將殘差轉換爲將來數據的原始值不是很直觀。

5預測時間序列

咱們看到了不一樣的技術，全部這些技術都至關好地使TS平穩。由於這是一種很是流行的技術，因此讓咱們在差分後在TS上創建模型。此外，在這種狀況下，將噪聲和季節性添加回預測殘差中相對容易。執行趨勢和季節性估計技術後，可能有兩種狀況：

1. 一個嚴格平穩的序列，值之間沒有依賴關係。這是一個簡單的例子，咱們能夠將殘差建模爲白噪聲。但這是很是罕見的。
2. 值之間有顯著相關性的一系列。在這種狀況下，咱們須要使用一些統計模型，如ARIMA來預測數據。

讓我簡單介紹一下ARIMA。我不會深刻討論技術細節，但若是您但願更有效地應用它們，您應該詳細瞭解這些概念。ARIMA表明自迴歸綜合移動平均線。平穩時間序列的ARIMA預測只不過是一個線性(相似於線性迴歸)方程。預測因子取決於ARIMA模型的參數(p,d,q):

1. AR(自迴歸)項數(p): AR項是因變量的滯後。例如，若是p爲5,x(t)的預測因子爲x(t-1)....x(t-5)。
2. 移動平均線項數(q):移動平均線項是預測方程中的滯後預測偏差。例如，若是q是5,x(t)的預測因子將是e(t-1)....e(t-5)，其中e(i)是移動平均在第一個瞬間和實際值之間的差值。
3. 差分的數量(d):這些是非季節性差分的數量，即在這種狀況下，咱們取一階差分。因此咱們能夠傳遞這個變量，讓d=0或者傳遞原始變量，讓d=1。二者都會產生相同的結果。

這裏的一個重要問題是如何肯定「p」和「q」的值。咱們先來討論一下。

1. 自相關函數（ACF）：它是TS與自身滯後之間相關性的度量。例如，在滯後5時，ACF會將「t1」…「t2」時刻的序列與「t1-5」…「t2-5」時刻的序列進行比較（t1-5和t2是終點）。
2. 偏自相關函數（PACF）：它測量了TS之間的相關性， 但在消除了已經解釋的變化以後。例如，在滯後5，它將檢查相關性，但消除已經由滯後1到4解釋的影響。

差分後TS的ACF和PACF圖能夠繪製爲：

#繪圖ACF：
plt.subplot(121) 
plt.plot(lag_acf)
plt.axhline(y=0,linestyle='--',color='gray')

#繪圖PACF：
plt.subplot(122)
plt.plot(lag_pacf)
plt.axhline(y=0,linestyle='--',color='gray')

plt.tight_layout()

在這個圖中，0兩邊的兩條虛線是置信區間。這些可用於肯定「p」和「q」值，以下所示：

1. p–PACF圖表第一次穿過置信區間上限的滯後值。若是你仔細觀察，在這種狀況下，p=2。
2. q–ACF圖表首次穿過置信區間上限的滯後值。若是你仔細觀察，在這種狀況下，q=2。

如今，讓咱們創建3種不一樣的ARIMA模型，既考慮個體效應，也考慮組合效應。我也會輸出RSS。請注意，這裏的RSS是殘差值，而不是實際序列。
咱們須要先加載ARIMA模型：

可使用ARIMA的order參數指定p、d、q值，該參數採用元組(p、d、q)。讓咱們模擬3個案例:

AR模型

model.fit(disp=-1)
plt.plot(ts_log_diff)

MA模型

results_MA = model.fit(disp=-1)
plt.plot(ts_log_diff)

組合模型

plt.plot(ts_log_diff)
plt.plot(results_ARIMA.fittedvalues, color='red')

在這裏，咱們能夠看到AR和MA模型有幾乎相同的RSS，可是結合起來會更好。如今，咱們剩下最後一步，即將這些值恢復到原始比例。

回到原來的比例

因爲組合模型給出了最佳結果，讓咱們將其縮放回原始值，看看它的預測表現如何。第一步是將預測結果存儲爲一個單獨的序列並觀察它。

請注意，這些開始於'1949-02-01'，而不是第一個月。爲何？這是由於咱們取了一個滯後1，第一個元素以前沒有任何東西能夠減去。將差分轉換爲對數刻度的方法是將這些差分連續地添加到基數上。一個簡單的方法是首先肯定索引處的累積和，而後將其添加到基數中。累積總和以下所示：

您可使用前面的輸出快速地作一些計算，以檢查這些是否正確。接下來咱們要把它們加到底數上。爲此，讓咱們建立一個以全部值做爲基數的序列，並將其差分相加。這能夠這樣作:

這裏的第一個元素是基數自己，而後從那裏累計加值。最後一步是取指數並與原數列進行比較。

plt.plot(predictions_ARIMA)

這些都是Python中的內容。咱們來學習一下如何在R中實現時間序列預測。

R時間序列預測

第一步：讀取數據，計算基本總結

#安裝包並調用庫
install.packages("tseries")

#讀取Airpaseengers數據
tsdata<-AirPassengers
#識別數據類別
class(tsdata)
#觀察時間序列數據
tsdata
#數據摘要
dfSummary(tsdata)

輸出

class(tsdata)
"ts"
> #觀察時間序列數據
> tsdata
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
1949 112 118 132 129 121 135 148 148 136 119 104 118
1950 115 126 141 135 125 149 170 170 158 133 114 140
1951 145 150 178 163 172 178 199 199 184 162 146 166
1952 171 180 193 181 183 218 230 242 209 191 172 194
1953 196 196 236 235 229 243 264 272 237 211 180 201
1954 204 188 235 227 234 264 302 293 259 229 203 229
1955 242 233 267 269 270 315 364 347 312 274 237 278
1956 284 277 317 313 318 374 413 405 355 306 271 306
1957 315 301 356 348 355 422 465 467 404 347 305 336
1958 340 318 362 348 363 435 491 505 404 359 310 337
1959 360 342 406 396 420 472 548 559 463 407 362 405
1960 417 391 419 461 472 535 622 606 508 461 390 432
> #

tsdata 
Dimensions: 144 x 1 
Duplicates: 26

----------------------------------------------------------------------------------------------------
No Variable Stats / Values Freqs (% of Valid) Graph Valid Missing 
---- ---------- -------------------------- ----------------------- --------------------- -------- --
1 tsdata Mean (sd) : 280.3 (120) 118 distinct values . : . 144 0 
[ts] min < med < max: Start: 1949-01 : : . . : (100%) (0%) 
104 < 265.5 < 622 End : 1960-12 : : : : : 
IQR (CV) : 180.5 (0.4) : : : : : : : 
: : : : : : : : . . 
----------------------------------------------------------------------------------------------------

第2步：檢查時間序列數據的週期並繪製原始數據

#檢查數據並繪製原始數據

plot(tsdata, ylab="乘客（1000人）", type="o")

輸出

cycle(tsdata)
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
1949 1 2 3 4 5 6 7 8 9 10 11 12
1950 1 2 3 4 5 6 7 8 9 10 11 12
1951 1 2 3 4 5 6 7 8 9 10 11 12
1952 1 2 3 4 5 6 7 8 9 10 11 12
1953 1 2 3 4 5 6 7 8 9 10 11 12
1954 1 2 3 4 5 6 7 8 9 10 11 12
1955 1 2 3 4 5 6 7 8 9 10 11 12
1956 1 2 3 4 5 6 7 8 9 10 11 12
1957 1 2 3 4 5 6 7 8 9 10 11 12
1958 1 2 3 4 5 6 7 8 9 10 11 12
1959 1 2 3 4 5 6 7 8 9 10 11 12
1960 1 2 3 4 5 6 7 8 9 10 11 12

步驟3:分解時間序列數據

#將數據分解爲趨勢、季節性和隨機偏差份量
plot(tsdata_decom)

輸出

第四步：檢驗數據的平穩性

#測試數據的平穩性
#單位根檢驗
adf(tsdata)

#自相關檢驗
plot(acf(tsdata,plot=FALSE))+ labs(title="航空旅客數據相關圖")

plot(acf(tsdata_decom$random,plot=FALSE))

輸出

Augmented Dickey-Fuller Test

data: tsdata
Dickey-Fuller = -7.3186, Lag order = 5, p-value = 0.01
alternative hypothesis: stationary

the p-value is 0.01 which ip值爲0.01，小於0.05，所以，咱們拒絕了零假設，所以時間序列是平穩的。s <0.05, therefore, we reject the null hypothesis and hence time series is stationary.

最大滯後爲1個月或12個月，代表與12個月週期正相關。
自動繪製7:138的隨機時間序列觀測值，不包括NA值

步驟5：擬合模型

#擬合模型
#線性模型
plot(tsdata) + smooth(method="lm")


#ARIMA 模型
arimats

Series: tsdata 
ARIMA(2,1,1)(0,1,0)[12]

Coefficients:
ar1 ar2 ma1
0.5960 0.2143 -0.9819
s.e. 0.0888 0.0880 0.0292

sigma^2 estimated as 132.3: log likelihood=-504.92
AIC=1017.85 AICc=1018.17 BIC=1029.35

第6步：預測

#Arima模型的預測
    fore(arimats, level = c(95))

最後咱們有一個原始數據的預測。

• • -

最受歡迎的看法

1.在python中使用lstm和pytorch進行時間序列預測

2.python中利用長短時間記憶模型lstm進行時間序列預測分析

3.使用r語言進行時間序列（arima，指數平滑）分析

4.r語言多元copula-garch-模型時間序列預測

5.r語言copulas和金融時間序列案例

6.使用r語言隨機波動模型sv處理時間序列中的隨機波動

7.r語言時間序列tar閾值自迴歸模型

8.r語言k-shape時間序列聚類方法對股票價格時間序列聚類

9.python3用arima模型進行時間序列預測