【精选】卡尔曼滤波及其在配对交易中的应用

# 以焦炭的收盘价数据作为x，螺纹的收盘价数据作为y
# 螺纹价格 = alpha + beta * 焦炭价格 + 随机误差 
from pykalman import KalmanFilter

#建立观测矩阵
observation_matrices = np.vstack(( np.ones(len(df[:'2013'])),
                                  df.loc[:'2013','焦炭'].values )).T
Shape = observation_matrices.shape
observation_matrices = observation_matrices.reshape(Shape[0],1,Shape[1])

#定义卡尔曼滤波的方程
kf = KalmanFilter(transition_matrices=np.array([[1,0],[0,1]]), #转移矩阵为单位阵
                  observation_matrices=observation_matrices)
np.random.seed(0)

# 使用2013年以前的数据，采用EM算法，估计出初始状态，
# 初始状态的协方差，观测方程和状态方程误差的协方差
kf.em(df.loc[:'2013','螺纹'])


#对2013年的数据做滤波
filter_mean,filter_cov = kf.filter(df.loc[:'2013','螺纹'])#观测值为螺纹

#从2014年开始滚动
start_index = np.where(df.index.year==2014)[0][0]

for i in range(start_index,len(df)):
    observation_matrix = np.array([[1,df['焦炭'].values[i]]])
    observation = df['螺纹'].values[i]

    #以上一个时刻的状态，状态的协方差以及当前的观测值，得到当前状态的估计
    next_filter_mean,next_filter_cov = kf.filter_update(
            filtered_state_mean = filter_mean[-1],
            filtered_state_covariance = filter_cov[-1],
            observation = observation,
            observation_matrix = observation_matrix)

    filter_mean = np.vstack((filter_mean,next_filter_mean))
    filter_cov = np.vstack((filter_cov,next_filter_cov.reshape(1,2,2)))

#得到alpha和beta
alpha = pd.Series(filter_mean[start_index:,0], index = df.index[start_index:])
beta = pd.Series(filter_mean[start_index:,1], index = df.index[start_index:])