We explore the structure and statistics of multiphase, magnetized ISM turbulence in the local Milky Way by means of driven periodic box numerical MHD simulations. Using the higher order-accurate piecewise-parabolic method on a local stencil (PPML), we carry out a small parameter survey varying the mean magnetic field strength and density while fixing the rms velocity to observed values. We quantify numerous characteristics of the transient and steady-state turbulence, including its thermodynamics and phase structure, kinetic and magnetic energy power spectra, structure functions, and distribution functions of density, column density, pressure, and magnetic field strength. The simulations reproduce many observables of the local ISM, including molecular clouds, such as the ratio of turbulent to mean magnetic field at 100 pc scale, the mass and volume fractions of thermally stable HI, the lognormal distribution of column densities, the mass-weighted distribution of thermal pressure, and the linewidth-size relationship for molecular clouds. Our models predict the shape of magnetic field probability density functions (PDFs), which are strongly non-Gaussian, and the relative alignment of magnetic field and density structures. Finally, our models show how the observed low rates of star formation per free-fall time are controlled by the multiphase thermodynamics and largescale turbulence.To explore this complexity, we develop self-consistent models based on very idealized numerical experiments, linking together scales relevant to molecular cloud formation and fragmentation. We exploit the concept of self-organization in nonequilibrium nonlinear dissipative systems [14] in application to the ISM, which is long overdue [15]. We treat interstellar turbulence as an agent that imposes 'order' in the form of coherent structures and correlations between flow fields emerging in a simple periodic box simulation when a statistically stationary state fully develops. In this case, the details of initial conditions are no longer important. Instead the steady state provides realistic initial conditions for star formation. Unlike various flavors of the popular converging flow scenario [16][17][18][19], this approach allows us to reproduce basic ISM observables, including properties of local molecular clouds, with only a few control parameters. The model can be readily extended to include more physics (e.g., radiative transfer and chemistry), to augment the range of resolved scales, and to close the star formation feedback loop.The paper is organized as follows. Section 2 contains the details of our numerical models, including equations, initial and boundary conditions, and parameters. Section 3 presents the results of statistical analysis and comparison with observations for a number of ISM diagnostics. Section 3.1 describes the transition to turbulence and global properties of the statistically stationary state. Section 3.2 summarizes the evolution of magnetic field under the action of random large-scale external forcing. In section 3...