A kinetic approach for the evolution of ultracold neutral plasmas including interionic correlations and the treatment of ionization/excitation and recombination/deexcitation by rate equations is described in detail. To assess the reliability of the approximations inherent in the kinetic model, we have developed a hybrid molecular dynamics method. Comparison of the results reveals that the kinetic model describes the atomic and ionic observables of the ultracold plasma surprisingly well, confirming our earlier findings concerning the role of ion-ion correlations [Phys. Rev. A {bf 68}, 010703]. In addition, the molecular dynamics approach allows one to study the relaxation of the ionic plasma component towards thermodynamical equilibrium.