A central challenge in the computational modeling of neural dynamics is the\ntrade-off between accuracy and simplicity. At the level of individual neurons,\nnonlinear dynamics are both experimentally established and essential for\nneuronal functioning. An implicit assumption has thus formed that an accurate\ncomputational model of whole-brain dynamics must also be highly nonlinear,\nwhereas linear models may provide a first-order approximation. Here, we provide\na rigorous and data-driven investigation of this hypothesis at the level of\nwhole-brain blood-oxygen-level-dependent (BOLD) and macroscopic field potential\ndynamics by leveraging the theory of system identification. Using functional\nMRI (fMRI) and intracranial EEG (iEEG), we model the resting state activity of\n700 subjects in the Human Connectome Project (HCP) and 122 subjects from the\nRestoring Active Memory (RAM) project using state-of-the-art linear and\nnonlinear model families. We assess relative model fit using predictive power,\ncomputational complexity, and the extent of residual dynamics unexplained by\nthe model. Contrary to our expectations, linear auto-regressive models achieve\nthe best measures across all three metrics, eliminating the trade-off between\naccuracy and simplicity. To understand and explain this linearity, we highlight\nfour properties of macroscopic neurodynamics which can counteract or mask\nmicroscopic nonlinear dynamics: averaging over space, averaging over time,\nobservation noise, and limited data samples. Whereas the latter two are\ntechnological limitations and can improve in the future, the former two are\ninherent to aggregated macroscopic brain activity. Our results, together with\nthe unparalleled interpretability of linear models, can greatly facilitate our\nunderstanding of macroscopic neural dynamics and the principled design of\nmodel-based interventions for the treatment of neuropsychiatric disorders.\n