Summary. The moment method is a well-known astronomical mode identification technique in asteroseismology which uses a time series of the first three moments of a spectral line to estimate the discrete oscillation mode parameters l and m. The method, in contrast with many other mode identification techniques, also provides estimates of other important continuous parameters such as the inclination angle α and the rotational velocity ve. We developed a statistical formalism for the moment method based on so-called generalized estimating equations. This formalism allows an estimation of the uncertainty of the continuous parameters, taking into account that the different moments of a line profile are correlated and that the uncertainty of the observed moments also depends on the model parameters. Furthermore, we set up a procedure to take into account the mode uncertainty, i.e. the fact that often several modes (l, m) can adequately describe the data. We also introduce a new lack-of-fit function which works at least as well as a previous discriminant function, and which in addition allows us to identify the sign of the azimuthal order m. We applied our method to star HD181558 by using several numerical methods, from which we learned that numerically solving the estimating equations is an intensive task. We report on the numerical results, from which we gain insight in the statistical uncertainties of the physical parameters that are involved in the moment method.