I have three time series of same length, all containing magnitude measurements of the same event "A". But each time series is using a different method of measurement.

My goal is to merge the three time series into a single one, so that I can more easily find the "upward" and "downward" phases. (Like a smooth sinus-curve or binary data-set 1=up, 0=down.)

For instance, as you can see the "raw" measurement data is very noisy and looks like this for measurement system number 3: The raw data from measurement system number 3

To get a better understanding of the data I have used a Moving Average (window length 2000) to smooth the time series data for each measurement system 1 - 3, which yields the following figure:

Each time series smoothed with a moving average

Since they are measurements of the same event, and often have their peaks and valleys at similar times, I would like to merge the three time series into a single time series, with as little noise as possible. What methods should I try? I'm open to try anything!

I have tried using Fourier analysis (FFT in the SciPy package), but I cannot find any significant frequencies in any part of the data.

[edit for Whuber's comment] I unfortunately only have access to the MA smoothed data for all three, which has the following statistics (using DataFrame.describe() function):

  • MA_measurement 1: mean 0.991957, std 0.156941
  • MA_measurement 2: mean -0.000003, std 0.000016
  • MA_measurement 3: mean -0.000800, std 0.000856

And using DataFrame.corr() to get the correlations between the (smoothed) measurement systems:

$$\begin{array}{c|c|c|c|} & \text{MA_meas3} & \text{MA_meas2} & \text{MA_meas1} \\ \hline \text{MA_meas3} & 1.000000 & 0.337050 & 0.297922\\ \hline \text{MA_meas2} & 0.337050 & 1.000000 & 0.807282 \\ \hline \text{MA_meas1} & 0.297922 & 0.807282 & 1.000000 \\ \hline \end{array}$$

For the unsmoothed (i.e. "raw") data I only have data for time series number 2 and 3:

  • measurement 2: mean -0.000003, std 0.000747
  • measurement 3: mean -0.000812, std 0.022399

And they have a correlation of Corr(meas2,meas3)=0.027199.

[edit 2] I have been able to get hold of the the MA data shown in the second graph (i.e. 3 subplots), hope this can be of use!

MA data measurement 1 MA data measurement 2 MA data measurement 3

[edit 3] To elaborate on Matt F.'s comment: I hope to find around 33 upward and downward phases in total (peaks + troughs) in each MA measurement series (see edit 2) of event A. In theory it should be cyclical, i.e. up -> down -> up -> down... etc.

  • 2
    $\begingroup$ What can you tell us about the statistical properties of these measurement systems? What we most need to know are their (a) biases and (b) second moments. The latter refers to possible correlations as well as different variances. $\endgroup$
    – whuber
    Dec 29, 2021 at 21:20
  • 1
    $\begingroup$ @whuber I am unsure how to compute the bias, but have added mean, standard deviation and correlations. $\endgroup$
    – litmus
    Dec 29, 2021 at 22:03
  • 1
    $\begingroup$ @litmus That's fine! You just have to use a lop-sided or even one-sided set of points of which to take the median when you are close to the edges of your overall time window. $\endgroup$ Jan 6, 2022 at 14:31
  • 1
    $\begingroup$ A lowpass filter is certainly going to do something similar. However, in my experience, there are several drawbacks: 1. You must know a fair amount about lowpass filters in order to get good performance. 2. The result of a lowpass filter does not contain actually measured values anymore. 3. In my experience, the signal-to-noise ratio gain from a lowpass filter isn't even close to a median filter. So what are the disadvantages of the median filter? The primary one is that the resulting data will not be as smooth as what you typically get from a lowpass filter. $\endgroup$ Jan 6, 2022 at 19:09
  • 1
    $\begingroup$ About how many up-periods and down-periods do you expect? As extremes, there might be three periods (if the event is a sine wave whose period is the length of the series) or hundreds (if the event is as spiky as the unsmoothed measurements), and those would suggest different models and techniques. $\endgroup$
    – Matt F.
    Jan 7, 2022 at 13:57

3 Answers 3


I would apply PCA to three series, then assess how much variance is explained by the first principal component. If it is very high, then use the PC1 as your signal. The PCA transforms your three signals into three orthogonal linear combinations of signals in decreasing order of variance. So, your first principal component can be seen as an average level of three signals. If PCA works very well, then PC1 will explain the most variance in the data

You may need to mess with phase differences, depending on how the measurements were done. For instance, if it's three microphones picking up the source, then sound waves reach with a lag. This can easily be taken care of with coherence/phase shift analysis. Once you identify the phase shifts, it's still the same PCA but with appropriate lags

  • $\begingroup$ Could you please elaborate? I looked up PCA on wikipedia but not sure what you mean. Also, I have attached links to the data if it is easy for you to try/show? $\endgroup$
    – litmus
    Jan 7, 2022 at 17:05
  • $\begingroup$ take a look at example in a. different domain: clarusft.com/… - interest rates modeling in finance $\endgroup$
    – Aksakal
    Jan 7, 2022 at 17:08

We have only few details given here, but in general, the problem of estimating the underlying signal from noisy measurements seems to me like a Kalman filtering problem (or structural time series more generally).

In the following, I follow very closely this textbook: https://www.stat.pitt.edu/stoffer/tsa4/ , Chapter 6 and in the code I'm using the R package accompanying the textbook.

(Here is how Kalman filter works in picutres: https://www.bzarg.com/p/how-a-kalman-filter-works-in-pictures/ )

To illustrate how one would approach such a problem, conisder $p=3$ noisy measurements $y_t$ of the same underlying phenomenon $x_t$. We can write this as

$y_t = A x_t + \epsilon_t$

$x_{t+1} = \delta + \Phi x_t + \eta_t$

$\epsilon_t \sim N(0,R)$, $\eta_t \sim N(0,Q)$, $x_1 \sim N(\mu_0,\Sigma_1)$

For simplicity, let's say that in our case $\Phi = 1$ and $A =\begin{bmatrix} 1 \\ 1 \\ 1\end{bmatrix}$, $\delta=0$, and that $R = \begin{bmatrix} r_{11} & 0 & 0 \\ 0 & r_{22} & 0 \\ 0 & 0 & r_{33} \\\end{bmatrix}$ is diagonal matrix (meaning that the three measurements $y_{t,1}, y_{t,2}, y_{t,3}$ are uncorrelated). The four unknown parameters that we want to estimate are variance of $\eta$ ,$Q$, and variances of $\epsilon$, $r_{ii}$. The following code shows how to do setup the above model and estimate it.

We can estimate this parameters by maximum likelihood, meaning the we will look for the values of these parameters for which the model-implied series of observations $\hat{y}_t$ resembles the actual one as much as possible.

enter image description here

enter image description here

We start by creating some artificial data

# library accompanying the book


# for reproducibility

# dimensions of the problem in our case
N <- 400; # number of periods
p <- 3; # no. of observation series


# we assumed that this holds
A <- matrix(c(1, 1, 1), p, 1);
Phi <- 1;drift<-0;

# generate shocks
shocks_eps1 <- rnorm(N,mean=0,sd=1);
shocks_eps2 <- rnorm(N,mean=0,sd=2);
shocks_eps3 <- rnorm(N,mean=0,sd=3);
shocks_eta <- rnorm(N,mean=0,sd=.5);


# create artificial data
y_t <- matrix(0, N, p) # Stores observations
x_t <- matrix(0, N, 1) # stores unobserved true state
x_t[1, 1] <- 0;
y_t[1, ] <- A * x_t[1, 1] + shocks_eps[1, ];

for (i in 2:N) {

x_t[i, 1] <- drift + Phi*x_t[i - 1, 1] + shocks_eta[i - 1];
y_t[i, ] <- A * x_t[i, 1] + shocks_eps[i, ];


Next, let's estimate the four parameters as described above

# now setup the estimation procedure
# see the textbook Ch. 6 for details

input <- rep(1,N)
mu0 <- 0; Sigma0 <- 1;

A <- array(A, dim=c(p,1,N))

# Function to Calculate Likelihood
Linn = function(para){
cQ <- para[1] # Q
cR1 <- para[2] # 11 element of R
cR2 <- para[3] # 22 element of R
cR3 <- para[4] # 22 element of R
cR <- diag(c(cR1,cR2, cR3)) # put the matrix together
kf <- Kfilter1(N,y_t,A,mu0,Sigma0,Phi,0,drift,cQ,cR,input)

# Estimation
init.par <- c(.5,1,2,3) # initial values of parameters
(est <- optim(init.par, Linn, NULL, method='Nelder-Mead', hessian=TRUE,
    control=list(trace=1,REPORT=1))) # output not shown
(est <- optim(est$par, Linn, NULL, method='BFGS', hessian=TRUE,
    control=list(trace=1,REPORT=1))) # output not shown

cQ <- est$par[1]
cR1 <- est$par[2]
cR2 <- est$par[3]
cR3 <- est$par[4]
cR <- diag(c(cR1,cR2, cR3)) # put the matrix together

ks <- Ksmooth1(N,y_t,A,mu0,Sigma0,Phi,0,drift,cQ,cR,input)

# extract estimated signal
xsm <- ts(as.vector(ks$xs), start=1)

# and correspodning uncertainty
rmse <- ts(sqrt(as.vector(ks$Ps)), start=1)
  • $\begingroup$ I have tried to use KF but it does not seem to work on the data I have provided (found in question, edit 2). I think it becomes much easier with artificial data, such as the type you created. Unfortunately it is my experiance so far that it is different from actual measured data. Would you like to try yourself with the data I provided in my question? $\endgroup$
    – litmus
    Jan 8, 2022 at 13:35

As a base reference you could try to directly merge the three data sets by interlacing them and then analyzing this new data set.

First you need to do a linear transformation on each data set individually so they have the same mean and variance (mean=0 and variance=1 is the simplest choice). Then you transform the x-scales of the three data sets by $x \mapsto 3x+1, 2, 3$ respectively. If you now look at the union of the three data sets you have created one new data set merged out of the three individual ones.

Next I would apply the moving average again (with a window length of now 6000) and compare the numbers and locations of peaks and throughs with the results on the three individual data sets and see whether this gives anything useful.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.