The GERDA experiment searches for the neutrinoless double beta decay of Ge-76 using high-purity germanium detectors enriched in Ge-76. The analysis of the signal time structure provides a powerful tool to identify neutrinoless double beta decay events and to discriminate them from gamma-ray induced backgrounds. Enhanced pulse shape discrimination capabilities of Broad Energy Germanium detectors with a small read-out electrode have been recently reported. This paper describes the full simulation of the response of such a detector, including the Monte Carlo modeling of radiation interaction and subsequent signal shape calculation. A pulse shape discrimination method based on the ratio between the maximum current signal amplitude and the event energy applied to the simulated data shows quantitative agreement with the experimental data acquired with calibration sources. The simulation has been used to study the survival probabilities of the decays which occur inside the detector volume and are difficult to assess experimentally. Such internal decay events are produced by the cosmogenic radio-isotopes Ge-68 and Co-60 and the neutrinoless double beta decay of Ge-76. Fixing the experimental acceptance of the double escape peak of the 2.614 MeV photon to 90%, the estimated survival probabilities at Qbb = 2.039 MeV are (86+-3)% for Ge-76 neutrinoless double beta decays, (4.5+-0.3)% for the Ge-68 daughter Ga-68, and (0.9+0.4-0.2)% for Co-60 decays.