We employ a new NARMAX (Nonlinear Auto-Regressive Moving Average with eXogenous inputs) code to disentangle the time-varying relationship between the solar wind and SYM-H. The NARMAX method has previously been used to formulate a Dst model, using a preselected solar wind coupling function. In this work, which uses the higher-resolution SYM-H in place of Dst, we are able to reveal the individual components of different solar wind-magnetosphere interaction processes as they contribute to the geomagnetic disturbance. This is achieved with a graphics processing unit (GPU)-based NARMAX code that is around 10 orders of magnitude faster than previous efforts from 2005, before general-purpose programming on GPUs was possible. The algorithm includes a composite cost function, to minimize overfitting, and iterative reorthogonalization, which reduces computational errors in the most critical calculations by a factor of ∼106. The results show that negative deviations in SYM-H following a southward interplanetary magnetic field (IMF) are first a measure of the increased magnetic flux in the geomagnetic tail, observed with a delay of 20–30 min from the time the solar wind hits the bow shock. Terms with longer delays are found which represent the dipolarization of the magnetotail, the injections of particles into the ring current, and their subsequent loss by flowout through the dayside magnetopause. Our results indicate that the contribution of magnetopause currents to the storm time indices increase with solar wind electric field, E = v × B. This is in agreement with previous studies that have shown that the magnetopause is closer to the Earth when the IMF is in the tangential direction.