Testing nonlinear structures to characterise their internal nonlinear forces is challenging. Often nonlinear structures are excited by harmonic forces and yield a multi-harmonic response. In many systems, particularly ones with strong nonlinearities, the effect of higher harmonics in the force and responses cannot be ignored. Even if the intended excitation is a single frequency sinusoidal force, the interaction of the shaker and the nonlinear structure can lead to harmonics in the applied force. The effects of these higher harmonics of the input force on nonlinear model identification in structural dynamics are often neglected. The objective of this study is to introduce an identification method, motivated by the alternating frequency/time approach using harmonic balance (AFTHB), which is able to consider both multi-harmonic forces and multi-harmonic responses of the system. The proposed AFTHB method can include all significant harmonics by selecting an appropriate time step and sampling frequency to guarantee the accuracy of the results. An analytical harmonic-balance-based (AHB) approach is also considered for comparison. However, the inclusion of all significant harmonics of the response in the analytical expansion of the nonlinear functions is often cumbersome. Furthermore, the AFTHB method can easily cope with complex nonlinearities such as Coulomb friction and with multi-degree of freedom nonlinear systems. Including the effect of higher harmonics in the identification process reduces the approximation error due to truncation and very accurate approximation of the balanced equations of each harmonic is obtained. The proposed identification method requires prior knowledge or an appropriate estimation of the type of system nonlinearities. However, the method of model selection may be used for a set of candidate models, and avoiding a dictionary of arbitrary candidate basis functions significantly reduces the computational costs. This paper highlights the important features of the AFTHB method to ensure accurate estimation using four simulated and two experimental examples. The effects of the number of harmonics considered, the modelling error, measurement noise and the frequency range on the quality of the estimated model are demonstrated.