To investigate the Brownian motion of individual particles suspended in viscoelastic fluids, the stochastic smoothed profile method (SPm) for direct numerical simulation is developed by extending deterministic SPm for suspensions in viscoelastic media. To simulate viscoelastic flow driven by thermal fluctuations in the suspending medium, the random stress in the fluid momentum equation as well as the random driving force for the conformation tensor in the Oldroyd-B model are incorporated according to the fluctuating hydrodynamics and fluctuating viscoelasticity formalisms. The thermal equilibrium and dynamical properties calculated by using numerical simulations successfully reproduce the analytic predictions, validating the direct simulation for the coupled fluctuating Navier–Stokes and Oldroyd-B equations and for the coupling between the stochastic viscoelastic medium and individual particles. As an application of the stochastic SPm, we investigate finite system-size effects under periodic boundary conditions (PBCs) on the passive microrheological relationship between the mean-square displacement (MSD) of a Brownian particle and the medium's dynamic modulus. Comparing the modulus that was microrheologically calculated from the MSD with the input modulus reveals that the effect of periodic image cell interaction appears not only in the long-time diffusive regime but also in the short-time region. A frequency-dependent finite system-size correction is implemented by phenomenologically extending the long-time diffusive regime correction, allowing passive microrheology analysis under PBCs. This result can be directly applied to other mesoscale numerical simulations including coarse-grained molecular dynamics and dissipative particle dynamics simulations.