  This paper has two interrelated foci: (i) obtaining stable and efficient
data-driven closure models by using a multivariate time series of partial
observations from a large-dimensional system; and (ii) comparing these closure
models with the optimal closures predicted by the Mori-Zwanzig (MZ) formalism
of statistical physics. Multilayer stochastic models (MSMs) are introduced as
both a generalization and a time-continuous limit of existing multilevel,
regression-based approaches to closure in a data-driven setting; these
approaches include empirical model reduction (EMR), as well as more recent
multi-layer modeling. It is shown that the multilayer structure of MSMs can
provide a natural Markov approximation to the generalized Langevin equation
(GLE) of the MZ formalism.
  A simple correlation-based stopping criterion for an EMR-MSM model is derived
to assess how well it approximates the GLE solution. Sufficient conditions are
derived on the structure of the nonlinear cross-interactions between the
constitutive layers of a given MSM to guarantee the existence of a global
random attractor. This existence ensures that no blow-up can occur for a very
broad class of MSM applications, a class that includes non-polynomial
predictors and nonlinearities that do not necessarily preserve quadratic energy
invariants.
  The approach is applied to conceptual nonlinear models borrowed from climate
dynamics and population dynamics. In both cases, it is shown that the resulting
closure models are able to capture the main statistical features of the
dynamics, even in presence of weak time-scale separation.
