A Monte Carlo program which simulates the response of SiPMs is presented. Input to the program are the mean number and the time distribution of Geiger discharges from light, as well as the dark-count rate. For every primary Geiger discharge from light and dark counts in an event, correlated Geiger discharges due to prompt and delayed cross-talk and after-pulses are simulated, and a table of the amplitudes and times of all Geiger discharges in a specified time window generated. A number of different physics-based models and statistical treatments for the simulation of correlated Geiger discharges can be selected. These lists for many events together with different options for the pulse shapes of single Geiger discharges are used to simulate charge spectra, as measured by a current-integrating charge-to-digital converter, or current transients convolved with an electronics response function, as recorded by a digital oscilloscope. The program can be used to compare simulations with different assumptions to experimental data, and thus find out which models are most appropriate for a given SiPM, optimise the operating conditions and readout for a given application or test programs which are used to extract SiPM parameters from experimental data.