Fractional stochastic integro-differential equations arise in many problems in mechanics, finance, biology, medical, social sciences. In many situation, we can not solve these equations analytically, therefore we try to get an approximate solution for these equations. In this article, we apply operational matrices method based on orthonormal Bernstein polynomials (OBPs) to solve fractional stochastic integro-differential equations. For this aim, we introduce Bernstein polynomials and OBPs and explain that how an arbitrary function can be expanded by using these polynomials. We get operational ma- trix and stochastic operational matrix based on OBPs. By applying these matrices fractional stochastic integro-differential equations convert to a linear system of algebraic equations which can be solved by a appropriate numerical method. Also, convergence and error analysis of the proposed method is discussed and an upper error bound of the numerical method are given. Finally, Accuracy and efficiency of the present method are shown with two examples.