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).

This signal, `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.