We can exclude integration by including additional equation in equations. Also we don't need to use ParametricNDSolve since there are random initial conditions. Actually we can reduce time for one order with the same option PrecisionGoal -> 2 (this is Automatic option for Method -> "MonteCarlo"). Initial code
It is often easier to use built-in functions instead of compiling and/or parallelizing. Here, a combination of ParametricNDSolve and concurrent integration could work.
As an example, assume you are starting from a differential equation
f''(t) + [a + b \cos(t)] f(t) = 0
which depends on two parameters $(a,b)$. We can solve this equation numerically for ...