時(shí)間序列的特征分析，一文幫你全Get！

【導(dǎo)語(yǔ)】時(shí)間序列是指以固定時(shí)間為間隔的序列值。本篇教程將教大家用 Python 對(duì)時(shí)間序列進(jìn)行特征分析。

1、什么是時(shí)間序列？

時(shí)間序列是指以固定時(shí)間為間隔的、由所觀察的值組成的序列。根據(jù)觀測(cè)值的不同頻率，可將時(shí)間序列分成小時(shí)、天、星期、月份、季度和年等時(shí)間形式的序列。有時(shí)候，你也可以將秒鐘和分鐘作為時(shí)間序列的間隔，如每分鐘的點(diǎn)擊次數(shù)和訪客數(shù)等等。

為什么我們要對(duì)時(shí)間序列進(jìn)行分析呢？

因?yàn)楫?dāng)你想對(duì)一個(gè)序列進(jìn)行預(yù)測(cè)時(shí)，首先要完成分析這個(gè)步驟。除此之外，時(shí)間序列的預(yù)測(cè)也具有極大商業(yè)價(jià)值，如企業(yè)的供求量、網(wǎng)站的訪客量以及股票價(jià)格等，都是極其重要的時(shí)間序列數(shù)據(jù)。

那么，時(shí)間序列分析都包括哪些內(nèi)容呢？

要做好時(shí)間序列分析，必須要理解序列的內(nèi)在屬性，這樣才能做出更有意義且精準(zhǔn)的預(yù)測(cè)。

2、如何在 Python 中引入時(shí)間序列？

關(guān)于時(shí)間序列的數(shù)據(jù)大都存儲(chǔ)在 csv 文件或其他形式的表格文件里，且都包含兩個(gè)列：日期和觀測(cè)值。

首先我們來(lái)看 panda 包里面的 read_csv() 函數(shù)，它可以將時(shí)間序列數(shù)據(jù)集（關(guān)于澳大利亞藥物銷售的 csv 文件）讀取為 pandas 數(shù)據(jù)框。增加一個(gè) parse_dates=['date'] 字段，可以把包含日期的數(shù)據(jù)列解析為日期字段。

from dateutil.parser import parse import matplotlib as mplimport matplotlib.pyplot as pltimport seaborn as snsimport numpy as npimport pandas as pdplt.rcParams.update({'figure.figsize': (10, 7), 'figure.dpi': 120})# Import as Dataframedf = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv', parse_dates=['date'])df.head()

時(shí)間序列數(shù)據(jù)框

此外，你也可以將文件讀取為 pandas 序列，把日期作為索引列，只需在 pd.read_csv() 中指定 index_col 參數(shù)。

ser = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv', parse_dates=['date'], index_col='date')ser.head()

pandas 序列

注意，在 pandas 序列中，'value' 列的位置高于 'date' 列，這表明它是一個(gè) pandas 序列而非數(shù)據(jù)框。

3、什么是面板數(shù)據(jù)？

面板數(shù)據(jù)同樣是基于時(shí)間的數(shù)據(jù)集。

不同之處是，除了時(shí)間序列，面板數(shù)據(jù)還包括一個(gè)或多個(gè)相關(guān)變量，這些變量也是在同個(gè)時(shí)間段內(nèi)測(cè)得的。

面板數(shù)據(jù)中的列包括有助于預(yù)測(cè) y 值的解釋變量，這些特征列可用于之后的預(yù)測(cè)。以下是關(guān)于面板數(shù)據(jù)的例子：

# dataset source: https://github.com/rouseguydf = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/MarketArrivals.csv')df = df.loc[df.market=='MUMBAI', :]df.head()

面板數(shù)據(jù)

4、時(shí)間序列可視化

接下來(lái)，我們用 matplotlib 對(duì)序列進(jìn)行可視化。

# Time series data source: fpp pacakge in R.import matplotlib.pyplot as pltdf = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv', parse_dates=['date'], index_col='date')#DrawPlotdef plot_df(df, x, y, title="", xlabel='Date', ylabel='Value', dpi=100): plt.figure(figsize=(16,5), dpi=dpi) plt.plot(x, y, color='tab:red') plt.gca().set(title=title, xlabel=xlabel, ylabel=ylabel)plt.show()plot_df(df,x=df.index,y=df.value,title='Monthlyanti-diabeticdrugsalesinAustraliafrom1992to2008.')

時(shí)間序列可視化

由于所有的值都為正數(shù)，無(wú)論使用 y 軸的哪一側(cè)都可以看到增長(zhǎng)的趨勢(shì)。

# Import datadf = pd.read_csv('datasets/AirPassengers.csv', parse_dates=['date'])x = df['date'].valuesy1 = df['value'].values# Plotfig, ax = plt.subplots(1, 1, figsize=(16,5), dpi= 120)plt.fill_between(x, y1=y1, y2=-y1, alpha=0.5, linewidth=2, color='seagreen')plt.ylim(-800, 800)plt.title('Air Passengers (Two Side View)', fontsize=16)plt.hlines(y=0, xmin=np.min(df.date), xmax=np.max(df.date), linewidth=.5)plt.show()

飛機(jī)乘客數(shù)據(jù) - 雙邊序列

