The variability of wind power production poses the greatest challenge in the integration of large scale wind power in power systems. Furthermore, larger scale penetration implies a wider geographical spreading of installed wind power, resulting in reduced variability and the smoothing effect of total power generation. Therefore, analysis of the impact of wind power variations on power system operation requires adequate modeling of aggregate power output from geographically dispersed wind farms. This paper analyzes different aspects of Markov chain Monte Carlo simulation methods for the synthetic generation of dependent wind power time series. However, testing indicates that these approaches do not adequately model the stochastic dependence between wind power time series in conjunction with individual persistence which is necessary to obtain realistic distributions of aggregate power output and total power variations. Consequently, a novel approach based on a modified second order Markov chain Monte Carlo simulation is proposed. Simulation results show that this method obtains synthetic time series of aggregate wind power which very closely fit the original data, with respect to both the cumulative density function of total output power and the probability density function of power variations.