We propose a data-driven stochastic method that allows the simulation of a complex system's long-term evolution. Given a large amount of historical data on trajectories in a multi-dimensional phase space, our method simulates the time evolution of a system based on a random selection of partial trajectories in the data without detailed knowledge of the system dynamics. We apply this method to a large data set of time evolution of approximately one million business firms for a quarter century. Accordingly, from simulations starting from arbitrary initial conditions, we obtain a stationary distribution in three-dimensional log-size phase space, which satisfies the allometric scaling laws of three variables. Furthermore, universal distributions of fluctuation around the scaling relations are consistent with the empirical data.