python時間序列科學計算基本知識(1.5萬字 19圖 閱讀需要37分鐘)

yuanzhoulvpi 2021-09-19 23:15:35 阅读数:752

python 序列 基本 需要

介紹

代碼、數據全部免費,都放在我的gitee倉庫裏面https://gitee.com/yuanzhoulvpi/time_series,想要使用的可以直接到我這個倉庫下載。

本文是繼python與時間序列(開篇)的第一篇,詳細介紹了一些時間序列的描述性統計知識,並且附有代碼、數據。百分百可以複現。

本文內容比較多【初略統計:所有文本超過1.5萬字;代碼超過600行】。

具體的jupyter notebook結構如下:

  1. 什麼是時間序列。
  2. 在python裏面如何導入時間序列數據。
  3. 什麼是面板數據。
  4. 可視化時間序列(包括時間序列的區域填充圖,季節性時間序列圖,箱線圖等)。
  5. 時間序列的幾種模式以及模式的分解。
  6. 時間序列的平穩性的介紹、原理、解釋等;以及讓時間序列平穩的方法。
  7. 白噪聲和平穩時間序列的區別。
  8. 時間序列的趨勢性、季節性如何提取、去除。
  9. 如何處理帶有缺失值的時間序列數據。
  10. 時間序列的自相關、偏自相關數據。
  11. 時間序列的滯後圖。
  12. 時間序列如何平滑。

1. 什麼是時間序列

同一統計指標數值按照時間先後順序排列而成的數據。本質上是反映一個變量隨時間序列變化的趨勢。

  1. 簡單的例子就像是學生每一年的身高數據,按照時間順序記錄下來的數據也是一個時間序列。變量是我們的身高;
  2. 我們支付寶或者微信的零錢,每一天每一個月都有一個實際的值。按照時間順序記錄下來,也是一個時間序列(沒做過,每一次打開支付寶都是在流淚)。

2. 在python裏面如何導入時間序列數據

一般python導入csv、excel格式的文件都是可以使用pandas包導入,如果時間格式需要被處理,可以使用datetime包來處理。

比如這樣:

import numpy as np
import pandas as pd
import datetime

df = pd.read_csv('../datasets/a10.csv')
df['date'] = df['date'].apply(lambda x: datetime.datetime.strptime(x, "%Y-%m-%d"))
df.head()

輸出的結果如下:

3. 什麼是面板數據

有多個時間節點,且每一個時間節點上有多個維度的數據。

舉一個簡單的例子,如果你是一個學生,每天都要考試,2021-01-01日你語文90、數學100、英語102;2021-01-02日你語文10、數學110、英語122 ……這樣的就是就是一個面板數據,有多個時間節點,且每一個時間節點,有多個維度(語文、數學、英語);

舉一個現實的例子;後來這個學生長大了,每天有自己的經濟來源了,他還是喜歡記錄數據,他會每天都會這麼記錄:2021-01-01日給女友花費1000元、兼職20元、日常開銷200元;2021-01-02日給女友10元、兼職200元、日常開銷10元;…… 這也是一個面板數據。(怎麼有點心酸呢)。

4. 可視化時間序列

  1. 最常見的就是時序圖:
import matplotlib.pyplot as plt

df = pd.read_csv('../datasets/a10.csv')
df['date'] = df['date'].apply(lambda x: datetime.datetime.strptime(x, "%Y-%m-%d"))

def plot_df(df, x, y, title="", xlabel="Date", ylabel="Value", dpi=100):
    plt.figure(figsize=(165), dpi=dpi)
    plt.plot(df[x], df[y], color='tab:red')
    plt.gca().set(title=title, xlabel=xlabel, ylabel=ylabel)
    plt.show()
plot_df(df, x="date", y="value")
  1. 區域填充圖:

