Velocity correlation functions and third order diffusion coefficients of ions moving in a buffer gas under the influence of an electrostatic field are determined via molecular dynamics simulation. For the closed shell system of K+ in Ar using a universal interaction model potential, the general form of the third order correlation functions is found to be monotonically decaying in time except in the cases of [Delta upsilon(Z)(0)Delta upsilon(X)(t)(2)], [Delta upsilon(Z)(0)Delta upsilon(Y)(t)(2)], and [Delta upsilon(Z)(0)Delta upsilon(Z)(t)(2)], with Delta upsilon(t)=upsilon(t)-[upsilon(t)] and the field in the z direction. These functions acquire positive slope at short times showing enhancement of correlations between instantaneous upsilon(z) components of the ions and their future kinetic energies or velocity measures. This feature is shown to quantify the dynamics of correlations between velocity components suggested in the past by Ong, Hogan, Lam and Viehland [Phys. Rev. A 45, 3997 (1992)] in order to explain the form of an ion velocity distribution function calculated through a Monte Carlo simulation method. In addition, within a stochastic analysis which establishes a relation between velocity correlation functions and third order diffusion coefficients, only two independent components of the diffusion tensor, Q(parallel to) and Q(perpendicular to), are predicted. We thereby calculate the Q(perpendicular to), component, which has not been determined so far, over a wide field range. The magnitudes of the resulting third order diffusion coefficients indicate that their contribution to the ion transport in usual drift-tube measurements should be very small.