The production of artificial light curves with known statistical and variability properties is of great importance in astrophysics. Consolidating the confidence levels during cross-correlation studies, understanding the artefacts induced by sampling irregularities, establishing detection limits for future observatories are just some of the applications of simulated data sets. Currently, the widely used methodology of amplitude and phase randomization is able to produce artificial light curves which have a given underlying power spectral density (PSD) but which are strictly Gaussian distributed. This restriction is a significant limitation, since the majority of the light curves, e.g. active galactic nuclei, X-ray binaries, gamma-ray bursts, show strong deviations from Gaussianity exhibiting ‘burst-like’ events in their light curves yielding longtailed probability density functions (PDFs). In this study, we propose a simple method which is able to precisely reproduce light curves which match both the PSD and the PDF of either an observed light curve or a theoretical model. The PDF can be representative of either the parent distribution or the actual distribution of the observed data, depending on the study to be conducted for a given source. The final artificial light curves contain all of the statistical and variability properties of the observed source or theoretical model, i.e. the same PDF and PSD, respectively. Within the framework of Reproducible Research, the code and the illustrative example used in this paper are both made publicly available in the form of an interactive MATHEMATICA notebook.