由于這是一個(gè)月份的時(shí)間序列，每年的走勢(shì)都遵循著特定重復(fù)的模式，你可以在同一個(gè)圖中畫出單獨(dú)每年的折線。這樣有助于對(duì)這幾年的趨勢(shì)走向進(jìn)行平行對(duì)比。

季節(jié)性時(shí)間序列的可視化：

#ImportDatadf = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv', parse_dates=['date'], index_col='date')df.reset_index(inplace=True)# Prepare datadf['year'] = [d.year for d in df.date]df['month'] = [d.strftime('%b') for d in df.date]years = df['year'].unique()# Prep Colorsnp.random.seed(100)mycolors = np.random.choice(list(mpl.colors.XKCD_COLORS.keys()), len(years), replace=False)# Draw Plotplt.figure(figsize=(16,12), dpi= 80)for i, y in enumerate(years): if i > 0: plt.plot('month', 'value', data=df.loc[df.year==y, :], color=mycolors[i], label=y) plt.text(df.loc[df.year==y, :].shape[0]-.9, df.loc[df.year==y, 'value'][-1:].values[0], y, fontsize=12, color=mycolors[i])# Decorationplt.gca().set(xlim=(-0.3, 11), ylim=(2, 30), ylabel='$Drug Sales$', xlabel='$Month$')plt.yticks(fontsize=12, alpha=.7)plt.title("Seasonal Plot of Drug Sales Time Series", fontsize=20)plt.sho

藥品銷售的季節(jié)性時(shí)間序列圖

每年的二月份，藥品銷售會(huì)有一次大幅度下跌，三月份會(huì)回漲，四月份再次下跌，以此規(guī)律循環(huán)。很明顯，該模式以年為單位，每年循環(huán)往復(fù)。

不過(guò)，隨著年份的變化，藥品銷售總體呈上升趨勢(shì)。你可以選擇使用箱型圖將這一趨勢(shì)進(jìn)行可視化，可以方便看出每一年的變化。同樣地，你也可以按月份繪制箱型圖，來(lái)觀察每個(gè)月的變化。

按月份（季節(jié)）和年份繪制箱型圖：你可以將數(shù)據(jù)處理成以季節(jié)為時(shí)間間隔，然后觀察特定年份內(nèi)值的分布，也可以將全部時(shí)間的數(shù)據(jù)進(jìn)行對(duì)比。

# Import Datadf = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv', parse_dates=['date'], index_col='date')df.reset_index(inplace=True)# Prepare datadf['year'] = [d.year for d in df.date]df['month'] = [d.strftime('%b') for d in df.date]years = df['year'].unique()# Draw Plotfig, axes = plt.subplots(1, 2, figsize=(20,7), dpi= 80)sns.boxplot(x='year', y='value', data=df, ax=axes[0])sns.boxplot(x='month', y='value', data=df.loc[~df.year.isin([1991, 2008]), :])# Set Titleaxes[0].set_title('Year-wise Box Plot\n(The Trend)', fontsize=18);axes[1].set_title('Month-wise Box Plot\n(The Seasonality)', fontsize=18)plt.show

年份序列和月份序列的箱型圖

上面的箱型圖可以使年份和月份的序列更易于觀察。同樣，在月份的箱線圖中，十二月和一月的藥品銷售額明顯更高，這是因?yàn)樘幱诩倨谡劭奂尽?/p>

到目前為止，我們了解了通過(guò)相似性來(lái)觀察時(shí)間序列的模式，接下來(lái)，如何發(fā)現(xiàn)這些常見(jiàn)模式中的差異呢？

5、時(shí)間序列的模式

任何一個(gè)時(shí)間序列都可以被分解成以下幾個(gè)部分：基準(zhǔn) + 趨勢(shì)成分 + 季節(jié)成分 + 殘差成分。

趨勢(shì)是指時(shí)間序列中上升或下降的傾斜程度。但季節(jié)成分是由于受季節(jié)因素影響而產(chǎn)生的周期性模式循環(huán)，也可能受每年內(nèi)不同月份、每月內(nèi)不同日期、工作日或周末，甚至每天內(nèi)不同時(shí)間的影響。

然而，不一定所有的時(shí)間序列都具備趨勢(shì)或季節(jié)性。一個(gè)時(shí)間序列也可能不存在趨勢(shì)，但具有季節(jié)性。反之亦然。

因此，一個(gè)時(shí)間序列可以被想象成趨勢(shì)、季節(jié)性和殘差項(xiàng)的組合。

fig, axes = plt.subplots(1,3, figsize=(20,4), dpi=100)pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/guinearice.csv', parse_dates=['date'], index_col='date').plot(title='Trend Only', legend=False, ax=axes[0])pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/sunspotarea.csv', parse_dates=['date'], index_col='date').plot(title='Seasonality Only', legend=False, ax=axes[1])pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/AirPassengers.csv',parse_dates=['date'],index_col='date').plot(title='TrendandSeasonality',legend=False,ax=axes[2]

時(shí)間序列的模式

另一個(gè)需要考慮的方面是周期性模式。當(dāng)序列中的上升和下降，不是按日歷中的特定時(shí)間間隔發(fā)生時(shí)，就會(huì)出現(xiàn)這種情況。注意不要把“周期”作用和“季節(jié)”作用混淆。

那么，如何區(qū)分“周期”和“季節(jié)”呢？

