回顾

我们之前做了什么：

• 我们观察到时间序列的 基本成分 有 Trend, Seasonal, Residual，并由此发展了 滑动窗口 方法和 STL 分解；
• 我们观察到有的时间序列具有 自相关性，建立了 AR 模型，并为此发展了 ACF/PACF 试验，引入了 平稳性，进行了 ADF 检验 和 平稳化；
• 我们沿用了 线性回归；
• 我们组合基本成分、AR 模型、MA 滑动窗口的想法，最终得到 SARIMAX 模型

这次我们希望从另一个角度处理问题：数据的不靠谱性；真实世界的数据总是包含因为传感器质量等因素造成的误差甚至错误甚至 NaN，我们希望在时间序列预测中能滤去这部分噪音，为此，需要利用 Kalman 滤波。

控制领域的 Kalman 滤波

我们的目的是多快好省地把它应用在时间序列预测中，因而略去了公式推导，这不太好，但属无奈之举。

Kalman 滤波是一种能对抗测量和模型的不确定性（测量噪声、过程噪声）的预测算法，思想在于：

• 利用观测值修正预测值，得到一个更接近真实值的估计
• 通过 预测 -> 修正 -> 预测 -> 修正 ... 不断递推

建立模型

我们设定过程模型 Statement Equation：

其中

• \mathbf{x}_k 是 k 时刻的真实值 (未知的)
• {\displaystyle \mathbf {F} _{k}} 是状态转移矩阵，{\displaystyle \mathbf {B} _{k}} 是控制矩阵，{\displaystyle \mathbf {w} _{k}\sim {\mathcal {N}}\left(0,\mathbf {Q} _{k}\right)} 是服从多元正态分布的过程噪声

设定观测模型 Observation Equation：

其中

• \mathbf{z}_k 是 k 时刻的观测值
• {\displaystyle \mathbf {H} _{k}} 是观测转移矩阵
• {\displaystyle \mathbf {v} _{k}\sim {\mathcal {N}}\left(0,\mathbf {R} _{k}\right)} 是服从正态分布的测量噪声

上面的转移矩阵会根据我们先验地对系统建模而确立。

Step 1. Prediction 预测

根据上一时刻的系统状态，用过程模型得到这一时刻系统状态的预测：

Step 2. Correction 修正

我们想知道观测值应该以多大的权重影响预测值？一个合理的想法是方差越大，越不靠谱，影响越小；具体而言，我们想使得这一步更新后的预测误差协方差的 MSE 最小，从而利用先验误差协方差 \mathbf{P}_{k|k-1} 和观测噪声协方差 \mathbf{R}，计算 Kalman Gain

并使用 Kalman Gain 修正预测值（会被传递到下一时刻，作为估计的这时状态真实值）

并更新先验误差协方差矩阵

不断进行下去即可。

时间序列预测的 Kalman 滤波

考虑 Kaggle 的一个例子

Daily Power Production of Solar Panels

In oktober 2011 we installed solar pannels (or Photovoltaic Modules) on our roof. The total power of the modules is 5kWp.

It may seem strange but we have the habbit of making daily recordings of our electricity usage and so it was obvious to make record of the powerproduction of the solar pannels. I am trying to predict the date of the next 1000kWh produced.

这个数据集有很明显的噪声和异常值

我们试着用 Kalman 滤波处理它，这里的问题是要向滤波器提交状态空间模型，所以首先，我们试着把一个 ARMA(p, q) \displaystyle y_t = \sum_{i=1}^p \phi_i y_{t-i} + \sum_{j=1}^q \theta_j \varepsilon_{t-j} + \varepsilon_t 改写为状态空间模型 \displaystyle \left\{ \begin{aligned} & \alpha_t = T\alpha_{t-1} + R\mu_t \\ & y_t = Z\alpha_t + \varepsilon_t \end{aligned} \right.（引入隐状态 \alpha_t，过程噪声和观测噪声）

使用 Harvey 方法可以完成这件事，这里给出结果（字母用法和前面不同）：

并取 Q=\begin{bmatrix} \sigma^2 & 0 & \cdots \\ 0 & 0 & \cdots \\ \vdots & \vdots & \ddots\end{bmatrix}, R=[0]

将上面的状态空间模型输入 Kalman 滤波器即可。

附录

Kalman + ARMA(p, q)

import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
from statsmodels.tsa.arima.model import ARIMA
from pmdarima import auto_arima

class KalmanFilter:
    def __init__(self, F, B, H, Q, R, x0, P0):
        self.F = F
        self.B = B  
        self.H = H  
        self.Q = Q  
        self.R = R  
        self.x = x0  
        self.P = P0  
    def predict(self, u):
        self.x = np.dot(self.F, self.x) + np.dot(self.B, u)
        self.P = np.dot(np.dot(self.F, self.P), self.F.T) + self.Q
        return self.x
    def update(self, z):
        S = np.dot(np.dot(self.H, self.P), self.H.T) + self.R
        K = np.dot(np.dot(self.P, self.H.T), np.linalg.inv(S))
        y = z - np.dot(self.H, self.x)  
        self.x = self.x + np.dot(K, y)
        I = np.eye(self.P.shape[0])
        self.P = np.dot(np.dot(I - np.dot(K, self.H), self.P), (I - np.dot(K, self.H)).T) + np.dot(np.dot(K, self.R), K.T)
        return self.x

def param(data):
    model = auto_arima(data, d=0, seasonal=False, trace=True)
    print(f"ARIMA(p, d=0, q) Choice: {model.order}")
    return model.order

def predict(data, kf):
    predictions = []
    for i in range(len(data)):
        x = kf.predict(np.array([[0]]))
        x = kf.update(data.iloc[i, 0])
        predictions.append(kf.x[0, 0])
        
    return predictions

if __name__ == "__main__":
    original = pd.read_csv("PV_Elec_Gas3.csv")
    original['date'] = pd.to_datetime(original['date'])
    original = original.set_index('date')
    data = original.iloc[:, 1:2]
    data_diff = data.diff().dropna()

    p, d, q = param(data_diff)
    model = ARIMA(data_diff, order=(p, d, q)).fit()
    print(model.summary())

    ar_params = model.arparams
    ma_params = model.maparams
    sigma2 = model.mse

    r = max(p, q+1)
    F = np.zeros((r, r))
    F[0, :p] = ar_params
    F[1:, :-1] = np.eye(r-1)
    B = np.zeros((r, 1))
    H = np.zeros((1, r))
    H[0, 0] = 1
    H[0, 1:] = ma_params
    Q = np.zeros((r, r))
    Q[0, 0] = sigma2
    R = np.zeros((1, 1))
    X0 = np.zeros((r, 1))
    P0 = np.eye(r) * sigma2

    kf = KalmanFilter(F, B, H, Q, R, X0, P0)
    predictions = predict(data_diff, kf)
    data['predictions'] = pd.Series(predictions, index=data_diff.index).cumsum() + data.iloc[0, 0]

    plt.figure()
    plt.plot(data.index, data.iloc[:, 0], label='Original', color='blue')
    plt.plot(data.index, data.iloc[:, 1], label='Filtered', color='red')

    plt.legend()
    plt.show()

记录本笔记的过程中参考了这些资料：

1. https://www.kalmanfilter.net/
2. https://zhuanlan.zhihu.com/p/649997859
3. https://en.wikipedia.org/wiki/Kalman_filter
4. https://www.geeksforgeeks.org/kalman-filter-in-python/