fig, ax = plt.subplots(11, figsize=(125), dpi=100)
ax.fill_between(x=df['date'], y1=df['value'], y2=-df['value'], alpha=0.5, linewidth=2, color='seagreen')
ax.set_title("demo data", fontsize=16)
ax.hlines(y=0, xmin=np.min(df['date']), xmax=np.max(df['date']), linewidth=0.5)
plt.show()
  1. 季節性圖
df['year'] = df['date'].apply(lambda x: x.year)
df['month'] = df['date'].apply(lambda x: int(x.strftime('%m')))

fig, ax = plt.subplots(figsize=(106))
for index, tempyear in enumerate(df['year'].unique().tolist()):
    tempdata = df.loc[df['year'] == tempyear].sort_values(by=['date'])
    ax.plot(tempdata['month'], tempdata['value'], label=tempyear)
    ax.text(tempdata.tail(1)['month'], tempdata.tail(1)['value'], tempyear)

font = {'family''serif',
        'color''darkred',
        'weight''normal',
        'size'16,
        }
ax.set_ylabel("$Drug Sales$", fontdict=font)
ax.set_xlabel('$Month$', fontdict=font)
ax.set_title("Seasibak Plot of Drug Sales Time Series",
             fontdict=font)
x_axis_ticks = np.arange(start=1, stop=13)
ax.set_xticks(x_axis_ticks)
ax.set_xticklabels(x_axis_ticks)
  1. 對季節性做下鑽,可以使用boxplot
fig, ax = plt.subplots(ncols=2, nrows=1, figsize=(165))
sns.boxplot(x='year', y='value', data=df, ax=ax[0])
sns.boxplot(x='month', y='value', data=df, ax=ax[1])

ax[0].set_title("Year-wise BoxPlot \n(The Trend)", fontsize=18)
xlabel = ax[0].get_xticklabels()
ax[0].set_xticklabels(xlabel, rotation=45)
ax[1].set_title("Month-wise BoxPlot \n(The Trend)", fontsize=18)

5. 時間序列中的模式

任何時間序列可以可以被拆分為3個部分:

  1. 趨勢:趨勢是比較明顯的,比如極速的上昇或者迅速下跌。
  2. 季節性:可以在數據中看到明顯的周期性,並且這個周期性和時間周期有關。這個周期可能是月,可能是季度,也可能是年。
  3. 誤差項。

但是不是說所有的時間序列都必須具備這3個部分。時間序列可能沒有明顯的趨勢、可能沒有明顯的周期性。或者兩個都沒有。 因此,可以將時間序列看成趨勢、周期性、誤差項的組合。

  1. 幾種模式的對比
fig, ax = plt.subplots(nrows=1, ncols=3, figsize=(204), dpi=100)
pd.read_csv("../datasets/guinearice.csv", parse_dates=['date'], index_col=['date']).plot(title="Trend Only", legend=False,
                                                                                      ax=ax[0])
pd.read_csv("../datasets/sunspotarea.csv", parse_dates=['date'], index_col=['date']).plot(title="Seasonality Only",
                                                                                       legend=False, ax=ax[1])
pd.read_csv("../datasets/AirPassengers.csv", parse_dates=['date'], index_col=['date']).plot(title="Trend and Seasonality",
                                                                                         legend=False, ax=ax[2])

5.1 時間序列中的模式分解

基於原始的時間序列的趨勢項和季節項,時間序列模型可以被分為加法模型和乘法模型.

  1. 乘法型:時間序列值 = 趨勢項 * 季節項 * 誤差項。
  2. 加法型:時間序列值 = 趨勢項 + 季節項 + 誤差項。

對時間序列的分解,可以從趨勢項、季節項、 誤差項的乘法和加法的角度進行。 代碼方面,我們使用statsmodels包的seasonal_decompose函數。

from statsmodels.tsa.seasonal import seasonal_decompose
from dateutil.parser import parse

# 導入數據
df = pd.read_csv("../datasets/a10.csv", parse_dates=['date'], index_col=['date'])

