Extracting Average of Several Timeseries

Example of taking a bunch of equal length timeseries and finding the “average” of them by treating each timeseries as a single draw from a multivariate normal distribution. The “fitted” mean vector of this distribution can be seen as the average.

library(mclust)
library(MASS)
library(dplyr)
library(ggplot2)
library(reshape2)

# Create a multivariate normal dataset with 25 timeseries, each having 100 values
ni <- 25; nt <- 100
d <- mvrnorm(
  ni,                        # Specify that we want ni seperate samples drawn
  diffinv(rnorm(nt-1)),       # Create a nonstationary mean vector
  cov_autocorrelation(nt, .7) # Create an auto-correlated covariance matrix (function omitted)
)

# Fit the multivariate normal distribution created randomly above
avg.est <- mvn('Spherical', d)$parameters$mean %>% melt %>% setNames(c('Time', 'Sample', 'Temp')) 

# Plot the random timeseries and their fitted "average"
d %>% melt %>% setNames(c('Sample', 'Time', 'Temp')) %>%
  ggplot(aes(x=Time, y=Temp, color=factor(Sample))) + 
  geom_line(alpha=.2) + geom_line(data=avg.est, color='black', width=3) + guides(color=F) +
  theme_bw() + ggtitle('Temp vs Time\n(Timeseries are light colored lines, fitted mean is black)')