如果序列中的模式不是以日歷中特定間隔循環(huán)出現(xiàn)的，那么就是周期。因?yàn)榕c季節(jié)性不同，周期作用通常受到商業(yè)或社會(huì)經(jīng)濟(jì)等因素的影響。

6、加法與乘法時(shí)間序列

根據(jù)趨勢(shì)和季節(jié)的固有屬性，一個(gè)時(shí)間序列可以被建模為加法模型或乘法模型，也就是說(shuō)，序列中的值可以用各個(gè)成分的加和或乘積來(lái)表示：

加法時(shí)間序列：

值 = 基準(zhǔn) + 趨勢(shì) + 季節(jié) + 殘差

乘法時(shí)間序列：

值 = 基準(zhǔn) x 趨勢(shì) x 季節(jié) x 殘差

7、如何將時(shí)間序列的成分分解出來(lái)？

通過(guò)將一個(gè)時(shí)間序列視為基準(zhǔn)、趨勢(shì)、季節(jié)指數(shù)及殘差的加法或乘法組合，你可以對(duì)時(shí)間序列進(jìn)行經(jīng)典分解。

statsmodels 的 seasonal_decompose 函數(shù)可以使這一過(guò)程非常容易。from statsmodels.tsa.seasonal import seasonal_decomposefrom dateutil.parser import parse# Import Datadf = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv', parse_dates=['date'], index_col='date')# Multiplicative Decompositionresult_mul = seasonal_decompose(df['value'], model='multiplicative', extrapolate_trend='freq')# Additive Decompositionresult_add = seasonal_decompose(df['value'], model='additive', extrapolate_trend='freq')# Plotplt.rcParams.update({'figure.figsize': (10,10)})result_mul.plot().suptitle('Multiplicative Decompose', fontsize=22)result_add.plot().suptitle('Additive Decompose', fontsize=22)plt.sho

加法和乘法分解

設(shè)置 extrapolate_trend='freq' 有助于處理序列首部趨勢(shì)和殘差中的空值。

如果你仔細(xì)觀察加法分解中的殘差項(xiàng)，會(huì)發(fā)現(xiàn)其中仍保留了一些模式。然而，乘法分解中的殘差項(xiàng)看起來(lái)更具有隨機(jī)性。因此，對(duì)于這一特定序列來(lái)說(shuō)，采用乘法分解更合適。

趨勢(shì)、季節(jié)和殘差成分的數(shù)值輸出均存儲(chǔ)在 result_mul 里，下面我們對(duì)其進(jìn)行提取，并放入數(shù)據(jù)框中：

# Extract the Components ----# Actual Values = Product of (Seasonal * Trend * Resid)df_reconstructed = pd.concat([result_mul.seasonal, result_mul.trend, result_mul.resid, result_mul.observed], axis=1)df_reconstructed.columns = ['seas', 'trend', 'resid', 'actual_values']df_reconstructed.head()

仔細(xì)觀察，會(huì)發(fā)現(xiàn) seas、trend 和 resid 三列的乘積正好等于 actual_values。

8、平穩(wěn)和非平穩(wěn)時(shí)間序列

平穩(wěn)是時(shí)間序列的屬性之一。平穩(wěn)序列中的值不是時(shí)間的函數(shù)。

也就是說(shuō)，平穩(wěn)序列的平均值、方差和自相關(guān)性等統(tǒng)計(jì)特征始終為常數(shù)。序列的自相關(guān)性是指該序列與之前的值間的相關(guān)性。

平穩(wěn)的時(shí)間序列也不包括季節(jié)因素的影響。那么如何分辨一個(gè)序列是否平穩(wěn)呢？讓我們來(lái)看幾個(gè)例子：

平穩(wěn)和非平穩(wěn)時(shí)間序列

上圖來(lái)自 R 語(yǔ)言的 TSTutorial 教程。

為什么說(shuō)序列平穩(wěn)很重要呢？

我會(huì)對(duì)此稍微做一些解釋，但要清楚一點(diǎn)，通過(guò)采用合適的變換，我們近乎可以將任何時(shí)間序列變得平穩(wěn)。大多數(shù)統(tǒng)計(jì)預(yù)測(cè)方法最初都是為平穩(wěn)時(shí)間序列而設(shè)計(jì)的。預(yù)測(cè)過(guò)程的第一步是通過(guò)一些變換，來(lái)將非平穩(wěn)序列變成平穩(wěn)序列。

9、如何將時(shí)間序列變平穩(wěn)？

你可以通過(guò)以下幾種方式得到平穩(wěn)序列：

求序列的差分

求序列的 log 值

求序列的 n 次方根

把上面三種方法相結(jié)合

將時(shí)間序列平穩(wěn)化最普遍且便捷的方法是對(duì)序列進(jìn)行差分運(yùn)算，至少執(zhí)行一次，直到序列趨于平穩(wěn)。

那么什么是差分呢？

如果 Y_t 為時(shí)間 't' 對(duì)應(yīng)的值，那么第一個(gè)差分值為 Y = Yt – Yt-1。簡(jiǎn)單來(lái)說(shuō)，對(duì)序列進(jìn)行差分運(yùn)算就是用下一個(gè)值減去當(dāng)前值。如果第一次差分不能使序列平穩(wěn)，你可以嘗試做第二次差分，直到符合要求。

