The movable antenna (MA) technology has attracted great attention recently due to its promising capability in improving wireless channel conditions by flexibly adjusting antenna positions. To reap maximal performance gains of MA systems, existing works mainly focus on MA position optimization to cater to the instantaneous channel state information (CSI). However, the resulting real-time antenna movement may face challenges in practical implementation due to the additional time overhead and energy consumption required, especially in fast time-varying channel scenarios. To address this issue, we propose in this paper a new approach to optimize the MA positions based on the users' statistical CSI over a large timescale. In particular, we propose a general field response based statistical channel model to characterize the random channel variations caused by the local movement of users. Based on this model, a two-timescale optimization problem is formulated to maximize the ergodic sum rate of multiple users, where the precoding matrix and the positions of MAs at the base station (BS) are optimized based on the instantaneous and statistical CSI, respectively. To solve this non-convex optimization problem, a log-barrier penalized gradient ascent algorithm is developed to optimize the MA positions, where two methods are proposed to approximate the ergodic sum rate and its gradients with different complexities. Finally, we present simulation results to evaluate the performance of the proposed design and algorithms based on practical channels generated by ray-tracing. The results verify the performance advantages of MA systems compared to their fixed-position antenna (FPA) counterparts in terms of long-term rate improvement, especially for scenarios with more diverse channel power distributions in the angular domain.