Two large ensembles (LEs) of historical climate simulations are used to compare how various statistical methods estimate the sea surface temperature (SST) changes due to anthropogenic and other external forcing, and how their removal affects the internally generated Atlantic multidecadal oscillation (AMO), Pacific decadal oscillation (PDO), and the SST footprint of the Atlantic meridional overturning circulation (AMOC). Removing the forced SST signal by subtracting the global mean SST (GM) or a linear regression on it (REGR) leads to large errors in the Pacific. Multidimensional ensemble empirical mode decomposition (MEEMD) and quadratic detrending only efficiently remove the forced SST signal in one LE, and cannot separate the short-term response to volcanic eruptions from natural SST variations. Removing a linear trend works poorly. Two methods based on linear inverse modeling (LIM), one where the leading LIM mode represents the forced signal and another using an optimal perturbation filter (LIMopt), perform consistently well. However, the first two LIM modes are sometimes needed to represent the forced signal, so the more robust LIMopt is recommended. In both LEs, the natural AMO variability seems largely driven by the AMOC in the subpolar North Atlantic, but not in the subtropics and tropics, and the scatter in the AMOC–AMO correlation is large between individual ensemble members. In three observational SST reconstructions for 1900–2015, linear and quadratic detrending, MEEMD, and GM yield somewhat different AMO behavior, and REGR yields smaller PDO amplitudes. Based on LIMopt, only about 30% of the AMO variability is internally generated, as opposed to more than 90% for the PDO. The natural SST variability contribution to global warming hiatus is discussed.