舉個(gè)例子，有這樣一個(gè)序列：[1, 5, 2, 12, 20]

第一次差分運(yùn)算的結(jié)果為：[5-1, 2-5, 12-2, 20-12] = [4, -3, 10, 8]

第二次差分運(yùn)算的結(jié)果為：[-3-4, -10-3, 8-10] = [-7, -13, -2]

10、為什么要在預(yù)測(cè)前將非平穩(wěn)序列變平穩(wěn)？

對(duì)平穩(wěn)序列進(jìn)行預(yù)測(cè)要相對(duì)容易一些，預(yù)測(cè)結(jié)果也更可信。其中一個(gè)重要原因是，自回歸預(yù)測(cè)模型本質(zhì)上是線性回歸模型，將序列自身的滯后作為預(yù)測(cè)因子。

如果預(yù)測(cè)因子之間互不相關(guān)，線性回歸的效果最優(yōu)。那么將序列平穩(wěn)化就可以解決這一問(wèn)題，因?yàn)樗梢匀コ魏纬志玫淖韵嚓P(guān)性，所以可以使預(yù)測(cè)模型中的預(yù)測(cè)因子近乎獨(dú)立。

現(xiàn)在我們知道了使序列平穩(wěn)化的重要性，那么應(yīng)該如何檢查一個(gè)序列是否平穩(wěn)呢？

11、如何測(cè)試平穩(wěn)性？

我們可以像之前那樣，通過(guò)繪制序列圖來(lái)觀察一個(gè)序列是否平穩(wěn)。

另一種方法是將序列分解成兩個(gè)或多個(gè)連續(xù)的部分，并求其統(tǒng)計(jì)值，如平均值、方差和自相關(guān)系數(shù)。如果這些統(tǒng)計(jì)值間的差異很大，那么該序列大概率不是平穩(wěn)序列。

盡管如此，你仍需要一種方法來(lái)定量地判斷某個(gè)序列是否平穩(wěn)。一個(gè)名為“單位根檢驗(yàn)”的統(tǒng)計(jì)檢驗(yàn)方法可以解決這一問(wèn)題。

單位根檢驗(yàn)有如下幾種方法：

ADF Test (Augmented Dickey Fuller test)

KPSS Test (Kwiatkowski-Phillips-Schmidt-Shin) — 趨勢(shì)平穩(wěn)

PP Test (Philips Perron test)

最常用的方法是 ADF test，零假設(shè)為：用時(shí)間序列對(duì)單位根進(jìn)行處理，結(jié)果是不平穩(wěn)。如果 P 值小于顯著性水平 0.05，則拒絕零假設(shè)，即不成立。

另一方面，KPSS test 可用來(lái)檢測(cè)趨勢(shì)平穩(wěn)性。零假設(shè)與 P 值含義都與 ADH test 相反。以下代碼基于 Python 的 statsmodels 包執(zhí)行了兩種檢測(cè)方法：

from statsmodels.tsa.stattools import adfuller, kpssdf = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv', parse_dates=['date'])# ADF Testresult = adfuller(df.value.values, autolag='AIC')print(f'ADF Statistic: {result[0]}')print(f'p-value: {result[1]}')for key, value in result[4].items(): print('Critial Values:') print(f' {key}, {value}')# KPSS Testresult = kpss(df.value.values, regression='c')print('\nKPSS Statistic: %f' % result[0])print('p-value: %f' % result[1])for key, value in result[3].items(): print('Critial Values:') print(f' {key}, {value}'

from statsmodels.tsa.stattools import adfuller, kpssdf = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv', parse_dates=['date'])# ADF Testresult = adfuller(df.value.values, autolag='AIC')print(f'ADF Statistic: {result[0]}')print(f'p-value: {result[1]}')for key, value in result[4].items(): print('Critial Values:') print(f' {key}, {value}')# KPSS Testresult = kpss(df.value.values, regression='c')print('\nKPSS Statistic: %f' % result[0])print('p-value: %f' % result[1])for key, value in result[3].items(): print('Critial Values:') print(f' {key}, {value}'

12、白噪聲與平穩(wěn)序列的差別是什么？

和平穩(wěn)序列一樣，白噪聲也是關(guān)于時(shí)間的函數(shù)，它的均值和方差始終不變。但區(qū)別在于，白噪聲是完全隨機(jī)的，且均值為零。

白噪聲沒(méi)有任何模式。如果你把調(diào)頻廣播的聲波訊號(hào)想象成時(shí)間序列，調(diào)頻道時(shí)的空白聲音就是白噪聲。

從數(shù)學(xué)上來(lái)講，一個(gè)完全隨機(jī)且均值為零的序列就是白噪聲。

randvals = np.random.randn(1000)pd.Series(randvals).plot(title='RandomWhiteNoise',color='k')

隨機(jī)白噪聲

13、如何對(duì)時(shí)間序列去趨勢(shì)？

對(duì)時(shí)間序列去趨勢(shì)，是指去除序列中的趨勢(shì)成分。但要如何提取趨勢(shì)成分呢？有以下幾種方法：

