A method of studying nonlinear random fields of surface gravity waves by direct numerical simulation is developed by combining the “high-order spectral method” of West et al. (J. Geophys. Res. 92 (1987) 11,803–11,824) and the complex amplitude function b(k) introduced by Zakharov (J. Appl. Mech. Tech. Phys. (Engl. Transl.) 2 (1968) 190–194) in his Hamiltonian formalism of the water wave problem. By computing the rate of change of b(k) for a given wave field by both the high-order spectral method and Zakharov's equation, it is confirmed that the high-order spectral method, when the order of nonlinearity is appropriately chosen, gives exactly the same result as Zakharov's equation within a much shorter computation time. Methods of evaluating the frequency spectrum, directional spectrum and the component phase speed are given. Results of some preliminary computations are also presented to prove the great potentialities of the present approach as a tool for the examination of various characteristics of nonlinear random wave fields with continuous spectra, such as ocean waves.