We deal with identification of Hammerstein-Wiener systems, or NLN systems, in which a linear subsystem is sandwiched by two nonlinearities. First we identify an approximate linear state space model of the NLN system by using the ORT (orthogonal projection) subspace method or PO-MOESP method. Then, initialized by the estimated state space model, we optimize the output error (OE) model, which is derived based on the basis functions expansion of nonlinearities, by using a gradient-based optimization method. Numerical results are included to show the applicability of the present approach.