# 乘法模型
result_mul = seasonal_decompose(df['value'], model="multilicative", extrapolate_trend='freq')
result_add = seasonal_decompose(df['value'], model="additive", extrapolate_trend='freq')

# 畫圖
fig, ax = plt.subplots(ncols=2, nrows=4, figsize=(147), sharex=True)


def plot_decompose(result, ax, index, title, fontdict=font):
    ax[0, index].set_title(title, fontdict=fontdict)
    result.observed.plot(ax=ax[0, index])
    ax[0, index].set_ylabel("Observed")

    result.trend.plot(ax=ax[1, index])
    ax[1, index].set_ylabel("Trend")

    result.seasonal.plot(ax=ax[2, index])
    ax[2, index].set_ylabel("Seasonal")

    result.resid.plot(ax=ax[3, index])
    ax[3, index].set_ylabel("Resid")


plot_decompose(result=result_add, ax=ax, index=0, title="Additive Decompose", fontdict=font)
plot_decompose(result=result_mul, ax=ax, index=1, title="Multiplicative Decompose", fontdict=font)

  1. 在上面的代碼中,設置 extrapolate_trend='freq'是為了填充 季節項、誤差項中開頭的缺失值。
  2. 在上圖中,可以發現加法模型裏面的誤差項還有一部分時間序列的模式沒有被提取完。但是乘法模型的誤差項(或者叫殘差項)基本上看不出來任何信息了,說明乘法模型對這個數據來說,有更强提取數據信息的能力。
  3. 乘法模型的季節項、趨勢項、誤差項數據都保存在 result_mul裏面。我們將這幾個數據提取出來。
df_reconstructed = pd.concat([result_mul.seasonal, result_mul.trend, result_mul.resid, result_mul.observed], axis=1)
df_reconstructed.columns = ['Seasonal''Trend''Resid''Actual_value']
df_reconstructed.head()

可以檢查一下上面的數據;基本上可以確定df_reconstructed['Seasonal'] * df_reconstructed['Trend'] * df_reconstructed['Resid'] = df_reconstructed['Actual_value'], 或者我們用均方誤差算一下:


value = np.sqrt(np.sum((df_reconstructed['Seasonal'] * df_reconstructed['Trend'] * df_reconstructed['Resid'] -
                        df_reconstructed['Actual_value']) ** 2))
value

6. 時間序列的平穩性

  1. 平穩性是時間序列的一個屬性,一個平穩的時間序列指的是這個時間序列和時間無關,也就是說這個時間序列不是時間的一個函數。
  2. 也就是說一個時間序列是平穩的,就是這個時間序列的幾個統計量:均值、方差、自相關系數都是一個常數,和時間無關。

6.1 如何讓一個時間序列平穩

  1. 對一個時間序列做一次或者多次差分。
  2. 季節性延遲(或者叫多階差分)。
  3. 取時間序列的第N個根。
  4. 上面三種方法混合使用。

最方便最常用的平穩一個時間序列,就是多次使用一階差分,直到時間序列平穩。

6.2在做預測之前為什麼需要對非平穩時間序列做平穩

  1. 預測一個平穩的時間序列相對容易,且結果更加可靠。
  2. 一個更加重要的原因是:自回歸預測模型本質上是一個線性回歸模型,它利用時間序列的滯後項做為預測變量。 我們都知道,在線性回歸裏面,如果預測變量的X都是互不相關的,那麼線性回歸預測的效果最好。因此對序列進行平穩化就解决了變量之間的相關性問題, 從而消除來時間序列的自相關,使得預測模型中的預測變量幾乎相互獨立。

現在已經知道的時間序列的平穩性的重要性,如何判斷一個時間序列是否平穩呢?

6.3檢驗時間序列是否平穩

