Methods to extract signal from three very noisy time series of same event 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:

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:

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.
 A: 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
A: 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.




We start by creating some artificial data
# library accompanying the book

library("astsa")

# for reproducibility
set.seed(1);

# 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);

shocks_eps=cbind(shocks_eps1,shocks_eps2,shocks_eps3)

# 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)
return(kf$like)
}

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

A: 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.
