Python面板时间序列数据预测：格兰杰因果关系检验Granger causality test药品销售实例与可视化
最近我們被客戶要求撰寫關于時間序列的研究報告，包括一些圖形和統計輸出。
時間序列是以固定時間區間記錄的觀察序列。本指南帶你完成在Python中分析一個給定的時間序列的特征的過程。
內容
1. 什么是時間序列？
時間序列是以固定時間區間記錄的觀察序列。
根據觀察的頻率，一個時間序列通常可能是每小時、每天、每周、每月、每季度和每年。有時，你也可能有以秒為單位的時間序列，比如，每分鐘的點擊量和用戶訪問量等等。
為什么要分析一個時間序列？
因為這是你對該序列進行預測前的準備步驟。
此外，時間序列預測具有巨大的商業意義，因為對企業來說很重要的東西，如需求和銷售，網站的訪問量，股票價格等基本上都是時間序列數據。
那么，分析一個時間序列涉及什么呢？
時間序列分析涉及到對序列性質的各個方面的理解，這樣你就能更好地了解創造有意義和準確的預測。
2. 如何在Python中導入時間序列？
那么，如何導入時間序列數據呢？
時間序列的數據通常存儲在.csv文件或其他電子表格格式中，包含兩列：日期和測量值。
我們使用pandas包中的read_csv()來讀取時間序列數據集（一個關于藥品銷售的csv文件）作為pandas數據框。添加parse_dates=['date']參數將使日期列被解析為一個日期字段。
import pandas as pd# 導入數據 df = pd.read_csv('10.csv ', parse_dates=['date']) df.head()數據框架時間序列
另外，你也可以把它導入為一個以日期為索引的pandas序列。你只需要在pd.read_csv()中指定index_col參數就可以了。
pd.read_csv('10.csv', parse_dates=['date'], index_col='date')3. 什么是面板數據？
面板數據也是一種基于時間的數據集。
不同的是，除了時間序列之外，它還包含一個或多個在相同時間段內測量的相關變量。
通常情況下，面板數據中存在的列包含了有助于預測Y的解釋變量，前提是這些列在未來的預測期是可用的。
下面是一個面板數據的例子。
df.head()面板序列
4. 時間序列的可視化
讓我們用matplotlib來可視化這個序列。
# 時間序列數據源：R中的fpp pacakge。 import matplotlib.pyplot as plt# 繪制圖表 def plot_df(df, x, y, title="", xlabel='日期', dpi=100):plt.show()時間序列的可視化
由于所有的值都是正數，你可以在Y軸的兩邊顯示，以強調增長。
# 導入數據 x = df['date'].values# 繪圖 fig, ax = plt.subplots(1, 1, figsize=（16,5）, dpi= 120) plt.fill(x, y1=y1, y2=-y1, alpha=0.5)航空乘客數據--兩面序列
由于它是一個月度的時間序列，并且每年都遵循一定的重復模式，你可以在同一張圖中把每年的情況作為一個單獨的線條來繪制。這讓你可以并排比較每年的模式。
時間序列的季節圖
# 導入數據df.reset_index(inplace=True)# 準備好數據years = df['year'].unique()# 預備顏色np.random.choice(list(mpl.color), len(year), # 繪制圖表plt.text(df.loc[df.year==y, :].shape[0]-.9] plt.gca().set(xlim=(-0.3, 11)plt.title("藥品銷售時間序列的季節圖", fontsize=20)藥品銷售的季節性圖譜
每年2月，藥品銷售量急劇下降，3月再次上升，4月再次下降，如此反復。顯然，這種模式每年都會在某一年內重復出現。
然而，隨著時間的推移，藥品銷售量總體上有所增加。你可以用一個漂亮的年度圖表很好地展示這一趨勢以及它每年的變化情況。同樣地，你也可以做一個按月排列的boxplot來顯示每月的分布情況。
逐月（季節性）和逐年（趨勢）分布的箱線圖
你可以按季節性對數據進行分組，看看數值在某年或某月是如何分布的，以及它在不同時期的對比情況。
# 導入數據df.reset_index(inplace=True)# 準備好數據 df['年'] = [d.year for d in df.date] df['月'] = [d.strftime('%b') for d in df.date]# 繪制圖表 sns.boxplot(x='年', y='值', data=df, ax=axes[0]) sns.boxplot(x='月', y='值', data=df.loc[~df.year.isin([1991, 2008]), :] )按年和按月排列的箱線圖
箱線圖使年度和月份的分布變得明顯。另外，在按月排列的圖表中，12月和1月的藥品銷售量明顯較高，這可歸因于假日折扣季節。
到目前為止，我們已經看到了識別模式的相似性。現在，如何找出與通常模式的任何偏差？
5. 時間序列中的模式
任何時間序列都可以被分成以下幾個部分。?基礎水平+趨勢+季節性+誤差
當在時間序列中觀察到有一個增加或減少的斜率時，就可以觀察到趨勢。而季節性是指由于季節性因素，在定期區間之間觀察到明顯的重復模式。這可能是由于一年中的哪個月，哪個月的哪一天，工作日或甚至一天中的哪個時間。
然而，并非所有時間序列都必須有趨勢和/或季節性。一個時間序列可能沒有一個明顯的趨勢，但有一個季節性。反之，也可以是真實的。
因此，一個時間序列可以被想象為趨勢、季節性和誤差項的組合。
? fig, axes = plt.subplots(1,3, figsize=(20,4), dpi=100)pd.read_csv.plot( legend=False, ax=axes[2])?時間序列中的模式
另一個需要考慮的方面是周期性行為。當序列中的上升和下降模式不發生在固定的基于日歷的時間區間內時，就會發生這種情況。應注意不要將 "周期性 "效應與 "季節性 "效應混淆。
那么，如何區分 "周期性 "和 "季節性 "模式？
如果這些模式不是基于固定的日歷頻率，那么它就是周期性的。因為，與季節性不同，周期性效應通常受到商業和其他社會經濟因素的影響。
6. 加法和乘法的時間序列
根據趨勢和季節性的性質，一個時間序列可以被建模為加法或乘法，其中，序列中的每個觀測值可以表示為各組成部分的和或積。
加法時間序列:值=基礎+趨勢+季節性+誤差
?
乘法時間序列:值=基礎x趨勢x季節性x誤差
7. 如何將一個時間序列分解成其組成部分？
你可以對一個時間序列進行經典的分解，將該序列視為基數、趨勢、季節性指數和殘差的加法或乘法組合。
statsmodels中的 seasonal_decompose 可以方便地實現這一點。
# 乘法分解 decompose(df['value'], model='multiplicative')# 加法分解 decompose(df['value'], model='additive')# 繪圖 result_mul.plot().suptitle( fontsize=22)加法和乘法分解
設置extrapolate_trend='freq'可以看到序列開始時趨勢和殘差中的任何缺失值。
如果你仔細看一下加法分解的殘差，它有一些模式殘留。然而，乘法分解看起來相當隨機。因此，理想情況下，對于這個特定的序列，乘法分解應該是首選。
趨勢、季節性和殘差成分的數字輸出存儲在result_mul輸出本身。讓我們把它們提取出來，放在一個數據框中。
# 提取成分---- # 實際值=（季節性*趨勢*殘差）的乘積constructed.columns = ['seas', 'trendance', 'resid', 'actual_values'] 。 constructed.head()如果你檢查一下，seas、trend和resid列的乘積應該正好等于actual_values。
8. 平穩和非平穩的時間序列
平穩性是時間序列的一個屬性。一個平穩的序列是指該序列的值不是時間的函數。
也就是說，序列的統計屬性，如平均數、方差和自相關，隨著時間的推移是不變的。序列的自相關只不過是序列與它以前的值的相關性，更多的是關于這一點。
一個平穩的時間序列也沒有季節性影響。
那么，如何識別一個序列是否是平穩的？讓我們繪制一些例子來說明。
平穩和非平穩的時間序列
那么，為什么平穩的序列很重要呢？
我稍后會說到這個問題，但要知道，通過應用適當的轉換，幾乎可以使任何時間序列成為平穩的。大多數統計預測方法都被設計為在平穩的時間序列上工作。預測過程的第一步通常是做一些轉換，把非平穩序列轉換成平穩的。
9. 如何使一個時間序列平穩？
你可以通過以下方式使序列平穩
使序列平穩的最常見和最方便的方法是對數列至少進行一次差分，直到它成為近似平穩的。
那么，什么是差分？
如果Y_t是時間't'的值，那么Y的第一個差值=Yt-Yt-1。更簡單地說，序列只不過是用當前值減去下一個值。
如果第一次差分不能使一個序列平穩，你可以進行第二次差分。以此類推。
例如，考慮以下序列。[1, 5, 2, 12, 20]
第一次差分可以得到。[5-1, 2-5, 12-2, 20-12] = [4, -3, 10, 8]
二次差分得出。[-3-4, -10-3, 8-10] = [-7, -13, -2]
9. 為什么在預測前要使非穩態序列成為穩態？
預測一個平穩的序列是相對容易的，而且預測結果也更可靠。
一個重要的原因是，自回歸預測模型本質上是線性回歸模型，利用序列本身的滯后期作為預測因子。
我們知道，如果預測因子（X變量）不相互關聯，線性回歸效果最好。因此，序列的平穩化解決了這個問題，因為它消除了任何持續的自相關，從而使預測模型中的預測因子（序列的滯后）幾乎是獨立的。
現在我們已經確定了序列平穩化的重要性，那么你如何檢查一個給定的序列是否是平穩的？
10. 如何檢驗平穩性？
一個序列的平穩性可以通過觀察序列的圖表來確定，就像我們之前做的那樣。
另一種方法是將序列分成2個或更多的連續部分，并計算平均數、方差和自相關等匯總統計數據。如果統計數字有很大差異，那么這個序列就不可能是平穩的。
然而，你需要一種方法來定量地確定一個給定的序列是否是平穩的。這可以通過稱為 "單位根測試 "的統計測試來完成。這方面有多種變化，測試檢查一個時間序列是否是非穩態的并擁有單位根。
單位根測試有多種實現方式，例如。
最常用的是ADF檢驗，無效假設是時間序列擁有單位根，并且是非平穩的。所以，如果ADH檢驗中的P值小于顯著性水平（0.05），你就拒絕無效假設。
另一方面，KPSS檢驗是用來檢驗趨勢平穩性的。空假設和P值的解釋與ADH檢驗正好相反。下面的代碼使用python中的statsmodels包來實現這兩個檢驗。
# ADF測試 result = adfuller(df.value.values, autolag='AIC')# KPSS測試 result = kpss(df.value.values, regression='c') print('\nKPSS Statistic: %f' % result[0])?
11. 白噪聲和平穩數列之間有什么區別？
與平穩序列一樣，白噪聲也不是時間的函數，即它的平均值和方差不隨時間變化。但不同的是，白噪聲是完全隨機的，平均值為0。
在白噪聲中，沒有任何模式。如果你把調頻收音機中的聲音信號看作是一個時間序列，你在各頻道之間聽到的空白聲音就是白噪聲。
在數學上，一個平均數為0的完全隨機的數字序列就是白噪聲。
np.random.randn(1000)隨機白噪聲
12. 如何對時間序列進行去趨勢處理？
去時間序列的趨勢是指從一個時間序列中去除趨勢成分。但如何提取趨勢呢？有多種方法。
去除我們前面看到的從時間序列分解中得到的趨勢成分。
去除平均數
應用像Baxter-King濾波器或Hodrick-Prescott濾波器這樣的濾波器來去除移動平均趨勢線或周期成分。
讓我們來實現前兩種方法。
# 使用scipy。減去最佳擬合線signal.detrend(df.value.values) plt.plot(detrended)通過減去最小二乘法來解讀時間序列的趨勢
# 使用statmodels。減去趨勢成分。 decompose(df['value'], model='multiplicative', detrended = df.value.value - trend plt.plot(detrended)通過減去趨勢成分進行去勢
13. 如何對一個時間序列進行去季節化？
有多種方法可以使時間序列去季候性化。下面是幾個例子。
- 1.取一個以季節性窗口為長度的移動平均線。這將在這個過程中使序列變得平滑。- 2.序列的季節性差異（用當前值減去前一季的值）。- 3.用從STL分解得到的季節性指數除以該序列。如果除以季節性指數效果不好，可以嘗試取序列的對數，然后進行去季節性處理。之后你可以通過取指數來恢復到原來的規模。
# 時間序列分解 seasonal_decompose(df['value'], model='multiplicative'')# 去季節化 df. value.values / result_mul.seasonal # 繪圖 plt.plot(deseason)對時間序列進行反季節處理
14. 如何測試一個時間序列的季節性？
常見的方法是繪制序列圖，檢查平穩時間區間內的可重復模式。所以，季節性的類型是由時鐘或日歷決定的。
然而，如果你想對季節性有一個更明確的檢查，可以使用自相關函數（ACF）圖。更多關于ACF的內容將在接下來的章節中介紹。但是，當有強烈的季節性模式時，ACF圖通常會顯示出在季節性窗口的倍數上有明確的重復峰值。
例如，藥品銷售時間序列是一個月度序列，每年都有重復的模式。因此，你可以看到在第12、24、36......行的尖峰。
在真實的數據集中，這種強烈的模式很難被注意到，并可能被任何噪音所扭曲，所以你需要仔細觀察。
?# 畫圖 correlation_plot(df.value.tolist())?自相關圖
另外，如果你想進行統計測試，CHTest可以確定是否需要進行季節性差分以使序列平穩化。
15. 如何處理時間序列中的缺失值？
有時，你的時間序列會有缺失的日期/時間。這意味著，這些時期的數據沒有被捕獲或無法獲得。在這種情況下，你可以用零來填補這些時期。
其次，當涉及到時間序列時，你通常不應該用序列的平均值來替換缺失值，特別是在序列不是平穩的情況下。你可以做的是是向前填充前一個值。
然而，根據序列的性質，你要在得出結論之前嘗試多種方法。一些有效的替代歸因法的方法是。
- 后向填充
- 線性插值
- 二次插值
- 最近鄰平均
- 季節性差值平均
為了衡量歸因性能，我手動引入時間序列的缺失值，用上述方法進行歸因，然后衡量歸因與實際值的平均平方誤差。
# 生成數據集prcPupdate({'xtick.bottom' : False})## 1. Actual ------------------------------- df_or.plot(style=".-")## 2. 正向填充-------------------------- df_ffill = df.ffill()## 3. 后向填充------------------------- df_bfill = df.bfill()## 4. 線性插值 ------------------ df['rownum'] = np.arange(df.shape[0])## 5. 立體插值-------------------- f2 = interp1d(df_um'], df_nona['v kind='cubic')# 內插法參考。## 6. n "過去最近的鄰居的平均值 ------knnmean(df.value.values, 8)np.round(meansquarerr('), 2)## 7. 季節性平均值 ----------------------------"""計算相應季節性時期的平均值ts: 時間序列的一維數組式n: 時間序列的季節性窗口長度"""out = np.cpy(ts)for i, val in merate(ts):if np.isnan(val):ts_seas = ts[i-1::-n] # 只有以前的季節如果np.isan(np.naean(ts_seas))。ts_seas = np.concenate([ts[i-1::-n], ts[i::n]]) # 以前和以后的out[i] = np.nanan(ts_seas) * lrseasonal_mean(df.value, n=12, lr=1.25)缺失值處理方法
你也可以考慮以下方法，這取決于你希望推斷的準確程度。
16. 什么是自相關和偏自相關函數？
自相關只是一個序列與它自己的滯后期的相關關系。如果一個序列是顯著的自相關，那就意味著，該序列以前的值（滯后）可能有助于預測當前的值。
偏自相關也傳達了類似的信息，但它傳達的是一個序列與其滯后期的純相關，排除了中間滯后期的相關貢獻。
# 計算ACF和PACF直到50個滯后期 acf_50 = acf(value, nlags=50) pacf_50 = pacf(value, nlags=50)# 繪制圖表 fig, axes = plt.subplots(1,2,figsize=（16,3）)ACF和PACF
17. 如何計算偏自相關函數？
那么，如何計算偏自相關？
一個序列的滯后期（k）的偏自相關是該滯后期在Y的自回歸方程中的系數。Y的自回歸方程只不過是以其自身的滯后期為預測因素的Y的線性回歸。
例如，如果Y_t是當前序列，Y_t-1是Y的滯后1，那么滯后3的偏自相關（Y_t-3）是Y_t-3的系數$alpha_3$，在以下方程中。
自回歸方程
18. 滯后圖
滯后圖是一個時間序列與自身滯后的散點圖。它通常是用來檢查自相關的。如果在序列中存在任何像你下面看到的模式，該序列是自相關的。如果沒有這樣的模式，該序列可能是隨機白噪聲。
在下面關于太陽黑子區域時間序列的例子中，隨著n_lag的增加，圖變得越來越分散。
# ＃導入a10 = pd.read_csv('10.csv')# 繪圖 for i, ax in enumerate(axes.flatten()[:4]):fig.suptitle('藥物銷售的滯后圖', y=1.05)藥物銷售滯后曲線圖
滯后圖 太陽黑子
19. 如何估計一個時間序列的可預測性？
一個時間序列的規律性和可重復性越強，就越容易進行預測。近似熵 "可以用來量化一個時間序列中波動的規律性和不可預測性。
近似熵越高，預測就越困難。
另一個更好的替代方法是 "樣本熵"。
樣本熵與近似熵相似，但在估計復雜性方面更加一致，即使是較小的時間序列。例如，一個數據點較少的隨機時間序列的 "近似熵 "可能比一個較 "規則 "的時間序列低，而一個較長的隨機時間序列的 "近似熵 "會更高。
樣本熵很好地處理了這個問題。請看下面的演示。
def AEn(U, m, r):""計算Aproximate entropy""def _maxdit(x_i, x_j):returnmax([abs(u- va) for ua, va in zp(x_i, x_j)])def _phi(m):x = [[U[j]for in range(i, i ] for i in range(N - m + 1)]C = [len([1 for x_jn x if _maist(x_i, x_j) <= r]) / (N - m + 1.0) for xi in x] 。返回 (N - m +.0)**(-1) * sump.log(C))N = len(U) def Samp:""計算樣本熵"""maxstx_i, x_j):retur max([abs(ua - va) fo ua, vain zip(_i x_j)])phi(m)x = [[j]or j i ane(i, i m - 1 + 1] for iin age(N - m + 1)]C= [len([1 for j in rane(len(x)) if i != j and _mist(x[i, xj]) <= r]) for i in rane(ln(x)) ]= en()20. 為什么和如何對時間序列進行平滑處理？
對一個時間序列進行平滑處理可能在以下方面有用。
- 減少信號中噪聲的影響，得到一個經過噪聲過濾的序列的近似值。
- 平滑化后的序列可以作為解釋原始序列本身的一個特征。
- 更好地觀察基本趨勢
那么，如何對一個序列進行平滑處理？讓我們討論一下以下方法。
移動平均數只不過是定義寬度的滾動窗口的平均值。但你必須合理地選擇窗口寬度，因為大的窗口尺寸會使序列過度平滑。例如，窗口尺寸等于季節性持續時間（例如：12個月的序列），將有效地消除季節性效應。
LOESS是 "LOcalized regrESSion "的簡稱，它在每個點的局部附近進行多次回歸。它是在statsmodels軟件包中實現的，你可以用frac參數來控制平滑的程度，該參數指定了附近的數據點的百分比，應該被視為適合回歸模型。
# 導入 pd.read_csv('ele.csv', parse_dates=['date'], index_col='date')# 1. 移動平均數 dma = df_rg.vale.r.man()# 2. 平滑(5%和15%) pd.DtaFame(lowess(dfoig.alu, np.ane(len(d_origale)), fac=0.05)# 繪圖 fig, axes = plt.ubplos(4,1, figsiz=（7, 7, sharex=rue, dp=120) df_ori['aue'].pot(ax=axes0], color='k')平滑化時間序列
如何使用格蘭杰因果測試來了解一個時間序列是否有助于預測另一個時間序列？
格蘭杰因果檢驗是用來確定一個時間序列是否有助于預測另一個時間序列的。
格蘭杰因果關系測試是如何工作的？
它是基于這樣的想法：如果X導致Y，那么基于Y的前值和X的前值對Y的預測應該優于僅基于Y的前值的預測。
因此，理解格蘭杰因果關系不應該被用來測試Y的滯后期是否導致Y。
無效假設是：第二列中的序列不會導致第一列中的序列的格蘭杰。如果P值小于顯著性水平（0.05），那么你就拒絕無效假設，并得出結論：上述X的滯后期確實是有用的。
第二個參數maxlag說的是在測試中應該包括多少個Y的滯后期。
df = pd.rea_csv('a10.csv', parse_dates=['date']) df['moth'] = df.date.dt.nth gragecaslitess(df[['alue', 'moh']], maxag2)在上述情況下，所有檢驗的P值都是零。因此，"月 "數據確實可以用來預測航空乘客。
總結
以上是生活随笔為你收集整理的Python面板时间序列数据预测：格兰杰因果关系检验Granger causality test药品销售实例与可视化的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： Good Luck!
- 下一篇： 管网三维激光扫描建模_BIM建模_可视化