檢驗時間序列是否為平穩性時間序列可以有這幾種方法:

  1. 查看時間序列的時序圖;
  2. 將時間序列分為2個或者多個連續片段,然後計算對應的統計量(比如均值、方差、自相關系數等),如果片段的統計量都是明顯不相等的話,說明肯定不是平穩性。
  3. 不管怎麼樣,需要使用一個量化的指標去衡量一個時間序列是否平穩,可以使用單比特根檢驗方法來對時間序列做平穩性檢驗,檢驗時間序列是否平穩並且具有單比特根。

還有很多單比特根檢驗方法的增强版本:

  1. Augmented Dickey Fuller test(ADH Test)
  2. Kwiatkowski-Phillips-Schmidt-Shin – KPSS test (trend stationary)
  3. Philips Perron test (PP Test)

最廣泛使用的方法是ADF test;原假設是時間序列有單比特根並且是非平穩的。如果ADF test的P-value小於顯著水平(0.05),就可以拒絕原假設。 KPSS testADH test恰恰相反,是為了證明時間序列有一個確定的趨勢。

from statsmodels.tsa.stattools import adfuller, kpss

df = pd.read_csv("../datasets/a10.csv", parse_dates=['date'])

# ADF test
result = adfuller(df.value.values, autolag='AIC')

print('*' * 50)
print(f"ADF Statistic: {result[0]}; p-value: {result[1]}")

for key, value in result[4].items():
    print('Crital Values:')
    print(f"{key}{value}")

# KPSS Test

result = kpss(df.value.values, regression='c', nlags='auto')

print('*' * 50)
print(f'\nKPSS Stattistic: {result[0]:.10f}; p-value: {result[1]:.3f}')

for key, value in result[3].items():
    print('Crital Values:')
    print(f"{key}{value}")

結果:

6.4白噪聲和平穩時間序列的區別

和平穩性時間序列一樣,白噪聲時間序列的每一個值和時間也是無關的,也就是說白噪聲的均值和方差不隨著時間改變。唯一一個不同點就是白噪聲是完全隨機的、且均值為0的時間序列。

白噪聲中沒有任何時間序列模式,生活中遇到的白噪聲就像是FM收音機中空白聲音。數學上的白噪聲就是一串均值為0的完全隨機的隨機數。

fig, ax = plt.subplots(figsize=(104))
randval = np.random.randn(1000)
pd.Series(randval).plot(title="Random White Noise", color='k', ax=ax)

8.時間序列的去趨勢

去趨勢時間序列指的是從時間序列中去趨勢成分。去除時間序列趨勢有多種方法。

  1. 從時間序列中獲得最佳擬合線。可以使用線性回歸來擬合時間序列,有時候可能還用得到二次項(x^2);
  2. 將剛才我們提取的時間序列成分中趨勢項部分移除;
  3. 使用像是Baxter-King過濾器或者Hodrick-Prescott過濾器去除移動平均趨勢線或者周期性成分。

8.1去時間序列的趨勢項趨勢

# use scipy: Subtract the Line of best fit
from scipy import signal

df = pd.read_csv("../datasets/a10.csv", parse_dates=['date'])
detrend = signal.detrend(df.value.values)

fig, ax = plt.subplots(ncols=2, nrows=1, figsize=(154))
ax[0].plot(detrend)
ax[0].set_title("Drug Sales detrended by subtracting the least squares fit")

# Using  statmodels: Subtracting the Trend Component
from statsmodels.tsa.seasonal import seasonal_decompose

df = pd.read_csv('../datasets/a10.csv', parse_dates=['date'], index_col='date')
result_mul = seasonal_decompose(df['value'], model="multilicative", extrapolate_trend='freq')
detrended = df.value.values / result_mul.trend
ax[1].plot(detrended)
ax[1].set_title("Drug Sales detrended by subtracting the trend component")

8.2去時間序列的季節性趨勢

fig, ax = plt.subplots(figsize=(84))

