Identification of bilinear systems subject to white inputs is studied. Assuming bilinearity we use the crosscovariances between the output and the higher order Hermite polynomials of the input for estimating the coefficients of the Hermite series expansion of the output. The parameters of the bilinear model are obtained via a balanced factorization of the appropriately constructed Hankel matrix. The skewness of the singular values of the estimated Hankel matrix is applied for testing the order selection for the bilinear model.