ts1101

1 Exponential smoothing & Kalman filter

  • Time Update (“Predict”)

\[\hat{x}_k^- ~=~ \hat{x}_{k-1}\]

\[P_k^- ~=~ P_{k-1}+Q\]

  • Measurement Update (“Correct”)

\[K_{k_{~~2×1}}=P_k^- \mathbf{}\left[\begin{array}{rrr} 1 ~\\ Z_t \end{array}\right] \mathbf{}\left[\begin{array}{rrr} \mathbf{}\left[\begin{array}{rrr} 1 & Z_t\end{array}\right] ~P_k^- ~ \mathbf{}\left[\begin{array}{rrr} 1 ~\\ Z_t \end{array}\right]_{2×1} +R \end{array}\right]^{-1}\]

\[\hat{x}_k ~=~ \hat{x}_k^-+K_k(x_k-\mathbf{}\left[\begin{array}{rrr} 1 & Z_t \end{array}\right]\hat{x}_k^-)\]

\[P_k ~=~ (I_{2×2} - K_k \mathbf{}\left[\begin{array}{rrr} 1 & Z_t \end{array}\right] )P_k^-\]

  • Exponential smoonthing

\[y_t = \mathbf{}\left[\begin{array}{rrr} 1 & Z_t \end{array}\right] \hat{x}_k +R\]

library(tidyverse)
library(lubridate)
library(foreach)
library(forecast)

y = read.csv("/home/yuan/Desktop/meeting/處理後來店人數資料.csv")$x

date1 = foreach(i = 0:363, .combine = "c") %do% {
    (ymd_h("2007-01-01-12") + hours(0:6) + days(i))
}

data1 = data.frame(data = y, date = date1, week = weekdays(date1), zt = 0, stringsAsFactors = F)


data1$zt[which(data1$week == "Saturday")] = 1
data1$zt[which(data1$week == "Sunday")] = 1
KF <- function(alpha = 0.5, R = 0.5, lt = 1, pt = 1, P = matrix(c(1, 0, 0, 1), 
    nrow = 2, ncol = 2)) {
    # xhat & P 要起始值， alpha & R 自設
    
    # lt & pt state value
    
    # zt regressor (0,1) 1*2548
    
    # xk true data 1*2548
    alpha = alpha
    R = R
    A = matrix(c(1, 0, 0, 1), nrow = 2, ncol = 2)
    I = matrix(c(1, 0, 0, 1), nrow = 2, ncol = 2)
    Q = matrix(c(alpha^2 * R, 0, 0, 0), nrow = 2, ncol = 2)
    lt = lt
    pt = pt
    
    # zt regressor 1*2548
    zt = matrix(data1$zt, nrow = 1, ncol = 2548)
    # xk true value 1*2548
    xk = matrix(data1$data, nrow = 1, ncol = 2548)
    
    xhat_update = list()
    P_update = list()
    
    P = P
    xhat = matrix(c(lt, pt), nrow = 2, ncol = 1)
    
    for (i in 1:2548) {
        # time update
        xhat_minus = A %*% xhat
        Pminus = A %*% P %*% t(A) + Q
        
        H = matrix(c(1, zt[1, i]), nrow = 1, ncol = 2)  # H  1*2 matrix
        # measurement update
        K = Pminus %*% t(H) %*% (H %*% Pminus %*% t(H) + R)^(-1)
        
        xhat = xhat_minus + K %*% (xk[i] - H %*% xhat_minus)
        
        P = (matrix(c(1, 0, 0, 1), nrow = 2, ncol = 2) - K %*% H) %*% Pminus
        
        xhat_update[[i]] = xhat
        # P_update[[i]] = P
    }
    return(xhat_update)
}

1.1 R = 0.01, pt = 1, lt = 1, alpha = 0.5

a = KF(R = 0.01, pt = 1, lt = 1, alpha = 0.5)
test1 = unlist(sapply(a, "[[", 1)) + data1$zt * unlist(sapply(a, "[[", 2)) + 
    0.01
accuracy(test1, data1$data)
                ME     RMSE      MAE       MPE     MAPE
Test set 0.2822206 27.18353 18.88427 -9.999182 23.37479

1.2 R = 0.01, pt = 1, lt = 1, alpha = 0.9

a = KF(R = 0.01, pt = 1, lt = 1, alpha = 0.9)
test1 = unlist(sapply(a, "[[", 1)) + data1$zt * unlist(sapply(a, "[[", 2)) + 
    0.01
accuracy(test1, data1$data)
                ME     RMSE      MAE       MPE     MAPE
Test set 0.1443447 18.56933 12.66392 -6.616363 16.03923
matplot(data1$data, type = "l")
matplot(test1, type = "l", add = T, col = 2)

                       RMSE      MAE     MAPE
Holt.Winter        32.18209 21.82968 19.20088
double.seasonal.HW 23.43904 16.81247 16.16967
kalman_filter      18.56933 12.66392 16.03923

2 Kalman Filter

  • Time Update (“Predict”)

\[\hat{x}_k^- ~=~ \hat{x}_{k-1}\]

\[P_k^- ~=~ P_{k-1}+Q\]

  • Measurement Update (“Correct”)

\[K_k ~=~ \frac{P_k^-}{P_k^-+R}\]

\[\hat{x}_k ~=~ \hat{x}_k^-+K_k(Z_k-\hat{x}_k^-)\]

\[P_k ~=~ (1-K_k)P_k^-\]

kalman_filter = function(n_iter = 100, R = 0.01, means = -0.37727, sd = 0.1) {
    
    n_iter = n_iter
    x = means  # truth value
    z = rnorm(x, sd, n = n_iter)  # observations
    
    Q = 1e-05  # process variance
    
    xhat = c(rep(0, n_iter))  # a posteri estimate of x
    P = c(rep(0, n_iter))  # a posteri error estimate
    xhatminus = c(rep(0, n_iter))  # a priori estimate of x
    Pminus = c(rep(0, n_iter))  # a priori error estimate
    K = c(rep(0, n_iter))  # Kalman gain
    
    R = R  # estimate of measurement variance
    
    # intial
    xhat[1] = 0
    P[1] = 1
    
    for (k in 2:n_iter) {
        # time update
        xhatminus[k] = xhat[k - 1]
        Pminus[k] = P[k - 1] + Q
        
        # measurement update
        K[k] = Pminus[k]/(Pminus[k] + R)
        xhat[k] = xhatminus[k] + K[k] * (z[k] - xhatminus[k])
        P[k] = (1 - K[k]) * Pminus[k]
    }
    matplot(z, ylim = c(means - 3 * sd, means + 3 * sd), pch = 3)
    abline(a = means, b = 0, col = 2)
    matplot(xhat, type = "l", add = T)
}

kalman_filter(n_iter = 50, R = 1e-04, means = -0.37727, sd = 0.1)

2017-11-01