Ground-based magnetic observatory series are the main source of information for constructing time-dependent spherical harmonic geomagnetic field models from sub-annual to pluri-decadal time scales. Assessing the reliability of such models requires accurate estimation of the data errors. We propose an analysis of observatory monthly means over the period 1930 to 2010, where we sequentially isolate (i) a stochastic regression for the main field at every site, performed in the framework of Gaussian processes, (ii) a local fit to annual and semiannual signals, (iii) a month by month estimate of global, large length-scale external and induced fields. We then estimate the unmodeled signal level (UMSL, which refers to the instrumental noise plus extra signals not captured by the above data treatment) from the standard deviation of the residuals to the sequential analysis. This may be used to estimate data error covariances in future field modeling studies. Mainly a function of the geomagnetic latitude, the UMSL is larger towards auroral regions and carries the temporal signature of solar activity. While the UMSL shows rather similar magnitudes in all three components in recent epochs (typically a few nT), a significant decrease is found in the downward component of the field around 1960, which correlates with the introduction of proton magnetometers. We detail the geographic distribution of the periodic signals and confirm the variation of their amplitude at pluri-decadal time scales. From the spherical harmonic description of horizontal and vertical fields, we isolate the main patterns of the inducing field in Z. These are dominated by a zonal structure of degree 1 (and to a lesser extent, of degree 3) in dipole coordinates. We nevertheless isolate secondary, non-zonal sources that are most active during the 1960s and around 1990, periods of particularly large solar activity, denoting an unusual morphology of the inducing field.