# Subtracing the Trend Component
from statsmodels.tsa.seasonal import seasonal_decompose

df = pd.read_csv('../datasets/a10.csv', parse_dates=['date'], index_col='date')
result_mul = seasonal_decompose(df['value'], model="multilicative", extrapolate_trend='freq')
detrended = df.value.values / result_mul.seasonal
ax.plot(detrended)
ax.set_title("Drug Sales Deseasonlized")

9.檢驗時間序列的季節性

最常用的方式是將時間序列的圖畫出來,然後從圖中查看重複的模式和固定時間間隔,一般周期性都是由日期或者時間决定:

  1. 一天中的小時;
  2. 一個月中的每天;
  3. 每周的;
  4. 每月的;
  5. 每年的。

如果想對季節性做更加准確的季節性檢驗,就使用自相關函數圖(ACF plot);但是當存在相當明顯的季節性模式的時候,ACF圖會在季節性窗口長度的倍數處出現明顯的峰值。

比如,前面的毒品銷售是每月的系列,每年都有重複的模式,比如在下面的第12,第24,第36比特置,都可以看到線有個小峰值。但是必須注意:在實際數據集中,這些强累的模式幾乎不會被注意到,而且可能還會被噪音幹擾。

from pandas.plotting import autocorrelation_plot

df = pd.read_csv("../datasets/a10.csv")
fig, ax = plt.subplots(figsize=(104))
autocorrelation_plot(df.value, ax=ax)

ybound = ax.get_ybound()
for i in 12 * np.arange(start=1, stop=13):
    ax.vlines(x=i, ymin=ybound[0], ymax=ybound[1], colors='red', linestyles='--', linewidth=0.6)
    ax.text(x=i, y=ybound[1], s=f"{i:d}th", fontdict={'size'8})

ax.grid(False)

相應的,如果需要統計檢驗,可以使用CH-test來决定是否需要對時間序列做季節差分。

10.對帶有缺失值的時間序列處理

現實中的時間序列數據可能會缺少時間(比如少了幾天,少了幾年),這也就意味著,這些時間段內的數據沒有被捕捉到或者不能用,在這些缺失的時間段範圍內,測量值可能為0,這樣的情况下,可以使0來填充這個時間段內的數據。 對時間序列做填充的時候,不能簡單的對缺失值用均值填充,應該使用向前填充,用後一個數據去填前面的數據。 處理現實的數據,需要做很多嘗試,來確定哪一個填充效果最好,比如有:

  1. 向前\向後填充;
  2. 線性插值、二次插值、最近鄰插值、使用季節性插值。

為了衡量插值的效果,做了一次小實驗,手動對數據做缺失值處理,然後使用不同方法做插值。

# 生成隨機數
from scipy.interpolate import interp1d
from sklearn.metrics import mean_squared_error

df_orig = pd.read_csv("../datasets/a10.csv", parse_dates=['date'], index_col=['date'])

# 制作缺失值數據
np.random.seed(2021)
df = df_orig.copy()
df.iloc[np.random.randint(low=0, high=df_orig.shape[0], size=30)] = None

# 畫圖
fig, ax = plt.subplots(nrows=7, ncols=1, figsize=(1012), sharex=True)

# 1. 原始數據 -------------------------------
df_orig.plot(title="Acutal", ax=ax[0], label='Acutal', color='red', style='.-')
df.plot(ax=ax[0], label='Acutal', color='green', style='.-')
ax[0].legend(['Missing Data''Available Data'])

# 2. 向前填充 -------------------------------
df_ffill = df.ffill()
error = mean_squared_error(df_orig['value'], df_ffill['value'])
df_ffill['value'].plot(title=f"Forward Fill (MSE: {error:.3f})", ax=ax[1])

# 3. 向後填充 -------------------------------
df_bfill = df.bfill()
error = mean_squared_error(df_orig['value'], df_bfill['value'])
df_ffill['value'].plot(title=f"Forward Fill (MSE: {error:.3f})", ax=ax[2])

