We develop a rigorous asymptotic derivation of two mathematical models of water waves that capture the full nonlinearity of the Euler equations up to quadratic and cubic interactions, respectively. Specifically, letting denote an asymptotic parameter denoting the steepness of the water wave, we use a Stokes expansion in to derive a set of linear recursion relations for the tangential component of velocity, the stream function, and the water wave parameterization. The solution of the water wave system is obtained as an infinite sum of solutions to linear problems at each O(k) level, and truncation of this series leads to our two asymptotic models, which we call the quadratic and cubic h-models. These models are well posed in spaces of analytic functions. We prove error bounds for the difference between solutions of the h-models and the water wave system. We also show that the Craig-Sulem models of water waves can be obtained from our asymptotic procedure. We then develop a novel numerical algorithm to solve the quadratic and cubic h-models as well as the full water wave system. For three very different examples, we show that the agreement between the model equations and the water wave solution is excellent, even when the wave steepness is quite large. We also present a numerical example of corner formation for water waves.