This paper introduces a new operator relevant to input–output analysis of flows in a statistically steady regime far from the steady base flow: the mean resolvent ${{\boldsymbol{\mathsf{R}}}}_0$. It is defined as the operator predicting, in the frequency domain, the mean linear response to forcing of the time-varying base flow. As such, it provides the statistically optimal linear time-invariant approximation of the input–output dynamics, which may be useful, for instance, in flow control applications. Theory is developed for the periodic case. The poles of the operator are shown to correspond to the Floquet exponents of the system, including purely imaginary poles at multiples of the fundamental frequency. In general, evaluating mean transfer functions from data requires averaging the response to many realizations of the same input. However, in the specific case of harmonic forcings, we show that the mean transfer functions may be identified without averaging: an observation referred to as ‘dynamic linearity’ in the literature (Dahan et al., J. Fluid Mech., vol. 704, 2012, pp. 360–387). For incompressible flows in the weakly unsteady limit, i.e. when amplification of perturbations by the unsteady part of the periodic Jacobian is small compared to amplification by the mean Jacobian, the mean resolvent ${{\boldsymbol{\mathsf{R}}}}_0$ is well-approximated by the well-known resolvent operator about the mean flow. Although the theory presented in this paper extends only to quasi-periodic flows, the definition of ${{\boldsymbol{\mathsf{R}}}}_0$ remains meaningful for flows with continuous or mixed spectra, including turbulent flows. Numerical evidence supports the close connection between the two resolvent operators in quasi-periodic, chaotic and stochastic two-dimensional incompressible flows.