# 4. 線性插值 -------------------------------
df['rownum'] = np.arange(df.shape[0])
df_nana = df.dropna(subset=['value'])

f = interp1d(df_nana['rownum'], df_nana['value'])
df['linear_fill'] = f(df['rownum'])

error = mean_squared_error(df_orig['value'], df['linear_fill'])
df['linear_fill'].plot(title=f"Linear Fill (MSE: {error:.3f})", ax=ax[3])

# 5. 三次方插值 -------------------------------
# df['rownum'] = np.arange(df.shape[0])
# df_nana = df.dropna(subset=['value'])
f = interp1d(df_nana['rownum'], df_nana['value'], kind='cubic')
df['cubic_fill'] = f(df['rownum'])

error = mean_squared_error(df_orig['value'], df['cubic_fill'])
df['cubic_fill'].plot(title=f"cubic Fill (MSE: {error:.3f})", ax=ax[4])


# 6. 鄰近均值插值 -------------------------------
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 out


df['knn_mean'] = knn_mean(df.value.values, 8)
error = mean_squared_error(df_orig['value'], df['knn_mean'])
df['knn_mean'].plot(title=f"KNN Mean (MSE: {error:.3f})", ax=ax[5])


# 7.季節性均值插值 -------------------------------

def seasonal_mean(ts, n, lr=0.7):
    """
    計算時間序列每一個季節性的均值,填充到季節性裏面的缺失值。
    ts:一個時間序列;
    n: 季節性的周期
    """

    out = np.copy(ts)
    for i, val in enumerate(ts):
        if np.isnan(val):
            ts_seas = ts[i - 1::-n]
            if np.isnan(np.nanmean(ts_seas)):
                ts_seas = np.concatenate(ts[i - 1::-n], ts[i::n])
            out[i] = np.nanmean(ts_seas) * lr

    return out


df['seasonal_mean'] = seasonal_mean(df.value, n=12, lr=1.25)
error = mean_squared_error(df_orig['value'], df['seasonal_mean'])
df['seasonal_mean'].plot(title=f"Seasonal Mean (MSE: {error:.3f})", ax=ax[6])

填充缺失值還有很多方法,填充的效果是取决於你對填充精度的要求。

  1. 如果還有很多輔助變量,可以使用隨機森林或者k鄰近等模型做預測填充;
  2. 可以使用過去的數據填充未來的;可以使用現在的數據去填充過去的。
  3. 利用時間序列的季節性填充數據;
  4. 利用時間序列的周期性對數據u哦填充。

11.什麼是自相關和偏自相關

  1. 自相關是時間序列和它自己滯後n階時序數據的相關性。如果一個時間序列有顯著的自相關性,說明使用曆史數據做預測效果會非常好。
  2. 偏自相關性是和自相關系數差不多,但是不包括中間滯後的相關貢獻。
from statsmodels.tsa.stattools import acf, pacf
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

df = pd.read_csv('../datasets/a10.csv')

# Calculate ACF and PACF upto 50 lags
# acf_50 = acf(df.value, nlags=50)
# pacf_50 = pacf(df.value, nlags=50)

# Draw Plot
fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(163), dpi=100)
plot_acf(df.value.tolist(), lags=50, ax=ax[0])
ax[0].set_title("acf")
plot_pacf(df.value.tolist(), lags=50, ax=ax[1])
ax[1].set_title("pacf")

11.1偏自相關系數到底是什麼

有一個時間序列叫 ,他的一階滯後叫 ,他的二階滯後叫 ,他的三階滯後叫

這裏有一個公式: .

那麼這裏 就是第三階偏自相關系數。

12.滯後項畫圖

滯後項畫圖,就是一個時間序列和他自己滯後n階的散點圖,通常可以用來檢查自相關性,如果可以從圖中看到一些相關性,說明數據有一點相關性的。如果從滯後項圖看不出來什麼數據信息,很有可能說明數據是白噪聲。

