Turbulent flows span a wide range of spatial and temporal scales, creating a major computational bottleneck. Large eddy simulations (LES) address this challenge through coarse graining, where the effects of unresolved motions appear as an unclosed subgrid-scale (SGS) term in the coarse-grained equations. Recent efforts in computational fluid dynamics have explored numerous data-driven approaches to model the SGS term, typically aiming to reproduce the full coarse-grained flow field, a high-dimensional and computationally demanding task. We propose an alternative approach, grounded in the observation that many practical applications of LES focus on a limited set of quantities of interest (QoIs)-such as average energy or enstrophy—rather than the full flow field. By focusing on these QoIs, we reformulated the SGS modeling task as a low-dimensional learning problem, significantly reducing the computational complexity while improving model interpretability. The key to our approach lies in representing unresolved dynamics using a minimal set of scalar time series, one for each QoI. We present the tau-orthogonal (TO) method for 3D turbulence, which captures QoI-state dependence and temporal correlations using regularized least-squares regression combined with a multivariate Gaussian residual model. This yields a simple yet effective stochastic time-series prediction model. The resulting stochastic time-series model accurately reproduces long-term QoI distributions, maintains robust performance across hyperparameter settings, and captures key flow features, such as kinetic energy spectra and coherent structures, despite being trained solely on QoI trajectories.