I’m doing research on multiple-model state estimators for a type of (discrete-time) dynamical system with a switching disturbance model. However, the estimators require estimates of the model parameters and I’m wondering if I can estimate these from time series data using a numerical method in PyMC. But I’m not sure which method or how to set up the problem in PyMC.
The system consists of a linear dynamic system with the disturbance at the input. The disturbance is modelled as an integrator with a special random variable at the input known as a ‘random shock’ which is sampled from either of two normal distributions depending on the system mode (a discrete random state of the system).
So in pseudo-code, the data generating model, which has three parameters, would be something like this:
def random_shock(sigma1, sigma2, epsilon): r = sample(dist.Bernoulli(epsilon)) if r == 0: wp = sample(dist.Normal(0, sigma1) else: wp = sample(dist.Normal(0, sigma2) return w_p
(It is a Gaussian mixture model with two modes).
w_p, is then fed through an integrator to produce the disturbance signal p at time k:
p(k) = p(k-1) + w_p(k)
p(k) is a ‘randomly-occurring step change’ process something like this:
This is then fed through a model of the process which for now assume is known. For simplicity, let’s assume the process model is a simple 1st-order stable system:
y(k) = 0.7 * y(k-1) + p(k)
The catch is, there is measurement error (noise) added to the output. So the measurement model is:
y_M(k) = y(k) + v(k)
where v(k) is the measurement noise.
Here is a diagram of the whole system:
My question is, if I had a time-series of measurements of
y_M(k) of length N, how could I estimate the random shock parameters, sigma1, sigma2, and epsilon?
I’m pretty sure this is a case of joint state and parameter estimation. I.e. We need to estimate r(k), k = 0, 1, …, N and sigma1, sigma2, and epsilon at the same time.