0
| 本文作者： skura | 2019-08-15 18:45 |
無論我們是想預(yù)測金融市場的趨勢還是用電量，時(shí)間都是我們模型中必須考慮的一個(gè)重要因素。例如，預(yù)測一天中什么時(shí)候會(huì)出現(xiàn)用電高峰是很有趣的，可以以此為依據(jù)調(diào)整電價(jià)或發(fā)電量。
輸入時(shí)間序列。時(shí)間序列只是按時(shí)間順序排列的一系列數(shù)據(jù)點(diǎn)。在時(shí)間序列中，時(shí)間往往是獨(dú)立變量，其目標(biāo)通常是預(yù)測未來。
然而，在處理時(shí)間序列時(shí)，還有一些其他因素會(huì)發(fā)揮作用。
它是靜止的嗎？
有季節(jié)性嗎？
目標(biāo)變量是否自相關(guān)？
在這篇文章中，我將介紹時(shí)間序列的不同特征，以及我們?nèi)绾螌?duì)它們進(jìn)行建模才能獲得準(zhǔn)確的預(yù)測。
預(yù)測未來是困難的
自相關(guān)
通俗地說，自相關(guān)是觀測值之間的相似度，它是觀測值之間時(shí)間滯后的函數(shù)。
自相關(guān)示例
上面是一個(gè)自相關(guān)的例子。仔細(xì)觀察，你會(huì)發(fā)現(xiàn)第一個(gè)值和第 24 個(gè)值具有很高的自相關(guān)性。同樣，第 12 個(gè)值和第 36 個(gè)觀測值也高度相關(guān)。這意味著我們將在每 24 個(gè)時(shí)間單位中找到一個(gè)非常相似的值。
注意，這個(gè)圖看起來像正弦函數(shù)。這是季節(jié)性的征兆，你可以通過在上面的圖中找到 24 小時(shí)的周期來找到它的價(jià)值。
季節(jié)性
季節(jié)性是指周期性波動(dòng)。例如，白天的用電量高，晚上的用電量低，或者圣誕節(jié)期間的在線銷售額增加，節(jié)后銷售再次放緩。
季節(jié)性示例
如你所見，每天都有明顯的季節(jié)性。每天晚上，你都會(huì)看到一個(gè)高峰，最低點(diǎn)出現(xiàn)在每天的開始和結(jié)束。
記住，如果季節(jié)性是滿足正弦函數(shù)的，它也可以從自相關(guān)圖中推導(dǎo)出來。簡單地看一下周期，它給出了季節(jié)的長度。
平穩(wěn)性
平穩(wěn)性是時(shí)間序列的一個(gè)重要特征。如果時(shí)間序列的統(tǒng)計(jì)性質(zhì)不隨時(shí)間變化，則稱其為平穩(wěn)的。換句話說，它有不變的均值和方差，協(xié)方差不隨時(shí)間變化。
平穩(wěn)過程示例
再看上面的圖，我們看到上面的過程是平穩(wěn)的，平均值和方差不會(huì)隨時(shí)間變化。
通常，股票價(jià)格不是一個(gè)平穩(wěn)的過程，因?yàn)槲覀兛赡軙?huì)看到一個(gè)增長的趨勢，或者，其波動(dòng)性可能會(huì)隨著時(shí)間的推移而增加（這意味著方差正在變化）。
理想情況下，我們需要一個(gè)用于建模的固定時(shí)間序列。當(dāng)然，不是所有的都是平穩(wěn)的，但是我們可以通過做不同的變換，使它們保持平穩(wěn)。
如何測試過程是否平穩(wěn)
你可能已經(jīng)注意到在上圖的標(biāo)題「Dickey-Fuller」。這是我們用來確定時(shí)間序列是否穩(wěn)定的統(tǒng)計(jì)測試。
在不討論 Dickey-Fuller 測試的技術(shù)特性的情況下，它檢測了單位根是否存在空假設(shè)。
如果是，則 P>0，并且過程不是平穩(wěn)的。
否則，p=0，無效假設(shè)被拒絕，過程被認(rèn)為是平穩(wěn)的。
例如，下面的過程不是平穩(wěn)的。請(qǐng)注意為什么平均值不隨時(shí)間變化。
非平穩(wěn)過程示例
時(shí)間序列建模
有很多方法可以模擬時(shí)間序列來進(jìn)行預(yù)測。在此，我將介紹：
移動(dòng)平均
指數(shù)平滑
ARIMA
移動(dòng)平均
移動(dòng)平均模型可能是最簡單的時(shí)間序列建模方法。這個(gè)模型簡單來說就是，下一個(gè)值是所有過去值的平均值。
雖然很簡單，但是這個(gè)模型的效果可能好到出乎意料，它代表了一個(gè)好的起點(diǎn)。
否則，移動(dòng)平均值可用于識(shí)別數(shù)據(jù)中有趣的趨勢。我們可以定義一個(gè)窗口來應(yīng)用移動(dòng)平均模型來平滑時(shí)間序列，并突出不同的趨勢。
24 小時(shí)窗口上的移動(dòng)平均值示例
在上面的圖中，我們將移動(dòng)平均模型應(yīng)用于一個(gè) 24 小時(shí)窗口。綠線平滑了時(shí)間序列，我們可以看到 24 小時(shí)內(nèi)有 2 個(gè)峰值。
當(dāng)然，窗口越長，趨勢就越平滑。下面是一個(gè)較小窗口上移動(dòng)平均值的示例。
12 小時(shí)窗口上的移動(dòng)平均值示例
指數(shù)平滑
指數(shù)平滑使用與移動(dòng)平均相似的邏輯，但這次，對(duì)每個(gè)觀測值分配了不同的遞減權(quán)重。換言之，離現(xiàn)在的時(shí)間距離越遠(yuǎn)，觀察結(jié)果的重要性就越低。
在數(shù)學(xué)上，指數(shù)平滑表示為：
指數(shù)平滑表達(dá)式
這里，alpha 是一個(gè)平滑因子，它的值介于 0 和 1 之間。它決定了之前觀測值的權(quán)重下降的速度。
指數(shù)平滑示例
在上面的圖中，深藍(lán)色線表示時(shí)間序列的指數(shù)平滑，平滑系數(shù)為 0.3，而橙色線表示平滑系數(shù)為 0.05。
如你所見，平滑因子越小，時(shí)間序列就越平滑。這是有意義的，因?yàn)楫?dāng)平滑因子接近 0 時(shí)，我們接近移動(dòng)平均模型。
雙指數(shù)平滑
當(dāng)時(shí)間序列中存在趨勢時(shí)，使用雙指數(shù)平滑。在這種情況下，我們使用這種技術(shù)，它只是指數(shù)平滑的兩次遞歸使用。
數(shù)學(xué)公式為：
雙指數(shù)平滑表達(dá)式
這里，beta 是趨勢平滑因子，它的值介于 0 和 1 之間。
下面，你可以看到 alpha 和 beta 的不同值如何影響時(shí)間序列的形狀。
雙指數(shù)平滑示例
三指數(shù)平滑
該方法通過添加季節(jié)平滑因子來擴(kuò)展雙指數(shù)平滑。當(dāng)然，如果你注意到時(shí)間序列中的季節(jié)性，這很有用。
在數(shù)學(xué)上，三指數(shù)平滑表示為：
三指數(shù)平滑表達(dá)式
其中 gamma 是季節(jié)平滑因子，L 是季節(jié)長度。
季節(jié)性差分自回歸滑動(dòng)平均模型（SARIMA）
SARIMA 實(shí)際上是簡單模型的組合，可以生成一個(gè)復(fù)雜的模型，該模型可以模擬具有非平穩(wěn)特性和季節(jié)性的時(shí)間序列。
首先，我們得到了自回歸模型 AR(p)。這基本上是時(shí)間序列對(duì)自身的回歸。在這里，我們假設(shè)當(dāng)前值依賴于它以前的值，并且有一定的滯后。它采用一個(gè)表示最大滯后的參數(shù) p。為了找到它，我們查看了部分自相關(guān)圖，在此之后大部分滯后并不顯著。
在下面的例子中，p 的值是 4。
部分自相關(guān)圖示例
然后，我們添加移動(dòng)平均模型 MA(q)。這需要一個(gè)參數(shù) q，它代表自相關(guān)圖上那些滯后不顯著的最大滯后。
下圖中，q 為 4。
自相關(guān)圖示例
之后，我們添加整合順序 I(d)。參數(shù) d 表示使序列平穩(wěn)所需的差異數(shù)。
最后，我們添加最后一部分：季節(jié)性 S(P, D, Q, s)，其中 S 只是季節(jié)的長度。此外，這里要求參數(shù) P 和 Q 與 p 和 q 相同，但用于季節(jié)部分。最后，D 是季節(jié)整合的順序，表示從系列中刪除季節(jié)性所需的差異數(shù)量。
綜合起來，我們得到了 SARIMA(p, d, q)(P, D, Q, s) 模型。
要注意的是：在用 SARIMA 建模之前，我們必須對(duì)時(shí)間序列進(jìn)行轉(zhuǎn)換，以消除季節(jié)性和任何非平穩(wěn)行為。
這是一個(gè)很好的理論！讓我們在第一個(gè)項(xiàng)目中應(yīng)用上面討論的技術(shù)。
我們將想辦法預(yù)測一家公司的股票價(jià)格?，F(xiàn)在，預(yù)測股票價(jià)格幾乎是不可能的。然而，這仍然是一個(gè)有趣的練習(xí)，它將是一個(gè)很好的來實(shí)踐我們所學(xué)到知識(shí)的方法。
項(xiàng)目 1：股票價(jià)格預(yù)測
我們將利用 New Germany Fund（GF）的歷史股價(jià)來預(yù)測未來五個(gè)交易日的收盤價(jià)。
你可以在這里獲取數(shù)據(jù)集和資料。
像往常一樣，我強(qiáng)烈推薦你動(dòng)手編碼！啟動(dòng)你的筆記本，我們開始吧！
你絕不會(huì)因?yàn)檫@個(gè)項(xiàng)目而發(fā)財(cái)
導(dǎo)入數(shù)據(jù)
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()
from sklearn.metrics import r2_score, median_absolute_error, mean_absolute_error
from sklearn.metrics import median_absolute_error, mean_squared_error, mean_squared_log_error
from scipy.optimize import minimize
import statsmodels.tsa.api as smt
import statsmodels.api as sm
from tqdm import tqdm_notebook
from itertools import product
def mean_absolute_percentage_error(y_true, y_pred):
return np.mean(np.abs((y_true - y_pred) / y_true)) * 100
import warnings
warnings.filterwarnings('ignore')
%matplotlib inline
DATAPATH = 'data/stock_prices_sample.csv'
data = pd.read_csv(DATAPATH, index_col=['DATE'], parse_dates=['DATE'])
data.head(10)
首先，我們導(dǎo)入一些庫，這些庫將在整個(gè)分析過程中都會(huì)用到。此外，我們用平均百分比誤差（MAPE）作為我們的誤差度量。
然后，我們導(dǎo)入數(shù)據(jù)集，排在前十的是：
數(shù)據(jù)集的前 10 個(gè)條目
正如你所看到的，我們有一些關(guān)于 New Germany Fund (GF) 不同股票的數(shù)據(jù)。此外，我們還有一個(gè)關(guān)于當(dāng)天信息的數(shù)據(jù)，但我們只需要當(dāng)天結(jié)束（EOD）時(shí)的股票信息。
數(shù)據(jù)清洗
data = data[data.TICKER != 'GEF']
data = data[data.TYPE != 'Intraday']
drop_cols = ['SPLIT_RATIO', 'EX_DIVIDEND', 'ADJ_FACTOR', 'ADJ_VOLUME', 'ADJ_CLOSE', 'ADJ_LOW', 'ADJ_HIGH', 'ADJ_OPEN', 'VOLUME', 'FREQUENCY', 'TYPE', 'FIGI']
data.drop(drop_cols, axis=1, inplace=True)
data.head()
首先，我們刪除不需要的條目。
然后，我們刪除不需要的列，因?yàn)槲覀冎幌腙P(guān)注股票的收盤價(jià)。
如果預(yù)覽數(shù)據(jù)集，則應(yīng)該看到的是：
清洗后的數(shù)據(jù)集
令人驚嘆！我們準(zhǔn)備好進(jìn)行探索性數(shù)據(jù)分析了！
探索性數(shù)據(jù)分析（EDA）
# Plot closing price
plt.figure(figsize=(17, 8))
plt.plot(data.CLOSE)
plt.title('Closing price of New Germany Fund Inc (GF)')
plt.ylabel('Closing price ($)')
plt.xlabel('Trading day')
plt.grid(False)
plt.show()
我們繪制數(shù)據(jù)集整個(gè)時(shí)間段的收盤價(jià)。
你將會(huì)得到：
New Germany Fund （GF）收盤價(jià)
很明顯，你看到的不是一個(gè)平穩(wěn)的過程，很難判斷是否有某種季節(jié)性。
移動(dòng)平均
讓我們使用移動(dòng)平均模型來平滑我們的時(shí)間序列。為此，我們將使用一個(gè)輔助函數(shù)，該函數(shù)將在指定的時(shí)間窗口上運(yùn)行移動(dòng)平均模型，并繪制結(jié)果平滑曲線：
def plot_moving_average(series, window, plot_intervals=False, scale=1.96):
rolling_mean = series.rolling(window=window).mean()
plt.figure(figsize=(17,8))
plt.title('Moving average\n window size = {}'.format(window))
plt.plot(rolling_mean, 'g', label='Rolling mean trend')
#Plot confidence intervals for smoothed values
if plot_intervals:
mae = mean_absolute_error(series[window:], rolling_mean[window:])
deviation = np.std(series[window:] - rolling_mean[window:])
lower_bound = rolling_mean - (mae + scale * deviation)
upper_bound = rolling_mean + (mae + scale * deviation)
plt.plot(upper_bound, 'r--', label='Upper bound / Lower bound')
plt.plot(lower_bound, 'r--')
plt.plot(series[window:], label='Actual values')
plt.legend(loc='best')
plt.grid(True)
#Smooth by the previous 5 days (by week)
plot_moving_average(data.CLOSE, 5)
#Smooth by the previous month (30 days)
plot_moving_average(data.CLOSE, 30)
#Smooth by previous quarter (90 days)
plot_moving_average(data.CLOSE, 90, plot_intervals=True)
使用5天的時(shí)間窗口，我們得到：
上一個(gè)交易周的平滑曲線
如你所見，我們幾乎看不到趨勢，因?yàn)樗咏鼘?shí)際曲線。讓我們看看上個(gè)月和上個(gè)季度的平滑處理結(jié)果。
上個(gè)月（30 天前）的平滑曲線
按上一季度（90 天）平滑
現(xiàn)在更容易發(fā)現(xiàn)趨勢。請(qǐng)注意，30 天和 90 天的趨勢圖在末尾顯示一條向下的曲線。這意味著股票可能在接下來的幾天內(nèi)會(huì)下跌。
指數(shù)平滑
現(xiàn)在，讓我們用指數(shù)平滑來看看它是否能獲得更好的趨勢。
def exponential_smoothing(series, alpha):
result = [series[0]] # first value is same as series
for n in range(1, len(series)):
result.append(alpha * series[n] + (1 - alpha) * result[n-1])
return result
def plot_exponential_smoothing(series, alphas):
plt.figure(figsize=(17, 8))
for alpha in alphas:
plt.plot(exponential_smoothing(series, alpha), label="Alpha {}".format(alpha))
plt.plot(series.values, "c", label = "Actual")
plt.legend(loc="best")
plt.axis('tight')
plt.title("Exponential Smoothing")
plt.grid(True);
plot_exponential_smoothing(data.CLOSE, [0.05, 0.3])
這里，我們使用 0.05 和 0.3 作為平滑因子的值。當(dāng)然你也可以嘗試其他值，看看結(jié)果如何。
指數(shù)平滑
如您所見，alpha 值 0.05 平滑了曲線，同時(shí)剔除了大部分向上和向下的趨勢。
現(xiàn)在，讓我們使用雙指數(shù)平滑。
雙指數(shù)平滑
def double_exponential_smoothing(series, alpha, beta):
result = [series[0]]
for n in range(1, len(series)+1):
if n == 1:
level, trend = series[0], series[1] - series[0]
if n >= len(series): # forecasting
value = result[-1]
else:
value = series[n]
last_level, level = level, alpha * value + (1 - alpha) * (level + trend)
trend = beta * (level - last_level) + (1 - beta) * trend
result.append(level + trend)
return result
def plot_double_exponential_smoothing(series, alphas, betas):
plt.figure(figsize=(17, 8))
for alpha in alphas:
for beta in betas:
plt.plot(double_exponential_smoothing(series, alpha, beta), label="Alpha {}, beta {}".format(alpha, beta))
plt.plot(series.values, label = "Actual")
plt.legend(loc="best")
plt.axis('tight')
plt.title("Double Exponential Smoothing")
plt.grid(True)
plot_double_exponential_smoothing(data.CLOSE, alphas=[0.9, 0.02], betas=[0.9, 0.02])
你將得到：
雙指數(shù)平滑
同樣，用不同的 α 和 β 組合進(jìn)行實(shí)驗(yàn)，以獲得更好的曲線。
建模
如前所述，我們必須將序列轉(zhuǎn)換為一個(gè)平穩(wěn)的過程，以便對(duì)其進(jìn)行建模。因此，讓我們應(yīng)用 Dickey-Fuller 測試來看看它是否是一個(gè)平穩(wěn)的過程：
def tsplot(y, lags=None, figsize=(12, 7), syle='bmh'):
if not isinstance(y, pd.Series):
y = pd.Series(y)
with plt.style.context(style='bmh'):
fig = plt.figure(figsize=figsize)
layout = (2,2)
ts_ax = plt.subplot2grid(layout, (0,0), colspan=2)
acf_ax = plt.subplot2grid(layout, (1,0))
pacf_ax = plt.subplot2grid(layout, (1,1))
y.plot(ax=ts_ax)
p_value = sm.tsa.stattools.adfuller(y)[1]
ts_ax.set_title('Time Series Analysis Plots\n Dickey-Fuller: p={0:.5f}'.format(p_value))
smt.graphics.plot_acf(y, lags=lags, ax=acf_ax)
smt.graphics.plot_pacf(y, lags=lags, ax=pacf_ax)
plt.tight_layout()
tsplot(data.CLOSE, lags=30)
# Take the first difference to remove to make the process stationary
data_diff = data.CLOSE - data.CLOSE.shift(1)
tsplot(data_diff[1:], lags=30)
你將看到：
通過 DickeyFuller 測試，時(shí)間序列是非平穩(wěn)的。另外，從自相關(guān)圖來看，我們發(fā)現(xiàn)它似乎沒有明顯的季節(jié)性。
因此，為了消除高度自相關(guān)并使過程穩(wěn)定，讓我們?nèi)〉谝粋€(gè)差異（代碼塊中的第 23 行）。我們簡單地用一天的滯后時(shí)間減去時(shí)間序列，得到：
令人驚嘆的！我們的序列現(xiàn)在是平穩(wěn)的，可以開始建模了！
SARIMA
#Set initial values and some bounds
ps = range(0, 5)
d = 1
qs = range(0, 5)
Ps = range(0, 5)
D = 1
Qs = range(0, 5)
s = 5
#Create a list with all possible combinations of parameters
parameters = product(ps, qs, Ps, Qs)
parameters_list = list(parameters)
len(parameters_list)
# Train many SARIMA models to find the best set of parameters
def optimize_SARIMA(parameters_list, d, D, s):
"""
Return dataframe with parameters and corresponding AIC
parameters_list - list with (p, q, P, Q) tuples
d - integration order
D - seasonal integration order
s - length of season
"""
results = []
best_aic = float('inf')
for param in tqdm_notebook(parameters_list):
try: model = sm.tsa.statespace.SARIMAX(data.CLOSE, order=(param[0], d, param[1]),
seasonal_order=(param[2], D, param[3], s)).fit(disp=-1)
except:
continue
aic = model.aic
#Save best model, AIC and parameters
if aic < best_aic:
best_model = model
best_aic = aic
best_param = param
results.append([param, model.aic])
result_table = pd.DataFrame(results)
result_table.columns = ['parameters', 'aic']
#Sort in ascending order, lower AIC is better
result_table = result_table.sort_values(by='aic', ascending=True).reset_index(drop=True)
return result_table
result_table = optimize_SARIMA(parameters_list, d, D, s)
#Set parameters that give the lowest AIC (Akaike Information Criteria)
p, q, P, Q = result_table.parameters[0]
best_model = sm.tsa.statespace.SARIMAX(data.CLOSE, order=(p, d, q),
seasonal_order=(P, D, Q, s)).fit(disp=-1)
print(best_model.summary())
現(xiàn)在，對(duì)于 SARIMA，我們首先定義一些參數(shù)值的范圍，以生成 p, q, d, P, Q, D, s 的所有可能組合的列表。
現(xiàn)在，在上面的代碼單元中，我們有 625 種不同的組合！我們將嘗試每種組合，并訓(xùn)練 SARIMA，以便找到性能最佳的模型。這可能需要一些時(shí)間，具體多長時(shí)間取決于計(jì)算機(jī)的處理能力。
完成后，我們將輸出最佳模型的摘要，你將看到：
令人驚嘆！最后，我們預(yù)測未來五個(gè)交易日的收盤價(jià)，并評(píng)估模型的 MAPE。
在這種情況下，有一個(gè) 0.79% 的 MAPE，這是非常好的！
將預(yù)測價(jià)格與實(shí)際數(shù)據(jù)進(jìn)行比較
# Make a dataframe containing actual and predicted prices
comparison = pd.DataFrame({'actual': [18.93, 19.23, 19.08, 19.17, 19.11, 19.12],
'predicted': [18.96, 18.97, 18.96, 18.92, 18.94, 18.92]},
index = pd.date_range(start='2018-06-05', periods=6,))
#Plot predicted vs actual price
plt.figure(figsize=(17, 8))
plt.plot(comparison.actual)
plt.plot(comparison.predicted)
plt.title('Predicted closing price of New Germany Fund Inc (GF)')
plt.ylabel('Closing price ($)')
plt.xlabel('Trading day')
plt.legend(loc='best')
plt.grid(False)
plt.show()
現(xiàn)在，為了將我們的預(yù)測與實(shí)際數(shù)據(jù)進(jìn)行比較，我們從雅虎財(cái)務(wù)（YahooFinance）獲取財(cái)務(wù)數(shù)據(jù)并創(chuàng)建一個(gè)數(shù)據(jù)框架。
然后，我們繪出曲線，看看我們與實(shí)際收盤價(jià)的差距有多大：
預(yù)計(jì)值和實(shí)際收盤價(jià)比較
我們的預(yù)測似乎有點(diǎn)偏離。事實(shí)上，預(yù)測價(jià)格很平穩(wěn)，這意味著我們的模型可能表現(xiàn)不佳。
當(dāng)然，這不是因?yàn)槲覀兊某绦?，而是因?yàn)轭A(yù)測股票價(jià)格基本上是不可能的。
從第一個(gè)項(xiàng)目開始，我們學(xué)習(xí)了在使用 SARIMA 建模之前平滑時(shí)間序列的整個(gè)過程。
現(xiàn)在，讓我們介紹一下 Facebook 的 Prophet。它是一個(gè)在 python 和 r 中都可用的預(yù)測工具。該工具幫助生成高質(zhì)量的預(yù)測。
讓我們看看如何在第二個(gè)項(xiàng)目中使用它！
項(xiàng)目2-使用 Prophet 預(yù)測空氣質(zhì)量
標(biāo)題說明了一切：我們將使用 Prophet 來幫助我們預(yù)測空氣質(zhì)量！
導(dǎo)入數(shù)據(jù)
import warnings
warnings.filterwarnings('ignore')
import numpy as np
import pandas as pd
from scipy import stats
import statsmodels.api as sm
import matplotlib.pyplot as plt
%matplotlib inline
DATAPATH = 'data/AirQualityUCI.csv'
data = pd.read_csv(DATAPATH, sep=';')
data.head()
打印出前五行：
如你所見，數(shù)據(jù)集包含有關(guān)不同氣體濃度的信息。每天隔一個(gè)小時(shí)記錄一次。
如果更深入地研究數(shù)據(jù)集，會(huì)發(fā)現(xiàn)有許多值 -200 的實(shí)例。當(dāng)然，負(fù)濃度是沒有意義的，所以我們需要在建模前清洗數(shù)據(jù)。
數(shù)據(jù)清洗與特征工程
# Make dates actual dates
data['Date'] = pd.to_datetime(data['Date'])
# Convert measurements to floats
for col in data.iloc[:,2:].columns:
if data[col].dtypes == object:
data[col] = data[col].str.replace(',', '.').astype('float')
# Compute the average considering only the positive values
def positive_average(num):
return num[num > -200].mean()
# Aggregate data
daily_data = data.drop('Time', axis=1).groupby('Date').apply(positive_average)
# Drop columns with more than 8 NaN
daily_data = daily_data.iloc[:,(daily_data.isna().sum() <= 8).values]
# Remove rows containing NaN values
daily_data = daily_data.dropna()
# Aggregate data by week
weekly_data = daily_data.resample('W').mean()
# Plot the weekly concentration of each gas
def plot_data(col):
plt.figure(figsize=(17, 8))
plt.plot(weekly_data[col])
plt.xlabel('Time')
plt.ylabel(col)
plt.grid(False)
plt.show()
for col in weekly_data.columns:
plot_data(col)
在這里，我們首先分析日期列，將其轉(zhuǎn)換為日期類型。
然后，我們把所有的測量值轉(zhuǎn)換成浮點(diǎn)數(shù)。
之后，我們用每天的平均值來匯總數(shù)據(jù)。
我們還有一些需要?jiǎng)h除的 NAN。
最后，我們按周匯總數(shù)據(jù)，這將提供一個(gè)更平滑的分析趨勢。
我們可以畫出每種化學(xué)物質(zhì)濃度的趨勢。這里，我們展示了 NOx。
NOx 濃度
氮氧化物是非常有害的，因?yàn)樗鼈儠?huì)形成煙霧和酸雨，同時(shí)也會(huì)形成細(xì)顆粒和臭氧。這些都會(huì)對(duì)健康產(chǎn)生不利影響，因此氮氧化物的濃度是空氣質(zhì)量的一個(gè)關(guān)鍵特征。
建模
# Drop irrelevant columns
cols_to_drop = ['PT08.S1(CO)', 'C6H6(GT)', 'PT08.S2(NMHC)', 'PT08.S4(NO2)', 'PT08.S5(O3)', 'T', 'RH', 'AH']
weekly_data = weekly_data.drop(cols_to_drop, axis=1)
# Import Prophet
from fbprophet import Prophet
import logging
logging.getLogger().setLevel(logging.ERROR)
# Change the column names according to Prophet's guidelines
df = weekly_data.reset_index()
df.columns = ['ds', 'y']
df.head()
# Split into a train/test set
prediction_size = 30
train_df = df[:-prediction_size]
# Initialize and train a model
m = Prophet()
m.fit(train_df)
# Make predictions
future = m.make_future_dataframe(periods=prediction_size)
forecast = m.predict(future)
forecast.head()
# Plot forecast
m.plot(forecast)
# Plot forecast's components
m.plot_components(forecast)
# Evaluate the model
def make_comparison_dataframe(historical, forecast):
return forecast.set_index('ds')[['yhat', 'yhat_lower', 'yhat_upper']].join(historical.set_index('ds'))
cmp_df = make_comparison_dataframe(df, forecast)
cmp_df.head()
def calculate_forecast_errors(df, prediction_size):
df = df.copy()
df['e'] = df['y'] - df['yhat']
df['p'] = 100 * df['e'] / df['y']
predicted_part = df[-prediction_size:]
error_mean = lambda error_name: np.mean(np.abs(predicted_part[error_name]))
return {'MAPE': error_mean('p'), 'MAE': error_mean('e')}
for err_name, err_value in calculate_forecast_errors(cmp_df, prediction_size).items():
print(err_name, err_value)
# Plot forecast with upper and lower bounds
plt.figure(figsize=(17, 8))
plt.plot(cmp_df['yhat'])
plt.plot(cmp_df['yhat_lower'])
plt.plot(cmp_df['yhat_upper'])
plt.plot(cmp_df['y'])
plt.xlabel('Time')
plt.ylabel('Average Weekly NOx Concentration')
plt.grid(False)
plt.show()
我們將只關(guān)注氮氧化物濃度。因此，我們刪除所有其他不相關(guān)的列。
然后，我們導(dǎo)入 Prophet。
Prophet 要求日期列命名為 ds，特征列命名為 y，因此我們進(jìn)行了適當(dāng)?shù)母摹?/p>
此時(shí)，我們的數(shù)據(jù)如下：
然后，我們定義一個(gè)訓(xùn)練集。為此，我們將保留最后 30 個(gè)條目進(jìn)行預(yù)測和驗(yàn)證。
之后，我們簡單地初始化 Prophet，將模型與數(shù)據(jù)匹配，并進(jìn)行預(yù)測！
你會(huì)看到：
這里，yhat 代表預(yù)測值，yhat_lower 和 yhat_upper 分別代表預(yù)測值的下限和上限。
Prophet 讓你可以輕松繪制預(yù)測圖，我們得到：
NOx 濃度預(yù)測
如你所見，Prophet 只是用一條直線來預(yù)測未來的 NOx 濃度。
然后，我們檢查時(shí)間序列是否具有某些有趣的特性，例如季節(jié)性：
在這里，Prophet 沒有發(fā)現(xiàn)季節(jié)性的趨勢。
通過計(jì)算模型的平均絕對(duì)百分誤差（MAPE）和平均絕對(duì)誤差（MAE）來評(píng)估模型的性能，我們發(fā)現(xiàn) MAPE 為 13.86%，MAE 為 109.32，這還不錯(cuò)！記住，我們根本沒有對(duì)模型進(jìn)行微調(diào)。
最后，我們只需繪制預(yù)測的上限和下限：
每周 NOx 平均濃度預(yù)測
恭喜你達(dá)到目的！這篇文章很長，但內(nèi)容豐富。你學(xué)會(huì)了如何強(qiáng)有力地分析和建模時(shí)間序列，并將你的知識(shí)應(yīng)用到兩個(gè)不同的項(xiàng)目中。我希望你覺得這篇文章有用。
雷鋒網(wǎng)雷鋒網(wǎng)雷鋒網(wǎng)
雷峰網(wǎng)版權(quán)文章，未經(jīng)授權(quán)禁止轉(zhuǎn)載。詳情見轉(zhuǎn)載須知。
