Laser shock peening (LSP) can be used to induce compressive residual stresses on the surface of a material, then to improve the mechanical properties such as performance of plasticity and fatigue. However, the residual stresses and their exact spatial distribution are very difficult to measure by experiment, especially for very small workpieces. In this paper, a finite-element model has been developed to numerically simulate the LSP process of bulk metallic glass (BMG) Zr41.2 Ti13.8Cu12.5Ni10Be22.5, and predict the stress distribution. The constitutive equation established in this work is hydrostatic-pressure sensitive and strain-rate dependent, it is based on the free volume model and Coulomb-Mohr yield criterion, and can describe such special deformation behaviors of BMG as strain softening. The simulated results show that, for one-side peening, along depth direction, the compressive residual stress gradually reduced to zero, then change to the tensile residual stress, but for two-side peening, the residual stress is from compressive to tensile and then to compressive along depth direction. These simulation results have a great significance to study the application of LSP in strengthening brittle amorphous alloys.