技術標籤：資料類序列資料python資料分析

移動視窗函式

移動視窗函式可以理解為時FIR濾波器，只不過這裡是濾波器在運動，而不是訊號在運動。但是從相對運動的角度來說，移動視窗函式就是FIR濾波器

資料樣例

載入一個時間序列，並且按照工作日頻率進行重新取樣。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
# %matplotlib inline

close_px_all = pd.read_csv("data/stock_px_2.csv") 
close_px_all.
 
head()

# 簡單檢查日期是否存在重複值, 確保資料彙總到day這個級別
_date_col = close_px_all.iloc[0,:]
if len(_date_col) == _date_col.nunique():
    print("沒有日期重複"
 
)
else:
    print("存在日期重複")

沒有日期重複

# 缺失值檢視: 全為零，沒有缺失值
close_px_all.isnull().sum()

Unnamed: 0    0
AAPL          0
MSFT          0
XOM           0
SPX           0
dtype: int64

close_px_all = pd.read_csv("data/stock_px_2.csv",
                          parse_dates=True,      # -> 自動解析日期，不需要專門指定日期的字串格式
 

                          index_col=0)           # -> 0列作為索引列
close_px = close_px_all[['AAPL', 'MSFT', 'XOM']] # -> AAPL：蘋果，MSFT：微軟，XOM：埃克森美孚

close_px = close_px.resample('B').ffill()        # -> 按照工作日頻率進行重新取樣

# 檢視蘋果股價
close_px.AAPL.plot()
plt.show()

圖1 蘋果股價

250日移動視窗平均值

均值濾波，濾波器的長度為250天

appl_mean250 = close_px.AAPL.rolling(250).mean()
close_px.AAPL.plot()
appl_mean250.plot()
plt.legend(['raw', 'mean-250'])
plt.show()

圖2 平均濾波和原始資料

仔細觀察上圖可以發現，橙色的部分有一段是沒有值的。這個是由於序列前249個值長度小於移動視窗（250），無法進行這個視窗下的均值計算，所以這種情形下的移動視窗計算結果的序列長度要小於原序列長度。

appl_mean250[:249]  # -> 前 249 個值為nan

2003-01-02   NaN
2003-01-03   NaN
2003-01-06   NaN
2003-01-07   NaN
2003-01-08   NaN
              ..
2003-12-10   NaN
2003-12-11   NaN
2003-12-12   NaN
2003-12-15   NaN
2003-12-16   NaN
Freq: B, Name: AAPL, Length: 249, dtype: float64

參考FIR濾波器的內容，可以通過序列頭部補零來確保移動視窗計算前後的序列長度保持一致。

# 補零
_tmp = pd.DataFrame({'AAPL':[0]*249,
                    'MSFT':[0]*249,
                    'XOM':[0]*249}, 
                   index=pd.date_range(end='2003-01-01', periods=249, freq='B'))

close_px_zero_fill = pd.concat([_tmp,close_px])

# 補零後求250移動視窗均值
appl_mean250_zero_fill = close_px_zero_fill.AAPL.rolling(250).mean()
close_px.AAPL.plot()
appl_mean250_zero_fill[249:].plot()     # -> 擷取原始資料部分
plt.legend(['raw', 'mean-250'])
plt.show()

圖3 補零後的均值濾波

觀察到上圖中，mean-250和raw這兩條線的長度一致了。

# 按照上面的補零做法，我們很容易得出結論：
# appl_mean250_zero_fill第1個值等於close_px.AAPL[0]/250
if appl_mean250_zero_fill[249] == close_px.AAPL[0]/250:
    print("appl_mean250_zero_fill第1個值等於close_px.AAPL[0]/250")
else:
    print("不相等")

print("appl_mean250_zero_fill[249]：{}, close_px.AAPL[0]:{}, close_px.AAPL[0]/250: {}".format( \
        appl_mean250_zero_fill[249], close_px.AAPL[0], close_px.AAPL[0]/250))

appl_mean250_zero_fill第1個值等於close_px.AAPL[0]/250
appl_mean250_zero_fill[249]：0.0296, close_px.AAPL[0]:7.4, close_px.AAPL[0]/250: 0.0296

# 對上述補零內容寫成函式
def fill_zero(df:pd.DataFrame, win:int, end_date='2003-01-01',freq='B'):
    """
    對輸入時間序列df頭部補win-1個零
    
    引數
    -----
    ：df：時間序列，index為時間
    ：win：視窗函式大小
    ：end_date：比時間序列df的起始時間早一個時間單位
    ：freq：時間單位，'B'表示時間單位是工作日的一天
    
    返回
    -----
    補零後的時間序列
    """
    if not isinstance(df,pd.DataFrame):
        df = pd.DataFrame(df)
        
    # 補零
    _tmp = pd.DataFrame(np.zeros((win-1, df.shape[1])), 
                        index=pd.date_range(end=end_date, periods=win-1, freq='B'))
    _tmp.columns = df.columns    
    df_fill_zero = pd.concat([_tmp,df])

    return df_fill_zero

