Hysteresis is a nonlinear dissipative phenomenon present in many structures, such as those assembled by bolted joints. However, the approaches for identification are still limited in the literature due to their complexity. Additionally, there are also uncertainties held in bolted structures caused to fluctuations in tightening torque and pressure distribution along the contact surface. Thus, this paper yields a methodology for identifying a stochastic Bouc-Wen model for bolted joints based on the harmonic balance method. Unfortunately, some challenges are encountered when applying conventional series approximation to hysteresis, caused mainly by non-smooth behavior, which induces abrupt transitions between different motion regimes. In this work, previous adaptations were made to split the hysteresis loop in smooth paths and then use a piecewise harmonic balance approach. In this way, it was possible to deal with a deterministic Identification problem based on minimizing the error between the Fourier amplitudes of an experimental signal and those obtained through harmonic balance applying the Cross-Entropy optimization method. So, the results were extended to a stochastic model by applying the Bayesian paradigm, in which the maximized likelihood function was also based on the harmonic balance amplitudes. This methodology was demonstrated to identify Bouc-Wen parameters capable of predicting hysteresis in the BERT benchmark, composed of two aluminum beams jointed by a bolted joint in a cantilever boundary condition. Evaluating the results in the time and frequency domain and the nonlinear behavior through the hysteresis loop, it can be concluded that the method was able to identify an accurate stochastic Bouc-Wen model in predicting the dynamics of bolted structures even taking into account the probable uncertainties of the system.