The dynamic behavior of isolated structures is strongly controlled by the force–deformation constitutive behavior of the isolators. Among the different types of existing isolation devices, High Damping Rubber Bearings (HDRBs) are commonly used in practice, which behavior is highly non-linear and difficult to model analytically. Consequently, this article proposes a simple, but sufficiently accurate, mathematical model for simulating the non-linear shear behavior of HDRBs under large deformations, and an estimation procedure for its parameter values using the geometrical features and mechanical characteristics of the device. First, we briefly describe the phenomena observed in the experimental test data, as well as other phenomena not observed within the range of experimental deformations. Then, the mathematical formulation is presented, which is based on the consideration of two components connected in parallel, a hyperelastic spring and a dissipative component. The governing equation for the former is derived from the expanded formulation of the Mooney–Rivlin model for isotropic hyperelastic materials, and the latter from a Bouc-Wen model with hardening. A novel model is included to account for stiffness degradation, including scragging and Mullins effects, which is developed from experimental data of 924 tested devices. The proposed model fits well the experimental test results of HDRBs with different geometric features and material properties. Based on the evolution laws for the different variables, the model can be successfully used in structural dynamic analysis. To facilitate model calibration, a statistical estimation procedure is proposed to reduce the 17 force–deformation constitutive model parameters of the isolator to 9 unknown parameters, which are computed from the geometric features of the device and mechanical characteristics of the rubber material. This makes the calibration of the force–deformation constitutive model parameters feasible. The estimation procedure successfully predicts the behavior of an average device within a batch of HDRBs, showing good agreement with two different experimental datasets.