0
$\begingroup$

I am trying to learn and implement the Kalman filter. Yesterday I successfully implemented a non linear kalman filter of the form:

$$ x_t = a(x_{t-1}) + u_t \\ y_t = Gy_{t-1} + v_t $$

$u_t$ and $v_t$ are distributed normally with mean 0, and covariance matrices $Q$ and $R$.

$a$ is a logistic function in my specific case.

I simulated the true model first and then tried estimating the Kalman filter based on the observables and got a pretty neat looking fit.

enter image description here

Now I tried to go beyond that and see if I could replicate the exercise except this time I pretended I didn't know Q, R, mu_0 (mean of the distribution from which the initial draw is taken), or the parameters of a. My method was to:

  1. Make a function which implements the kalman filtering and prediction algorithm for a given set of parameters
  2. Make that function output a log likelihood function
  3. Make an optimizer minimize that log likelihood and hence get my parameters

I implemented this in python. Unfortunately, while my output looks similar to the underlying simulation, its far from perfect. I have attached a picture of the the graph below.enter image description here

Here's the relevant python code for simulating the underlying process

import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt
from scipy import optimize

# define logistic function

def logistic_f(x_in, k = 100, dt = 0.1):

    r = float(x_in[0])
    p = float(x_in[1])

    return k * p * np.exp(r*dt) / (k + p * (np.exp(r*dt) - 1))

# = simulating a model of logistic growth = #

# setting params
p0 = 10.0
r = 0.2

# generating data
obs_var = 25
t = 250

# population at each point
nu = np.random.normal(loc = 0.0, scale = np.sqrt(obs_var), size = t)
p_init = p0
pop = np.empty(t)
x = np.empty(t)

# getting underlying true growth
for i in range(t):

    x[i] = logistic_f([r, p_init])
    p_init = x[i]

# getting the noisy measure
pop = x + nu

# plotting
fig, ax = plt.subplots(figsize = (9, 6))
ax.plot(range(t), pop, label = 'Noisy Measure')
ax.plot(range(t), x, label = 'True Growth')
ax.set_title(label = "True Population at each point")
ax.legend()
plt.show()

And here is the relevant bit of code for implementing the kalman filter and the mle optimization.

# define A function

def A_f(x_in, k = 100, dt = 0.1):

    r = float(x_in[0])
    p = float(x_in[1])

    return np.array([[r],
                    [k * p * np.exp(r*dt) / (k + p * (np.exp(r*dt) - 1))]])


# == kalman updater == #

def kalman_update(x_init, sigma_init, y,
                 logistic_f, Af, G, Q, R,
                 dt = 0.1, k = 100):

    # Af is the function who's Jacobian yields A

    x_0, sigma_0= x_init, sigma_init

    # filter

    x_f = x_0 + sigma_0 @ G.T @ np.linalg.inv(G @ sigma_0 @ G.T + R) @ (y - G @ x_0)
    sigma_f = sigma_0 - sigma_0 @ G.T @ np.linalg.inv(G @ sigma_0 @ G.T + R) @ G @ sigma_0

    # setting up A
    # A2 is the second row of matrix A. The first row will be the jacobian of the constant function in this case. 
    # So first row = [df(r)/dr, df(r)/dx] = [1, 0]. We bind this with A2

    eps = np.sqrt(np.finfo(float).eps) # defining a small epsilon to calculate gradient

    A2 = optimize.approx_fprime(x_f.flatten(), logistic_f, [eps, np.sqrt(200) * eps], k, dt)
    A = np.array([[1, 0],
                    A2])

    K_0 = A @ sigma_f @ G.T @ np.linalg.inv(G @ sigma_f @ G.T + R)
    x_1 = Af(x_f) + K_0 @ (y - G @ x_f)
    sigma_1 = A @ sigma_f @ A.T - K_0 @ G @ sigma_f @ A.T + Q

    # error and error covariance
    v = y - G @ x_1
    F = G @ sigma_1 @ G.T + R

    return x_1, sigma_1, v, F


def likelihood(params):

    # unpacking unknown params
    k = params[0]
    dt = params[1]
    Q = np.array([[params[2], params[3]],
                  [params[4], params[5]]])
    R = params[6]
    x_0 = np.array([[params[7], params[8]]]).T


    # known params
    G = np.array([[0, 1]])
    sigma_0 = np.array([[144, 0],
                    [0, 25]])
    y = pop

    # initializing empty matrices for storage
    kf_x = np.empty((2, t))
    kf_sigma = np.empty((2, t))
    loglike = np.empty(t)

    for i in range(t):

        kf_x[:, i] = x_0.flatten()

        y_i = y[i]
        x_1, sigma_1, v, F = kalman_update(x_init = x_0, 
                                     sigma_init = sigma_0, 
                                     y = y_i,
                                     logistic_f = logistic_f,
                                     Af = A_f,
                                     G = G,
                                     Q = Q,
                                     R = R)

        loglike[i] = 0.5 * (t * np.log(2 * np.pi) + np.log(np.linalg.det(F)) + v.T @ np.linalg.inv(F) @ v);
        x_0, sigma_0 = x_1, sigma_1

    L = np.sum(loglike)
    return L


sigma_0 = np.array([[144, 0],
                    [0, 25]])
y = pop


init = [1] * 9
params = optimize.minimize(likelihood, init, method='SLSQP')

And finally, here's the relevant bit of code for running the kalman filter with the parameters generated by the optimizer.

# == setting up kalman process == #
k = params['x'][0]
dt = params['x'][1]
Q = np.array([[params['x'][2], params['x'][3]],
              [params['x'][4], params['x'][5]]])
R = params['x'][6]
x_0 = np.array([[params['x'][7], params['x'][8]]]).T


sigma_0 = np.array([[144, 0],
                    [0, 25]])

kf_x = np.empty((2, t))
kf_sigma = np.empty((2, t))

loglike = np.empty(t)

y = pop # for ease

for i in range(t):

    kf_x[:, i] = x_0.flatten()

    y_i = y[i]
    x_1, sigma_1, v, F = kalman_update(x_init = x_0, 
                                 sigma_init = sigma_0, 
                                 y = y_i,
                                 logistic_f = logistic_f,
                                 Af = A_f,
                                 G = G,
                                 Q = Q,
                                 R = R)

    loglike[i] = -0.5 * (t * np.log(2 * np.pi) + np.log(np.linalg.det(F)) + v.T*np.linalg.inv(F)*v);
    x_0, sigma_0 = x_1, sigma_1


fig, ax = plt.subplots(figsize = (9, 6))
ax.plot(range(t), kf_x[1, :], label = 'Kalman Filter + MLE', marker='x', markersize=3, linestyle = 'None')
ax.plot(range(t), x, label = 'True Process', alpha = 0.7)
ax.plot(range(t), y, label = 'Observable')
ax.set_title('Population Growth')
ax.legend()
plt.show()

Sorry for the overly long post, but I can't seem to figure out what to do here. For additional information, the Kalman filter implemented with the true paramaters yields a log likelihood of 57964.658446252404. The MLE optimized set of parameters yields 58138.962169457925.

$\endgroup$

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.