In this paper, we propose a promising extension of the path integral molecular dynamics method to Bose-Einstein and Fermi-Dirac statistics. The partition function for the quantum statistics was rewritten in a form amenable to the molecular dynamics method with the aid of an idea of pseudopotential for the permutation of particles. Our pseudopotential, here, is a rigorous one describing the whole effect of Bose-Einstein and Fermi-Dirac statistics. For a model calculation, we chose a system consisting of three independent particles in a one-dimensional harmonic well. The calculation has been performed for the particles obeying Bose-Einstein and Fermi-Dirac statistics. The calculated kinetic and potential energies were in excellent agreement with the analytical results even near the ground state. It was found that the pseudopotential shows attractive and repulsive characters for the static properties of Bose-Einstein and Fermi-Dirac particles, respectively. For interacting model particle systems, we studied a bosonic triatomic cluster. The calculated thermodynamic quantities were in qualitative agreement with those obtained by Fourier path integral Monte Carlo calculation.