from pandas.plotting import lag_plot

# import data
ss = pd.read_csv("../datasets/sunspotarea.csv")
a10 = pd.read_csv("../datasets/a10.csv")

# plot Sun Spot Area
fig, ax = plt.subplots(ncols=4, nrows=1, figsize=(103), sharex=True, sharey=True, dpi=100)

for i, ax in enumerate(ax.flatten()):
    lag_plot(ss['value'], lag=i + 1, ax=ax, c='firebrick')
    ax.set_title(f"Lag ({i + 1})")
fig.suptitle("Lag Plot of Sun Spot Area\n(Points get wide and scattered with increasing lag -> lesser correlation)",
             fontdict=font, y=1.12)

# plot  Drugs Sales
fig, ax = plt.subplots(ncols=4, nrows=1, figsize=(103), sharex=True, sharey=True, dpi=100)

for i, ax in enumerate(ax.flatten()):
    lag_plot(a10['value'], lag=i + 1, ax=ax, c='firebrick')
    ax.set_title(f"Lag ({i + 1})")
fig.suptitle("Lag Plots of Drug Sales", fontdict=font, y=1.05)

從下圖可以觀察到,sun spot的滯後階數遠小於季節性,所以滯後1、2、3、4的數據和原始數據看起來沒什麼關系,可以說明時間序列前後可能沒有什麼相關性。 但是下面的毒品售賣滯後1、2、3、4階和原始數據看起來有很强的相關性,需要向下繼續研究。

13.如何對時間序列做平滑

對時間序列做平滑的優點:

  1. 减少噪聲對時間序列的影響;獲得除去噪聲後的時間序列數據;
  2. 時間序列平滑後的數據可以用來解釋原序列的一些特征;
  3. 更好的可視化潜在的趨勢。

如何平滑時間序列:

  1. 移動平均;
  2. 做一個局部平滑(LOESS);
  3. 做一個局部加權回歸(LOWESS )。

做移動平均的時候,需要注意窗口的大小。如果窗口過大,導致時間序列過於平滑,以至於還消除了一些其他信息。

LOESS是局部回歸的意思,取局部數據做回歸;LOWESS是局部加權回歸。

from statsmodels.nonparametric.smoothers_lowess import lowess

# 導入數據
df_orig = pd.read_csv("../datasets/elecequip.csv", parse_dates=['date'], index_col=['date'])

# 1. 移動平均
df_ma = df_orig.value.rolling(3, center=True, closed='both').mean()

# 2. loess 平滑 (5% and 15%)
# df_loess_5 = pd.DataFrame(lowess(df_orig.value, np.arange(len(df_orig.value))[:,1]), index=df_orig.index, columns=['value'])
df_loess_5 = pd.DataFrame(lowess(endog=df_orig.value, exog=np.arange(len(df_orig.value)), frac=0.05)[:,1],index=df_orig.index, columns=['value'])
df_loess_15 = pd.DataFrame(lowess(endog=df_orig.value, exog=np.arange(len(df_orig.value)), frac=0.15)[:,1],index=df_orig.index, columns=['value'])

# plot
fig, ax = plt.subplots(4,1, figsize=(77), sharex=True, dpi=100)
df_orig['value'].plot(ax=ax[0], color='k', title='Original Series')
df_loess_5['value'].plot(ax=ax[1], title='Loess Smoothed 5%')
df_loess_15['value'].plot(ax=ax[2], title='Loess Smoothed 15%')
df_ma.plot(ax=ax[3], title='Moving Average (3)')
fig.suptitle('How to Smoothen a Time Series', y=0.95, fontsize=14)

版权声明:本文为[yuanzhoulvpi]所创,转载请带上原文链接,感谢。 https://gsmany.com/2021/09/20210919231534635w.html