Abstract: Naturally occurring systems often have inherent uncertainties and often evolve on complicated smooth and differentiable hypersurfaces that are not necessarily Euclidean. The dynamics of such systems can be potentially described using a stochastic (opens in a new window)Hamiltonian formulation on differentiable manifold. Solution of these stochastic (opens in a new window)Hamiltonian systems using the current state-of-the-art methods, either fails to preserve the geometry of the manifold, or, ignore the stochastic nature of the differential equations. To address this, we develop a framework which preserves geometry while accounting for (opens in a new window)stochasticity in input and states. A new development for higher order geometry preserving Ito-Taylor expansion based stochastic integration scheme is demonstrated. Detailed mathematical development of the scheme is laid down and numerically demonstrated through the application of this method on a number of practical physical systems such as spherical pendulum on a cart and two quadrotors transporting a mass point. The accuracy of the method is established by comparing the global error of the solution with the existing geometric Euler-Maruyama integration scheme, geometric Milstein scheme and the non-geometric Ito-Taylor based stochastic integration schemes.