減去與時(shí)間序列擬合程度最好的曲線。這條最優(yōu)曲線可由線性回歸模型獲得，時(shí)間步長(zhǎng)作為預(yù)測(cè)因子。如需處理更復(fù)雜的趨勢(shì)，你可以嘗試在模型中使用二次項(xiàng) (x^2)。

減去由時(shí)間序列分解得到的趨勢(shì)成分。

減去均值。

使用過(guò)濾器，如 Baxter-King (statsmodels.tsa.filters.bkfilter) 或 Hodrick-Prescott (statsmodels.tsa.filters.hpfilter)，來(lái)去除移動(dòng)平均趨勢(shì)線或周期性成分。

下面我們來(lái)執(zhí)行這兩種方法：

# Using scipy: Subtract the line of best fitfrom scipy import signaldf = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv', parse_dates=['date'])detrended = signal.detrend(df.value.values)plt.plot(detrended)plt.title('DrugSalesdetrendedbysubtractingtheleastsquaresfit',fontsize=16)

通過(guò)減掉最小二乘擬合線對(duì)時(shí)間序列去趨勢(shì)

# Using statmodels: Subtracting the Trend Component.from statsmodels.tsa.seasonal import seasonal_decomposedf = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv', parse_dates=['date'], index_col='date')result_mul = seasonal_decompose(df['value'], model='multiplicative', extrapolate_trend='freq')detrended = df.value.values - result_mul.trendplt.plot(detrended)plt.title('Drug Sales detrended by subtracting the trend component', fontsize=16)

通過(guò)減掉趨勢(shì)成分對(duì)時(shí)間序列去趨勢(shì)

14、如何對(duì)時(shí)間序列去季節(jié)性？

對(duì)時(shí)間序列去季節(jié)性同樣有多種方法，如下：

把特定長(zhǎng)度的移動(dòng)平均值作為季節(jié)窗口。

對(duì)序列做季節(jié)性差分（用當(dāng)前值減去上個(gè)季度的值）。

用當(dāng)前序列除以由 STL 分解得到的季節(jié)指數(shù)。

若除以季節(jié)指數(shù)的效果不好，可以嘗試取序列的對(duì)數(shù)，然后對(duì)其去季節(jié)。之后你可以通過(guò)指數(shù)運(yùn)算來(lái)恢復(fù)原來(lái)的值。

# Subtracting the Trend Component.df = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv', parse_dates=['date'], index_col='date')# Time Series Decompositionresult_mul = seasonal_decompose(df['value'], model='multiplicative', extrapolate_trend='freq')# Deseasonalizedeseasonalized = df.value.values / result_mul.seasonal# Plotplt.plot(deseasonalized)plt.title('Drug Sales Deseasonalized', fontsize=16)plt.plot

對(duì)時(shí)間序列去季節(jié)

15、如何檢測(cè)時(shí)間序列的季節(jié)性？

一般方法是畫出序列圖，觀察固定的時(shí)間間隔是否有重復(fù)模式出現(xiàn)。因此，季節(jié)性的類型由時(shí)鐘或日歷決定：

一天中的小時(shí)

月份中的日期

星期

月份

年份

不過(guò)，如果你想對(duì)季節(jié)性做一個(gè)明確的檢驗(yàn)，可以使用自相關(guān)函數(shù) (ACF) 圖，接下來(lái)的部分會(huì)做相關(guān)詳細(xì)介紹。當(dāng)季節(jié)模式明顯時(shí)，ACF 圖中季節(jié)窗口的整數(shù)倍處會(huì)反復(fù)出現(xiàn)特定的尖峰。

例如，藥品銷售的時(shí)間序列是月份序列，每年會(huì)出現(xiàn)重復(fù)的模式，你會(huì)在第 12、24、36 個(gè)序列值處看到尖峰。

必須要提醒你的是，現(xiàn)實(shí)生活中的數(shù)據(jù)集很少會(huì)出現(xiàn)特別明顯的模式，可能會(huì)被一些噪聲污染，所以需要更加仔細(xì)觀察其中的模式。

from pandas.plotting import autocorrelation_plotdf = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv')# Draw Plotplt.rcParams.update({'figure.figsize':(9,5), 'figure.dpi':120})autocorrelation_plot(df.value.tolist())

自相關(guān)系數(shù)圖

16、如果處理時(shí)間序列中的缺失值？

有時(shí)候，時(shí)間序列中會(huì)出現(xiàn)缺失的值或日期。這意味著，某些數(shù)據(jù)沒(méi)有獲取到，或者無(wú)法對(duì)這些時(shí)間段進(jìn)行觀測(cè)。也可能那些時(shí)間的測(cè)量值本身為零，這種情況下你只需對(duì)其填充零。

第二種情況，你不應(yīng)該直接用序列的均值對(duì)缺失處進(jìn)行填充，尤其當(dāng)該序列不是平穩(wěn)序列時(shí)。比較暴力但有效的解決方法是用前一個(gè)值來(lái)填充缺失處。

根據(jù)序列的內(nèi)在屬性，你可以嘗試多種方法。以下是幾種比較有效的填充方法：

向后填充法

線性插值法

二次插值法

最近鄰均值法

季節(jié)均值法

為了評(píng)估缺失值的填充效果，我在時(shí)間序列中手動(dòng)加入缺失值，用以上幾種方法對(duì)其進(jìn)行填充，然后計(jì)算填充后的序列與原序列的均方誤差。

