  We present a novel method for stochastic interpolation of sparsely sampled
time signals based on a superstatistical random process generated from a
multivariate Gaussian scale mixture. In comparison to other stochastic
interpolation methods such as Gaussian process regression, our method possesses
strong multifractal properties and is thus applicable to a broad range of
real-world time series, e.g. from solar wind or atmospheric turbulence.
Furthermore, we provide a sampling algorithm in terms of a mixing procedure
that consists of generating a 1 + 1-dimensional field u(t, {\xi}), where each
Gaussian component u{\xi}(t) is synthesized with identical underlying noise but
different covariance function C{\xi}(t,s) parameterized by a log-normally
distributed parameter {\xi}. Due to the Gaussianity of each component
u{\xi}(t), we can exploit standard sampling alogrithms such as Fourier or
wavelet methods and, most importantly, methods to constrain the process on the
sparse measurement points. The scale mixture u(t) is then initialized by
assigning each point in time t a {\xi}(t) and therefore a specific value from
u(t, {\xi}), where the time-dependent parameter {\xi}(t) follows a log-normal
process with a large correlation time scale compared to the correlation time of
u(t, {\xi}). We juxtapose Fourier and wavelet methods and show that a
multiwavelet-based hierarchical approximation of the interpolating paths, which
produce a sparse covariance structure, provide an adequate method to locally
interpolate large and sparse datasets.