# 測試這個程式碼
df_fill_zero = fill_zero(close_px.AAPL, 250)
plt.figure()
plt.plot(close_px.AAPL)
plt.plot(df_fill_zero.rolling(250).mean()[249:])
plt.legend(['raw', 'mean-250'])
plt.show()

圖4 補零後的均值濾波

df_fill_zero = fill_zero(close_px, 250)
plt.figure()
plt.plot(close_px)
plt.plot(df_fill_zero.rolling(250).mean()[249:])
plt.legend(['raw-AAPL', 'raw-MSFT', 'raw-XOM', 'mean-250-AAPL','mean-250-MSFT', 'mean-250-XOM'])
plt.show()

圖5 多個股票補零後均值濾波結果

移動視窗標準差、方差、偏度、峭度

df_fill_zero = fill_zero(close_px, 250)
_std250 = df_fill_zero.rolling(250).std()[249:]
plt.figure()
plt.plot(_std250)
plt.legend(['std-250-AAPL','std-250-MSFT', 'std-250-XOM'])
plt.show()

圖6 標準差

# 方差
df_fill_zero = fill_zero(close_px, 250)
_var250 = df_fill_zero.rolling(250).var()[249:]
plt.figure()
plt.plot(_var250)
plt.legend(['var-250-AAPL','var-250-MSFT', 'var-250-XOM'])
plt.show()

圖7 方差

# 偏度：補零對偏度計算來說是沒有意義的
# 只從不會有缺失值的地方開始
# skew > 0， 正偏態，尾巴在右面
_skew365 = close_px.rolling(365).skew()[364:]
plt.figure()
plt.plot(_skew365)
plt.legend(['skew-365-AAPL','skew-365-MSFT', 'skew-365-XOM'])
plt.show()

圖8 偏度

# 峭度：補零對峭度計算來說是沒有意義的
# 只從不會有缺失值的地方開始
# 2005年出現峭度值大於3，說明此階段股價的pdf屬於厚尾分佈（尾巴比相應的正態分佈尾巴厚）
# 峭度很大的原因是序列中存在大於均值很多的值（超過正態分佈的情況），
# 即存在衝擊型值
_kurt365 = close_px.rolling(365).kurt()[364:]
plt.figure()
plt.plot(_kurt365)
plt.legend(['kurt-365-AAPL','kurt-365-MSFT', 'kurt-365-XOM'])
plt.show()

圖9 峭度

最大值、最小和中位數

# 最大值
_max365 = close_px.rolling(365).max()[364:]
plt.figure()
plt.plot(_max365)
plt.legend(['max-365-AAPL','max-365-MSFT', 'max-365-XOM'])
plt.show()

圖10 最大值

# 最小值
_min365 = close_px.rolling(365).min()[364:]
plt.figure()
plt.plot(_min365)
plt.legend(['min-365-AAPL','min-365-MSFT', 'min-365-XOM'])
plt.show()

圖11 最小值

# 中位數
_median365 = close_px.rolling(365).median()[364:]
plt.figure()
plt.plot(_median365)
plt.legend(['median-365-AAPL','median-365-MSFT', 'median-365-XOM'])
plt.show()

圖12 中位數

眾數

通過roling().apply()方式自定義實現。

# 眾數
_mode365 = close_px.rolling(365).apply(lambda r: r.mode()[0])[364:]
plt.figure()
plt.plot(_mode365)
plt.legend(['mode-365-AAPL','mode-365-MSFT', 'mode-365-XOM'])
plt.show()

圖13 眾數

指數加權函式：又叫指數平滑

對時間序列來說，過去太久的觀測和近期的觀測對未來的影響程度肯定是不同的。
而指數加權通過向近期的觀測值提供更高的權重，來區分近期和更遠的歷史觀測之間的區別。

aapl_px = close_px.AAPL['2006':'2007']
ma60 = aapl_px.rolling(30, min_periods=20).mean()
ewma60 = aapl_px.ewm(span=30).mean()
ma60.plot(style='r--', label = 'simple MA')
ewma60.plot(style='k-', label='EW MA')
plt.legend()
plt.show()

圖14 指數平滑和簡單的平均值濾波的比較

二元移動視窗函式：相關度分析

檢視兩個時間序列的相關度，例如，股票對基準指數的關聯性。

spx_px = close_px_all['SPX']
# 基準指數的百分比變化
spx_rets = spx_px.pct_change()
# 三隻股票的百分比變化
returns = close_px.pct_change()

# 滾動計算三隻公司股票和基準指數相關性
corr = returns.rolling(125, min_periods=100).corr(spx_rets)
corr.plot()
plt.show()

圖15 三隻股票和基準指數相關度

參考

• python for data analysis:Data Wrangling with Pandas, Numpy, IPython, Wes Mckinney