# # Generate datasetfrom scipy.interpolate import interp1dfrom sklearn.metrics import mean_squared_errordf_orig = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv', parse_dates=['date'], index_col='date').head(100)df = pd.read_csv('datasets/a10_missings.csv', parse_dates=['date'], index_col='date')fig, axes = plt.subplots(7, 1, sharex=True, figsize=(10, 12))plt.rcParams.update({'xtick.bottom' : False})## 1. Actual -------------------------------df_orig.plot(title='Actual', ax=axes[0], label='Actual', color='red', style=".-")df.plot(title='Actual', ax=axes[0], label='Actual', color='green', style=".-")axes[0].legend(["Missing Data", "Available Data"])## 2. Forward Fill --------------------------df_ffill = df.ffill()error = np.round(mean_squared_error(df_orig['value'], df_ffill['value']), 2)df_ffill['value'].plot(title='Forward Fill (MSE: ' + str(error) +")", ax=axes[1], label='Forward Fill', style=".-")## 3. Backward Fill -------------------------df_bfill = df.bfill()error = np.round(mean_squared_error(df_orig['value'], df_bfill['value']), 2)df_bfill['value'].plot(title="Backward Fill (MSE: " + str(error) +")", ax=axes[2], label='Back Fill', color='firebrick', style=".-")## 4. Linear Interpolation ------------------df['rownum'] = np.arange(df.shape[0])df_nona = df.dropna(subset = ['value'])f = interp1d(df_nona['rownum'], df_nona['value'])df['linear_fill'] = f(df['rownum'])error = np.round(mean_squared_error(df_orig['value'], df['linear_fill']), 2)df['linear_fill'].plot(title="Linear Fill (MSE: " + str(error) +")", ax=axes[3], label='Cubic Fill', color='brown', style=".-")## 5. Cubic Interpolation --------------------f2 = interp1d(df_nona['rownum'], df_nona['value'], kind='cubic')df['cubic_fill'] = f2(df['rownum'])error = np.round(mean_squared_error(df_orig['value'], df['cubic_fill']), 2)df['cubic_fill'].plot(title="Cubic Fill (MSE: " + str(error) +")", ax=axes[4], label='Cubic Fill', color='red', style=".-")# Interpolation References:# https://docs.scipy.org/doc/scipy/reference/tutorial/interpolate.html# https://docs.scipy.org/doc/scipy/reference/interpolate.html## 6. Mean of 'n' Nearest Past Neighbors ------def knn_mean(ts, n): out = np.copy(ts) for i, val in enumerate(ts): if np.isnan(val): n_by_2 = np.ceil(n/2) lower = np.max([0, int(i-n_by_2)]) upper = np.min([len(ts)+1, int(i+n_by_2)]) ts_near = np.concatenate([ts[lower:i], ts[i:upper]]) out[i] = np.nanmean(ts_near) return outdf['knn_mean'] = knn_mean(df.value.values, 8)error = np.round(mean_squared_error(df_orig['value'], df['knn_mean']), 2)df['knn_mean'].plot(title="KNN Mean (MSE: " + str(error) +")", ax=axes[5], label='KNN Mean', color='tomato', alpha=0.5, style=".-")## 7. Seasonal Mean ----------------------------def seasonal_mean(ts, n, lr=0.7): """ Compute the mean of corresponding seasonal periods ts: 1D array-like of the time series n: Seasonal window length of the time series """ out = np.copy(ts) for i, val in enumerate(ts): if np.isnan(val): ts_seas = ts[i-1::-n] # previous seasons only if np.isnan(np.nanmean(ts_seas)): ts_seas = np.concatenate([ts[i-1::-n], ts[i::n]]) # previous and forward out[i] = np.nanmean(ts_seas) * lr return outdf['seasonal_mean'] = seasonal_mean(df.value, n=12, lr=1.25)error = np.round(mean_squared_error(df_orig['value'], df['seasonal_mean']), 2)df['seasonal_mean'].plot(title="Seasonal Mean (MSE: " + str(error) +")", ax=axes[6], label='Seasonal Mean', color='blue', alpha=0.5, s

缺失值處理

17、什么是自相關(guān)和偏自相關(guān)函數(shù)？

簡(jiǎn)單來(lái)說(shuō)，自相關(guān)就是一個(gè)序列的值與自己本身具有相關(guān)性。如果一個(gè)序列呈現(xiàn)顯著自相關(guān)，意味著序列的前一個(gè)值可用于預(yù)測(cè)當(dāng)前值。

偏自相關(guān)也傳達(dá)了類似的信息，但更偏重于序列與自身滯后序列的相關(guān)性，消除了由于較短滯后所導(dǎo)致的任何相關(guān)性。

from statsmodels.tsa.stattools import acf, pacffrom statsmodels.graphics.tsaplots import plot_acf, plot_pacfdf = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv')# Calculate ACF and PACF upto 50 lags# acf_50 = acf(df.value, nlags=50)# pacf_50 = pacf(df.value, nlags=50)# Draw Plotfig, axes = plt.subplots(1,2,figsize=(16,3), dpi= 100)plot_acf(df.value.tolist(), lags=50, ax=axes[0])plot_pacf(df.value.tolist(),lags=50,ax=axes[1

ACF 和 PACF

