New approaches for numerical implementation of the path integration (PI) method are described. In essence the PI method is a stepwise calculation of the joint probability density function (PDF) of a set of state space variables describing a white noise excited nonlinear dynamic system. The basic idea behind the proposed procedure is to apply a splines interpolation method to the logarithm of the calculated PDF to obtain an accurate representation of the PDF over the whole domain and not only at the chosen grid points. This exploits the fact that the logarithm of the PDF shows a more polynomial behaviour than the PDF itself, and therefore is better adapted to a splines representation. It is demonstrated that the proposed techniques may lead to significantly improved performance in calculating the response statistics of large classes of nonlinear oscillators excited by white or coloured noise when compared to other available implementations of the PI method. An advantage of the new approaches is that they allow time-variant dynamic systems to be analysed without significant increase in computer time. Numerical results for both 2D and 3D problems are presented.