  We investigate the rotational evolution of young stars through Monte Carlo
simulations. We simulate 280,000 stars, each of which is assigned a mass, a
rotational period, and a mass accretion rate. The mass accretion rate depends
on mass and time, following power-laws indices 1.4 and -1.5, respectively. A
mass-dependent accretion threshold is defined below which a star is considered
as diskless, which results in a distribution of disk lifetimes that matches
observations. Stars are evolved at constant angular spin rate while accreting
and at constant angular momentum when they become diskless. We recover the
bimodal period distribution seen in several young clusters. The short period
peak consists mostly of diskless stars and the long period one is mainly
populated by accreting stars. Both distributions present a long tail towards
long periods and a population of slowly rotating diskless stars is observed at
all ages. We reproduce the observed correlations between disk fraction and spin
rate, as well as between IR excess and rotational period. The period-mass
relation we derive from the simulations exhibits the same global trend as
observed in young clusters only if we release the disk locking assumption for
the lowest mass stars. We find that the time evolution of median specific
angular momentum follows a power law index of -0.65 for accreting stars and of
-0.53 for diskless stars, a shallower slope that results from a wide
distribution of disk lifetimes. Using observationally-documented distributions
of disk lifetimes, mass accretion rates, and initial rotation periods, and
evolving an initial population from 1 to 12 Myr, we reproduce the main
characteristics of pre-main sequence angular momentum evolution, which supports
the disk locking hypothesis. (abridged)
