Block-oriented nonlinear models are popular in nonlinear modeling because of their advantages to be quite simple to understand and easy to use. To increase the flexibility of single branch block-oriented models, such as Hammerstein, Wiener, and Wiener-Hammerstein models, parallel block-oriented models can be considered. This paper presents a method to identify parallel Wiener-Hammerstein systems starting from input-output data only. In the first step, the best linear approximation is estimated for different input excitation levels. In the second step, the dynamics are decomposed over a number of parallel orthogonal branches. Next, the dynamics of each branch are partitioned into a linear time invariant subsystem at the input and a linear time invariant subsystem at the output. This is repeated for each branch of the model. The static nonlinear part of the model is also estimated during this step. The consistency of the proposed initialization-procedure is proven. The method is validated on real-world measurements using a custom built parallel Wiener-Hammerstein test system. (C) 2014 Elsevier Ltd. All rights reserved.