18、如何計(jì)算偏自相關(guān)系數(shù)？

序列滯后 k 處的偏自相關(guān)系數(shù)是 Y 的自回歸方程的滯后系數(shù)。Y 的自回歸方程是指 Y 以自己的滯后作為預(yù)測(cè)因子的線性回歸。

舉個(gè)例子，如果 Y_t 為當(dāng)前序列，Y_t-1 即為滯后期為 1 的 Y 值，那么滯后期為 3 處 (Y_t-3) 的偏自相關(guān)系數(shù)是下面方程中的 α3：

自回歸方程

19、滯后圖

滯后圖是時(shí)間序列與自身滯后的分布圖，通常用來(lái)檢驗(yàn)自相關(guān)性。如果序列中出現(xiàn)下圖中的模式，那么說(shuō)明該序列具有自相關(guān)性。如果沒(méi)有出現(xiàn)這些模型，該序列很可能為隨機(jī)白噪聲。

在下面太陽(yáng)黑子區(qū)的例子中，隨著滯后期的增長(zhǎng)，圖中的點(diǎn)越來(lái)越分散。

from pandas.plotting import lag_plotplt.rcParams.update({'ytick.left' : False, 'axes.titlepad':10})# Importss = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/sunspotarea.csv')a10 = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv')# Plotfig, axes = plt.subplots(1, 4, figsize=(10,3), sharex=True, sharey=True, dpi=100)for i, ax in enumerate(axes.flatten()[:4]): lag_plot(ss.value, lag=i+1, ax=ax, c='firebrick') ax.set_title('Lag ' + str(i+1))fig.suptitle('Lag Plots of Sun Spots Area \n(Points get wide and scattered with increasing lag -> lesser correlation)\n', y=1.15)fig, axes = plt.subplots(1, 4, figsize=(10,3), sharex=True, sharey=True, dpi=100)for i, ax in enumerate(axes.flatten()[:4]): lag_plot(a10.value, lag=i+1, ax=ax, c='firebrick') ax.set_title('Lag ' + str(i+1))fig.suptitle('Lag Plots of Drug Sales', y=1.05)plt.sh

藥品銷售滯后圖

太陽(yáng)黑子滯后圖

20、如何評(píng)估時(shí)間序列的可預(yù)測(cè)性？

一個(gè)時(shí)間序列的模式越有規(guī)律，就越容易預(yù)測(cè)。可以用近似熵來(lái)量化時(shí)間序列的規(guī)律性和波動(dòng)的不可預(yù)測(cè)性。

近似熵越高，意味著預(yù)測(cè)難度越大。

另一個(gè)選擇是樣本熵。

樣本熵與近似熵類似，但在不同的復(fù)雜度上更具有一致性，即使是小型時(shí)間序列。例如，相比于“有規(guī)律”的時(shí)間序列，一個(gè)數(shù)據(jù)點(diǎn)較少的隨機(jī)時(shí)間序列的近似熵較低，但一個(gè)較長(zhǎng)的隨機(jī)時(shí)間序列卻有較高的近似熵。

因此，樣本熵更適于解決該問(wèn)題。

# https://en.wikipedia.org/wiki/Approximate_entropyss = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/sunspotarea.csv')a10 = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv')rand_small = np.random.randint(0, 100, size=36)rand_big = np.random.randint(0, 100, size=136)def ApEn(U, m, r): """Compute Aproximate entropy""" def _maxdist(x_i, x_j): return max([abs(ua - va) for ua, va in zip(x_i, x_j)]) def _phi(m): x = [[U[j] for j in range(i, i + m - 1 + 1)] for i in range(N - m + 1)] C = [len([1 for x_j in x if _maxdist(x_i, x_j) <= r]) / (N - m + 1.0) for x_i in x] return (N - m + 1.0)**(-1) * sum(np.log(C)) N = len(U) return abs(_phi(m+1) - _phi(m))print(ApEn(ss.value, m=2, r=0.2*np.std(ss.value))) # 0.651print(ApEn(a10.value, m=2, r=0.2*np.std(a10.value))) # 0.537print(ApEn(rand_small, m=2, r=0.2*np.std(rand_small))) # 0.143print(ApEn(rand_big,?m=2,?r=0.2*np.std(rand_big)))?????#?0.?

0.65147049703335340.53747752249734890.08983769407988440.7369242960384561

# https://en.wikipedia.org/wiki/Sample_entropydef SampEn(U, m, r): """Compute Sample entropy""" def _maxdist(x_i, x_j): return max([abs(ua - va) for ua, va in zip(x_i, x_j)]) def _phi(m): x = [[U[j] for j in range(i, i + m - 1 + 1)] for i in range(N - m + 1)] C = [len([1 for j in range(len(x)) if i != j and _maxdist(x[i], x[j]) <= r]) for i in range(len(x))] return sum(C) N = len(U) return -np.log(_phi(m+1) / _phi(m))print(SampEn(ss.value, m=2, r=0.2*np.std(ss.value))) # 0.78print(SampEn(a10.value, m=2, r=0.2*np.std(a10.value))) # 0.41print(SampEn(rand_small, m=2, r=0.2*np.std(rand_small))) # 1.79print(SampEn(rand_big, m=2, r=0.2*np.std(rand_big))) # 2.

