# Kalman Smoother with Repeated Measurements

I have a problem that I thought should be easy but is turning out to be quite difficult.

Lets say I have 100 measurements of a time-varying processes at time points $t\in {1, \dots, 100}$ which are denoted $y = (y_1, \dots, y_{100})$. However, I also have 20 repeated measurements of the last time-point such that I will say that $(y_{101}, \dots, y_{120})$ are all taken at time point $t = 100$.

To summarize, I have 1 measurement from each time point $t \in {1, \dots, 99}$ and I have 21 measurements from time-point $t=100$.

I believe that the kalman filter can be used for both time-varying and non-time-varying processes. What I would like to do is figure out how to performing filtering and smoothing of this hypothetical time-series using these replicate measurements to improve my state estimation. In particular I am interested in Linear-Gaussian models of the form

$$y_t = F_t'\Theta_t+v_t$$ $$\Theta_t = G_t\Theta_{t-1}+w_t$$ $$v_t \sim N(0, V)$$ $$w_t \sim N(0, W)$$

My current Approach

My current approach is to simply pretend that I have a time-series with 120 observations $(y_1, \dots, y_{120})$ but that nothing is changing between the last 20 time-points. Specifically I set:

1. all $G_t$ matrices in $G_{101}, \dots, G_{120}$ are the identity matrix with dimension equal to the dimension of my parameter space $\Theta_t$.

2. I set $W = \mathbf{0}$ for $t \in \{101, \dots, 120\}$

This seems to work fine for the kalman filter (at least I think it is correct, I would appreciate it if someone confirmed that this makes sense); but it breaks down for my smoothing relations.

Specifically I am using following the filtering/smoothing notation of the Prado and West book. My smoothing Recursions are as follows:

$$\Theta_T \sim p(\Theta_T|D_T)$$ for $t = T-1, T-2, \dots$
$$B_t = C_tG'_{t+1}R^{-1}_{t+1}$$ $$\Theta_t \sim N(u_t, H_t)$$ where $$u_t = m_t+B_t(\theta_{t+1}-a_{t+1})$$ $$H_t = C_t-B_tR_{t+1}B_t'$$

Note that $a_{t+1} = G_tm_t$ and $R_{t} = G_tC_{t-1}G_t'+W_t$ are both calculate in the forward filtering recursions (which I have not written down here).

The issue I am having is that with my approach that I described above, $B_t = 1$ (take $t=120$ and work backwards: $B_{120} = C_t\mathbf{1}R_t^{-1} = C_t(\mathbf{1}C_t^{-1}\mathbf{1}'+\mathbf{0}) = 1$ and I end up getting that $H_t = \mathbf{0}$ (which then throws an error in my code because you can't draw from a normal distribution with zero variance)!

What I am doing wrong and is there a better way to handle this? (Note that I only want 1 sample from all the "time-points" in $t\in\{100,\dots, 120\}$, I don't need 21 if that helps)

I don't have the Prado & West book handy. I think, in your notation, the recurrence equations should be:

$$u_t = m_t + B_t( u_{t+1} - a_{t+1})$$ $$H_t = C_t + B_t( H_{t+1} - R_{t+1})B_t'$$

Using different notation for the retrospective distribution's mean and variance parameters, if you have

• West, M., & Harrison, J. (1997). Bayesian forecasting and dynamic models (2nd ed.). New York: Springer.

see Ch. 4.7 Filtering recurrences, Theorem 4.4 proof, page 115.

Otherwise, see pages 102-103 of

• Bayesian risk management: a guide to model risk and sequential learning in financial markets by Matt Sekerke, John Wiley & Sons, 2015.

Here, use notation $a_t,R_t$ for the (forward) filtering prior distribution parameters; and, $u_t,H_t$ for the (backward) smoothed state estimates.

Store $\{a_t, R_t,m_t, C_t\} \quad\forall~ t=1 \ldots T$ for use below.

For $t=T$, $$u_T = m_T \quad\text{and}\quad H_T = C_T\quad.$$

For $t=T-1, ~\ldots,~ 1$ $$u_t = m_t + B_t(u_{t+1} - a_{t+1})\quad\text{and}$$ $$H_t = C_t + B_t(H_{t+1} - R_{t+1}) B_t^{'}\quad.$$

The smoothed state distribution is $\Theta_t \sim N(u_t,H_t) \quad.$

• I don't see how this deals with repeated observations. In fact I believe these are the exact recurrence relationships that have been causing me problems (please reread my question). Am I incorrect?
– jds
Aug 4, 2017 at 22:02
• Is Prado&West 2010 the copy you are referring to? Your equations are not consistent with Eqns. (4.10) and (4.11). P&W is consistent with the equations I provided, substituting notation $u_t = a_T(t-T)$, and $H_t = R_T(t-T)$, and of course changing +B(...) to -B(...) and reversing the signs inside the parenthesis. In your equation for the mean $u_t$, you don't know $\theta_{t+1}$ (and never will); so $u_{t+1}$ is the best available estimate you have. In your equation for the variance, you are subtracting the whole $R_{t+1}$, not a small correction $H_{t+1} - R_{t+1}$ in my equation. Aug 5, 2017 at 19:48
• For repeated observations, just be careful on calculating $R_t$. $R_{100}$ is as usual. For $R_{101}$ to $R_{120}$, $R_t = C_{t-1}$; but note $C_{100} > C_{101} > C_{102} \ldots > C_{120}$ as you learn more / perform Bayesian updates with the observations; and, importantly, don't incur state evolution $w$. I am assuming you apply the observations 100 to 120 in order (which doesn't mean anything, since they occur at the same time.) Otherwise, the posteriors $C_t$ will contract in the order you apply your updates. Aug 5, 2017 at 19:55
• I am dealing with forward-filtering backwards sampling recursions. Thats how I know $\theta_{t+1}$.
– jds
Aug 6, 2017 at 22:24

All parameters can be time varying. This includes the variance of the observation error $v_t$. If you instead let $y_{100}$ denote the average of what is observations 100 to 120 in your current notation, then the observation error at time $t=100$ has variance $\text{Var}(v_{100})=V_{100}=V/21$ while all other observation errors have variances $V_t=V$. You then plug this time varying variance into the appropriate places in the forecast, filtering and smoothing recursions.

• I think I see why you are saying this but I am not yet convinced that this is correct. Could you please elaborate on why this would work? How does this compare to the approach I am taking by sequentially processing the observations?
– jds
Aug 2, 2017 at 17:34
• @jds If you go through the proofs of the forecast, filtering and smoothing recursions you will see that there is nothing that precludes any of the parameters from being time dependent, a point also made by e.g. Shumway & Stoffer, 2010, corollary 6.1. Aug 2, 2017 at 17:52
• @jds I guess your approach may also work but you would need to make the process noise "time"-dependent, that is, when forecasting $\Theta_{101}$ etc. you should not add anything to the variance from the filtering step after having conditioned on $y_{100}$ and so on. So it seems that you just end up having something else "time" dependent. Aug 2, 2017 at 17:56
• I would like to award you the bounty but I must say you last comment was somewhat unintelligible. Could you please make this more clear with more details?
– jds
Aug 4, 2017 at 22:03
• @jds I don't have access to Prado & West but your equation $H_t = C_t - B_t R_{t+1} B_t'$ should read $H_t = C_t - B_t (R_{t+1}-H_{t+1} ) B_t'$ instead as noted by krkeane. Fixing this bug, your currrent approach might work too. Aug 7, 2017 at 13:37