0
$\begingroup$

Following some examples on Chad Fulton's blog and in statsmodels' tests, I have tried to come up with an equivalent of a pykalman implementation.

The original question was deemed unclear and was requested to be edited.

In the meanwhile, I familiarised myself a bit more with Kalman fiters and partially answered the question.

The below is thus a mix bag of a reformulated question and a less-than-ideal answer.

First, a back-ground of the specific example.

What we are attempting here is sometimes called "Sensor fusion". We attempt to observe the same hidden state from several noisy sources and try to rebuild it.

Mathematically and with some common notation in the Kalman litterature, the problem is stated as follow:

  1. State propagation

    we make no prediction here. We just observe a scalar quantity and so, if my understanding is correct, the below is intended to say that the best prediction we can make is the previous observation: \alpha_{t} = 1 \alpha_{t-1} + \epsilon_{t} where:

    • \alphat_{t} is the hidden state at time t
    • \epsilon_{t} ~ N(0, H)
    • H: Covariance matrix of state
  2. Observation

    We have several sensors, all noisy, observing the same single state: y_{t} = Z_{t} \alpha_{t} + R_t \eta_{t}

    • z_t{} = column vector of 1s with dimension the number of observations.
    • \eta_{t} ~ N(0,Q)
    • where Q is the observation covariance matrix

Now for the code.
In order to give a workable example as requested by the diligent editors of Cross Validated, I need a few preliminaries.

Create a hidden state to be observed

For no special reason, we pick a Geometric Brownian Motion.

import pykalman
import statsmodels as sm
import statsmodels.tsa
import statsmodels.tsa.statespace
import statsmodels.tsa.statespace.mlemodel
import numpy as np
import pandas as pd
import pendulum as pdl


def brownian_motion(T = 1, N = 100, mu = 0.1, sigma = 0.1, S0 = 20):
    dt = float(T)/N
    t = np.linspace(0, T, N)
    W = np.random.standard_normal(size = N)
    W = np.cumsum(W)*np.sqrt(dt) # standard brownian motion
    X = (mu-0.5*sigma**2)*t + sigma*W
    S = S0*np.exp(X) # geometric brownian motion
    return S
# strftime is to get rid of the nanosecond issue between pendulum 
# and pandas (taken from Stackoverflow som
dates = pd.date_range(pdl.today().strftime("%Y-%m-%d %H:%M:%S %z"), periods=365, freq='d')
T = (dates.max()-dates.min()).days / 365
N = dates.size
start_value = 100
y = pd.Series(brownian_motion(T, N, sigma=0.1, S0=start_value), index=dates)

Simulate observations with noise

nb_obs = 8
stdevs = np.arange(nb_obs)*2+2

dobs = {}
for i in range(nb_obs):
    dobs[i]=y+np.random.normal(0,stdevs[i], y.shape[0])
dbm_obs={'BM': y, 'obs': obs}
brownian_and_obs = pd.concat([y,obs],axis=1)
brownian_and_obs.columns = ['BM']+list(map(lambda s: 'obs'+str(s), range(nb_obs)))

Two implementations of the Kalman filter

nb_srcs = 2 #nb_obs
states_initial = brownian_and_obs.filter(regex='obs.*',axis=1).iloc[0,:nb_srcs].values
states_all =brownian_and_obs.filter(regex='obs.*',axis=1).iloc[:,:nb_srcs].values

#
# pykalman
#
kf_pyk = pykalman.KalmanFilter(transition_matrices=[1.],

                                   observation_matrices=np.ones(nb_srcs).reshape(nb_srcs,1),

                                   initial_state_mean=np.mean(states_initial),

                                   initial_state_covariance=[1.],

                                   observation_covariance=np.eye(nb_srcs))



states_kfpk, _ = kf_pyk.filter(states_all)
states_kfpk = states_kfpk.reshape(states_kfpk.shape[0])
states_kfpk = pd.Series(states_kfpk,index=y.index)


#
# State model
#

mod_current = sm.tsa.statespace.mlemodel.MLEModel(states_all, k_states=1,initialization='stationary', k_posdef=1)

# y_{t} = Zt αt+ dt + εt

# Z : design (k_endog×k_states×nobs)

# d : obs_intercept (k_endog×nobs)

# H : obs_cov (k_endog×k_endog×nobs)

mod_current['design'] = np.ones(nb_srcs).reshape(nb_srcs,1,1)

mod_current['obs_cov'] = 100*np.eye(nb_srcs)



# αt = Tt α{t−1}+ct+Rt ηt

# T : transition (k_states×k_states×nobs)

# c : state_intercept (k_states×nobs)

# R : selection (k_states×k_posdef×nobs)

# Q : state_cov (k_posdef×k_posdef×nobs)

#mod_current['transition'] = np.ones(1) GOT AN ERROR WITH THIS UNCOMMENTED

# mod_current['state_intercept'] = np.zeros(a.shape[0]).reshape(1,a.shape[0]) WORKS BUT USELESS

mod_current['selection'] =np.ones(1)#.reshape(1,a.shape[0])#np.ones(1)

mod_current['state_cov'] = [10000]  # Variance set high initially if unknown

#mod_current.initialize_known(1, [[0.02]])

#mod_current['state_intercept'] = [[0]]#np.r_[controls[1:], np.nan][None,:]

res=mod_current.filter(states_initial)

states_kfsm = res.filtered_state
states_kfsm = pd.Series(states_kfsm[0],index=y.index)

Visual results

plt.figure(figsize=(15,10))
estimates = pd.concat([y,states_kfsm, states_kfpk],axis=1)
estimates.columns=['BM',  'KalmanSM', 'KalmanPyK']
sns.lineplot(data = estimates.loc(axis=1)['BM','KalmanSM'],dashes=False)
sns.scatterplot(data=brownian_and_obs.filter(regex='obs.*',axis=1).iloc[:,:nb_srcs]);

Two observations (markers), Kalman estimate (orange line) and true (hidden) state (blue line) Comparison of the two Kalman filters

As said, I believe I partly answered my original question: the pykalman and the statsmodels implementation seem to yield similar results.

Now, it has raised others, which I list in order to turn this post into a question:

  • I cannot un-comment the state transition in the statsmodels version. It yields an error I am not sure why.
  • the Kalman seems sometimes far away (see February 2020 for example). Besides, increasing the number of sources increases the standard deviation of the estimate and I was under the impression that the algorithm has the potential not to degrade when more noisy sources are added. Is that a correct implementation of the problem described above ?
$\endgroup$
  • 1
    $\begingroup$ your code isn't reproducible, but I find it odd that in the first example you're setting the 'obs_cov' to np.eye(nb_srcs_current) while in the second example you're setting the observation matrices to the same thing. We also might be able to say more if you $\LaTeX$ up your model $\endgroup$ – Taylor Apr 6 at 15:56
  • $\begingroup$ Hi Taylor, the intention was to have a transition of the form: d y_t = d \alpha_t an. d d \alpha_t = d W_t where d W_t is a nb_src_current brownian motion. So a very basic model actually. The pykalman seems to be a ok job at filtering the sources when I look at the result. $\endgroup$ – Trad Dog Apr 6 at 17:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.