In [1]:
import os, sys
sys.path.append(os.path.abspath('../../main/python'))
import datetime as dt
import numpy as np
import numpy.testing as npt
import matplotlib.pyplot as plt
from thalesians.tsa.distrs import NormalDistr as N
import thalesians.tsa.filtering as filtering
import thalesians.tsa.filtering.kalman as kalman
import thalesians.tsa.filtering.particle as particle
import thalesians.tsa.numpyutils as npu
import thalesians.tsa.processes as proc
import importlib
importlib.reload(particle)
importlib.reload(proc)
Out[1]:
In [7]:
process = proc.WienerProcess.create_from_cov(mean=3., cov=0.0001)
This we pass to a newly created particle filter, along with the initial time and initial state. The latter takes the form of a normal distribution. We have chosen to use Python datetimes as our data type for time, but we could have chosen ints or something else.
In [8]:
t0 = dt.datetime(2017, 5, 12, 16, 18, 25, 204000)
pf = particle.ParticleFilter(t0, state_distr=N(mean=100., cov=0.0000000000001), process=process)
Next we create an observable, which incorporates a particular observation model. In this case, the observation model is particularly simple, since we are observing the entire state of the particle filter. Our observation model is a 1x1 identity:
In [10]:
observable = pf.create_observable(kalman.LinearGaussianObsModel.create(1.), process)
Let's roll forward the time by one hour:
In [11]:
t1 = t0 + dt.timedelta(hours=1)
What is our predicted observation at this time? Since we haven't observed any actual information, this is our prior observation estimate:
In [12]:
prior_predicted_obs1 = observable.predict(t1)
prior_predicted_obs1
Out[12]:
We confirm that this is consistent with how our (linear-Gaussian) process model scales over time:
In [13]:
np.mean(pf._prior_particles), 100. + 3./24.
Out[13]:
In [14]:
prior_predicted_obs1
Out[14]:
In [15]:
prior_predicted_obs1 = observable.predict(t1)
npt.assert_almost_equal(prior_predicted_obs1.distr.mean, 100. + 3./24.)
npt.assert_almost_equal(prior_predicted_obs1.distr.cov, 250. + 25./24.)
npt.assert_almost_equal(prior_predicted_obs1.cross_cov, prior_predicted_obs1.distr.cov)
Let us now actually observe our observation. Say, the observation is 100.35 and the observation noise covariance is 100.0:
In [ ]:
observable.observe(time=t1, obs=N(mean=100.35, cov=100.0))
Having seen an actual observation, let us obtain the posterior observation estimate:
In [ ]:
posterior_predicted_obs1 = observable.predict(t1); posterior_predicted_obs1
We can now fast-forward the time, by two hours, say, and repeat the process:
In [ ]:
t2 = t1 + dt.timedelta(hours=2)
prior_predicted_obs2 = observable.predict(t2)
npt.assert_almost_equal(prior_predicted_obs2.distr.mean, 100.28590504 + 2.*3./24.)
npt.assert_almost_equal(prior_predicted_obs2.distr.cov, 71.513353115 + 2.*25./24.)
npt.assert_almost_equal(prior_predicted_obs2.cross_cov, prior_predicted_obs2.distr.cov)
observable.observe(time=t2, obs=N(mean=100.35, cov=100.0))
posterior_predicted_obs2 = observable.predict(t2)
npt.assert_almost_equal(posterior_predicted_obs2.distr.mean, 100.45709020)
npt.assert_almost_equal(posterior_predicted_obs2.distr.cov, 42.395213845)
npt.assert_almost_equal(posterior_predicted_obs2.cross_cov, posterior_predicted_obs2.distr.cov)
The real power of our particle filter interface is demonstrated for process models consisting of several (independent) stochastic processes:
In [ ]:
process1 = proc.WienerProcess.create_from_cov(mean=3., cov=25.)
process2 = proc.WienerProcess.create_from_cov(mean=[1., 4.], cov=[[36.0, -9.0], [-9.0, 25.0]])
Such models are common in finance, where, for example, the dynamics of a yield curve may be represented by a (multivariate) stochastic process, whereas the idiosyncratic spread for each bond may be an independent stochastic process.
Let us pass process1 and process2 as a (compound) process model to our particle filter, along with the initial time and state:
In [ ]:
t0 = dt.datetime(2017, 5, 12, 16, 18, 25, 204000)
kf = kalman.KalmanFilter(
t0,
state_distr=N(
mean=[100.0, 120.0, 130.0],
cov=[[250.0, 0.0, 0.0],
[0.0, 360.0, 0.0],
[0.0, 0.0, 250.0]]),
process=(process1, process2))
We shall now create several observables, each corresponding to a distinct observation model. The first one will observe the entire state:
In [ ]:
state_observable = kf.create_observable(
kalman.KalmanFilterObsModel.create(1.0, np.eye(2)),
process1, process2)
The second observable will observe the first coordinate of the first process:
In [ ]:
coord0_observable = kf.create_observable(
kalman.KalmanFilterObsModel.create(1.),
process1)
The third, the first coordinate of the second process:
In [ ]:
coord1_observable = kf.create_observable(
kalman.KalmanFilterObsModel.create(npu.row(1., 0.)),
process2)
The fourth, the second coordinate of the second process:
In [ ]:
coord2_observable = kf.create_observable(
kalman.KalmanFilterObsModel.create(npu.row(0., 1.)),
process2)
The fifth will observe the sum of the entire state (across the two processes):
In [ ]:
sum_observable = kf.create_observable(
kalman.KalmanFilterObsModel.create(npu.row(1., 1., 1.)),
process1, process2)
And the sixth a certain linear combination thereof:
In [ ]:
lin_comb_observable = kf.create_observable(
kalman.KalmanFilterObsModel.create(npu.row(2., 0., -3.)),
process1, process2)
Fast-forward the time by one hour:
In [ ]:
t1 = t0 + dt.timedelta(hours=1)
Let's predict the state at this time...
In [ ]:
predicted_obs1_prior = state_observable.predict(t1)
predicted_obs1_prior
And check that it is consistent with the scaling of the (multivariate) Wiener process with time:
In [ ]:
npt.assert_almost_equal(predicted_obs1_prior.distr.mean,
npu.col(100.0 + 3.0/24.0, 120.0 + 1.0/24.0, 130.0 + 4.0/24.0))
npt.assert_almost_equal(predicted_obs1_prior.distr.cov,
[[250.0 + 25.0/24.0, 0.0, 0.0],
[0.0, 360.0 + 36.0/24.0, -9.0/24.0],
[0.0, -9.0/24.0, 250 + 25.0/24.0]])
npt.assert_almost_equal(predicted_obs1_prior.cross_cov, predicted_obs1_prior.distr.cov)
Suppose that a new observation arrives, and we observe each of the three coordinates individually:
In [ ]:
state_observable.observe(time=t1, obs=N(mean=[100.35, 121.0, 135.0],
cov=[[100.0, 0.0, 0.0],
[0.0, 400.0, 0.0],
[0.0, 0.0, 100.0]]));
Let's look at our (posterior) predicted state:
In [ ]:
state_observable.predict(t1)
Let's also look at the predictions for the individual coordinates:
In [ ]:
coord0_observable.predict(t1)
In [ ]:
coord1_observable.predict(t1)
In [ ]:
coord2_observable.predict(t1)
The predicted sum:
In [ ]:
sum_observable.predict(t1)
And the predicted linear combination:
In [ ]:
lin_comb_observable.predict(t1)
Let's now go 30 minutes into the future:
In [ ]:
t2 = t1 + dt.timedelta(minutes=30)
And observe only the first coordinate of the second process, with a pretty high confidence:
In [ ]:
coord1_observable.observe(time=t2, obs=N(mean=125.25, cov=4.))
How does our predicted state change?
In [ ]:
state_observable.predict(t2)
Thirty minutes later...
In [ ]:
t3 = t2 + dt.timedelta(minutes=30)
We observe the sum of the three coordinates, rather than the individual coordinates:
In [ ]:
sum_observable.observe(time=t3, obs=N(mean=365.00, cov=9.))
How has our prediction of the state changed?
In [ ]:
state_observable.predict(t3)
And what is its predicted sum?
In [ ]:
sum_observable.predict(t3)