We investigate direct estimation of step-response models from closed-loop data using subspace identification. Necessary information concerning impulse-response coefficients is embedded in subspace matrices. Therefore, the step-response coefficients can be directly obtained from this matrix by integration without the need of extracting state space models first, as the conventional subspace identification does. Since the estimated subspace matrix contains more than one set of impulse-response coefficients, a question arises about how to best synthesize them to obtain an optimal estimate of the impulse-response coefficients and subsequently the step-response coefficients. For this purpose, a reformulation of the subspace identification problem is required for which the variance of all elements in the related subspace matrix can be evaluated. The calculated variances are then used to perform a weighted averaging on the estimated impulse-response coefficients to attenuate the noise influence on the final step-response model estimation. Monte Carlo simulations and pilot-scale experiments are provided to illustrate the proposed method.