0.78533113663800390.41887013457621214inf2.181224235989778del sys.path[0]

21、為什么要對(duì)時(shí)間序列做平滑處理？如果操作？

對(duì)時(shí)間序列做平滑處理有以下幾個(gè)用處：

減少噪聲影響，從而得到過(guò)濾掉噪聲的、更加真實(shí)的序列。

平滑處理后的序列可用作特征，更好地解釋序列本身。

可以更好地觀察序列本身的趨勢(shì)。

那么如果進(jìn)行平滑處理呢？現(xiàn)討論以下幾種方法：

取移動(dòng)平均線

做 LOESS 平滑（局部回歸）

做 LOWESS 平滑（局部加權(quán)回歸）

移動(dòng)平均是指對(duì)一個(gè)滾動(dòng)的窗口計(jì)算其平均值，該窗口的寬度固定不變。但你必須謹(jǐn)慎選擇窗口寬度，因?yàn)榇翱谶^(guò)寬會(huì)導(dǎo)致序列平滑過(guò)度。例如，如果窗口寬度等于季節(jié)長(zhǎng)度，就會(huì)消除掉季節(jié)因素的作用。

LOESS 是 LOcalized regrESSion 的縮寫，該方法會(huì)在每個(gè)點(diǎn)的局部近鄰點(diǎn)做多次回歸擬合。此處可以使用 statsmodels 包，你可以使用參數(shù) frac 控制平滑程度，即決定周圍多少個(gè)點(diǎn)參與回歸模型的擬合。

from statsmodels.nonparametric.smoothers_lowess import lowessplt.rcParams.update({'xtick.bottom' : False, 'axes.titlepad':5})# Importdf_orig = pd.read_csv('datasets/elecequip.csv', parse_dates=['date'], index_col='date')# 1. Moving Averagedf_ma = df_orig.value.rolling(3, center=True, closed='both').mean()# 2. Loess Smoothing (5% and 15%)df_loess_5 = pd.DataFrame(lowess(df_orig.value, np.arange(len(df_orig.value)), frac=0.05)[:, 1], index=df_orig.index, columns=['value'])df_loess_15 = pd.DataFrame(lowess(df_orig.value, np.arange(len(df_orig.value)), frac=0.15)[:, 1], index=df_orig.index, columns=['value'])# Plotfig, axes = plt.subplots(4,1, figsize=(7, 7), sharex=True, dpi=120)df_orig['value'].plot(ax=axes[0], color='k', title='Original Series')df_loess_5['value'].plot(ax=axes[1], title='Loess Smoothed 5%')df_loess_15['value'].plot(ax=axes[2], title='Loess Smoothed 15%')df_ma.plot(ax=axes[3], title='Moving Average (3)')fig.suptitle('How to Smoothen a Time Series', y=0.95, fontsize=14)plt.sho

時(shí)間序列的平滑處理

22、如何用格蘭杰因果關(guān)系檢驗(yàn)判斷一個(gè)時(shí)間序列是否可以預(yù)測(cè)另一個(gè)？

格蘭杰因果關(guān)系檢驗(yàn)可用作檢測(cè)一個(gè)時(shí)間序列是否可以用來(lái)預(yù)測(cè)另一個(gè)序列。那么格蘭杰因果關(guān)系檢驗(yàn)是如何進(jìn)行運(yùn)算的呢？

該檢驗(yàn)基于一個(gè)假設(shè)，即 X 導(dǎo)致了 Y 的發(fā)生，那么基于 Y 的先前值與 X 的先前值得到的 Y 的預(yù)測(cè)值，要優(yōu)于僅根據(jù) Y 本身得到的預(yù)測(cè)值。

statsmodels 包提供了便捷的檢驗(yàn)方法。它可以接受一個(gè)二維數(shù)組，其中第一列為值，第二列為預(yù)測(cè)因子。

零假設(shè)為：第二列的序列與第一列不存在格蘭杰因果關(guān)系。如果 P 值小于顯著性水平 0.05，就可以拒絕零假設(shè)，從而知道 X 的滯后是起作用的。

第二個(gè)參數(shù) maxlag 表示該檢驗(yàn)中涉及了 Y 的多少個(gè)滯后期。

from statsmodels.tsa.stattools import grangercausalitytestsdf = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv', parse_dates=['date'])df['month'] = df.date.dt.monthgrangercausalitytests(df[['value','month']],maxlag=2)

Granger Causalitynumber of lags (no zero) 1ssr based F test: F=54.7797 , p=0.0000 , df_denom=200, df_num=1ssr based chi2 test: chi2=55.6014 , p=0.0000 , df=1likelihood ratio test: chi2=49.1426 , p=0.0000 , df=1parameter F test: F=54.7797 , p=0.0000 , df_denom=200, df_num=1Granger Causalitynumber of lags (no zero) 2ssr based F test: F=162.6989, p=0.0000 , df_denom=197, df_num=2ssr based chi2 test: chi2=333.6567, p=0.0000 , df=2likelihood ratio test: chi2=196.9956, p=0.0000 , df=2parameterFtest:F=162.6989,p=0.0000,df_denom=197,df_num=2

在上面的例子中，全部測(cè)試結(jié)果中的 P 值都為零，說(shuō)明 'month' 可用作預(yù)測(cè)航班的乘客數(shù)量。