Suppose we have a collection of bivariate random variables $X_{1i}$ and $X_{2i}$ indexed by a continuous variable $t$ such that, for the vector ${\bf{X}} = (X_{1i} \ X_{2i})^T$ we can assume
\begin{equation*} {\bf{X}} \sim {\bf{N(\mu}}(t_{i}), {\bf{\Sigma}}(t_{i})) \end{equation*}
for some unknown 2-dimensional vector ${\bf{\mu}}(t_{i})$ and a 2x2 covariance matrix ${\bf{\Sigma}}(t_{i})$, which must be estimated.
How can we estimate the form of the covariance matrix ${\bf{\Sigma}}(t)$ using independent observations $(X_{1i} \ X_{2i})^T$ measured at $t_1, ..., t_n, \ 1 \leq i \leq n$?
One approach I have considered is data binning, i.e. diving the range of values of $t$ into $m$ intervals and calculating $\bf{\Sigma}$ for data with $t$ in each interval $[a_j,b_j], \ 1 \leq j \leq m, \ 1 < m < n$. This is undesirable given the specific relationship I am dealing with but it could work.
The specific problem I am looking at is TLE covariance. The TLE is a set of numbers sent out by a satellite that can be used to approximate a satellite's position and velocity at a time at, before, or after the TLE is generated. Therefore, given a TLE and a time $t$ (measured from the TLE epoch, i.e. TLE generation time), the expected position of the satellite in 3D Euclidean space depends on $t$ as does the error ellipsis (which expands with increasing $t$). The position uncertainty (i.e., information about the error ellipsis) is not public but attempts have been made to determine this from public data in papers like this one:
H. Yurong, L. Zhi and H. Lei, "Covariance propagation of two-line element data," 2016 Chinese Control and Decision Conference (CCDC), Yinchuan, China, 2016, pp. 3836-3841, doi: 10.1109/CCDC.2016.7531654.
This paper specifically uses data binning and mentions using "the quadratic polynomial fitting error standard deviation of propagation time function in least squares sense" to estimate the error evolution function (i.e. how the error estimate of the estimated position varies with time). I'm interested in ways to generate a function $\Sigma(t)$ to capture this evolution with time and investigate the behaviour of error ellipses at different times. I have stated the problem in 2 dimensions and with a Gaussian PDF as I think this would be a good starting point for understanding the theory.