Natural and social systems exhibit complex behavior reflecting their rich dynamics, whose governing laws are not fully known. This study develops a unified data-driven approach to estimate predictability of such systems when several independent realizations of the system's evolution are available. If the underlying dynamics are quasi-linear, the signal associated with the variable external factors, or forcings, can be estimated as the ensemble mean; this estimation can be optimized by filtering out the part of the variability with a low ensemble-mean-signal-to-residual-noise ratio. The dynamics of the residual internal variability is then encapsulated in an optimal, in a Bayesian sense, linear stochastic model able to predict the observed behavior. This model's self-forecast covariance matrices define a basis of patterns (directions) associated with the maximum forecast skill. Projecting the observed evolution onto these patterns produces the corresponding component time series. These ideas are illustrated by applying the proposed analysis technique to (1) ensemble data of regional sea-surface temperature evolution in the tropical Pacific generated by a state-of-the-art climate model and (2) consumer-spending records across multiple regions of the Russian Federation. These examples map out a range of possible solutions-from a solution characterized by a low-dimensional forced signal and a rich spectrum of predictable internal modes (1)-to the one in which the forced signal is extremely complex, but the number of predictable internal modes is limited (2). In each case, the proposed decompositions offer clues into the underlying dynamical processes, underscoring the usefulness of the proposed framework.