diff --git a/doc/developers.md b/doc/developers.md
index 27bc3b145d..05e8b0af52 100644
--- a/doc/developers.md
+++ b/doc/developers.md
@@ -59,6 +59,9 @@ The following contains some general guidelines for developers.
   #include ...
   ```
 
+### File structure
+- Test files in `src/test/` end with the suffix `*Test.cpp`.
+
 ### Output
 - We provide custom macros for output and logging.
   The use of `std::cout` should be avoided and instead, macros such as `STORM_LOG_DEBUG`, `STORM_LOG_INFO` or `STORM_PRINT_AND_LOG` should be used.
diff --git a/src/storm-dft/simulator/DFTTraceSimulator.cpp b/src/storm-dft/simulator/DFTTraceSimulator.cpp
index f5cbb34015..7c6f94dd88 100644
--- a/src/storm-dft/simulator/DFTTraceSimulator.cpp
+++ b/src/storm-dft/simulator/DFTTraceSimulator.cpp
@@ -5,15 +5,16 @@ namespace simulator {
 
 template<typename ValueType>
 DFTTraceSimulator<ValueType>::DFTTraceSimulator(storm::dft::storage::DFT<ValueType> const& dft,
-                                                storm::dft::storage::DFTStateGenerationInfo const& stateGenerationInfo, boost::mt19937& randomGenerator)
+                                                storm::dft::storage::DFTStateGenerationInfo const& stateGenerationInfo,
+                                                storm::utility::RandomProbabilityGenerator<ValueType> randomGenerator)
     : dft(dft), stateGenerationInfo(stateGenerationInfo), generator(dft, stateGenerationInfo), randomGenerator(randomGenerator) {
     // Set initial state
     resetToInitial();
 }
 
 template<typename ValueType>
-void DFTTraceSimulator<ValueType>::setRandomNumberGenerator(boost::mt19937& randomNumberGenerator) {
-    this->randomGenerator = randomNumberGenerator;
+void DFTTraceSimulator<ValueType>::setRandomGenerator(storm::utility::RandomProbabilityGenerator<ValueType> randomGenerator) {
+    this->randomGenerator = randomGenerator;
 }
 
 template<typename ValueType>
@@ -65,8 +66,7 @@ std::tuple<storm::dft::storage::FailableElements::const_iterator, double, bool>
         bool successful = true;
         if (!dependency->isFDEP()) {
             // Flip a coin whether the PDEP is successful
-            storm::utility::BernoulliDistributionGenerator probGenerator(dependency->probability());
-            successful = probGenerator.random(randomGenerator);
+            successful = randomGenerator.randomProbability() <= dependency->probability();
         }
         STORM_LOG_TRACE("Let dependency " << *dependency << " " << (successful ? "successfully" : "unsuccessfully") << " fail");
         return std::make_tuple(iterFailable, 0, successful);
@@ -75,15 +75,13 @@ std::tuple<storm::dft::storage::FailableElements::const_iterator, double, bool>
         // Initialize with first BE
         storm::dft::storage::FailableElements::const_iterator nextFail = iterFailable;
         double rate = state->getBERate(nextFail.asBE(dft)->id());
-        storm::utility::ExponentialDistributionGenerator rateGenerator(rate);
-        double smallestTimebound = rateGenerator.random(randomGenerator);
+        double smallestTimebound = randomGenerator.randomExponential(rate);
         ++iterFailable;
 
         // Consider all other BEs and find the one which fails first
         for (; iterFailable != state->getFailableElements().end(); ++iterFailable) {
             rate = state->getBERate(iterFailable.asBE(dft)->id());
-            rateGenerator = storm::utility::ExponentialDistributionGenerator(rate);
-            double timebound = rateGenerator.random(randomGenerator);
+            double timebound = randomGenerator.randomExponential(rate);
             if (timebound < smallestTimebound) {
                 // BE fails earlier -> use as nextFail
                 nextFail = iterFailable;
diff --git a/src/storm-dft/simulator/DFTTraceSimulator.h b/src/storm-dft/simulator/DFTTraceSimulator.h
index 4fa0b92002..f3118871d8 100644
--- a/src/storm-dft/simulator/DFTTraceSimulator.h
+++ b/src/storm-dft/simulator/DFTTraceSimulator.h
@@ -4,7 +4,6 @@
 #include "storm-dft/storage/DFT.h"
 #include "storm-dft/storage/DFTState.h"
 #include "storm-dft/storage/FailableElements.h"
-
 #include "storm/utility/random.h"
 
 namespace storm::dft {
@@ -41,14 +40,14 @@ class DFTTraceSimulator {
      * @param randomGenerator Random number generator.
      */
     DFTTraceSimulator(storm::dft::storage::DFT<ValueType> const& dft, storm::dft::storage::DFTStateGenerationInfo const& stateGenerationInfo,
-                      boost::mt19937& randomGenerator);
+                      storm::utility::RandomProbabilityGenerator<ValueType> randomGenerator);
 
     /*!
      * Set the random number generator.
      *
-     * @param randomNumberGenerator Random number generator.
+     * @param randomGenerator Random number generator.
      */
-    void setRandomNumberGenerator(boost::mt19937& randomNumberGenerator);
+    void setRandomGenerator(storm::utility::RandomProbabilityGenerator<ValueType> randomGenerator);
 
     /*!
      * Set the current state back to the initial state in order to start a new simulation.
@@ -144,8 +143,8 @@ class DFTTraceSimulator {
     // Currently elapsed time
     double time;
 
-    // Random number generator
-    boost::mt19937& randomGenerator;
+    // Random generator
+    storm::utility::RandomProbabilityGenerator<ValueType> randomGenerator;
 };
 
 }  // namespace simulator
diff --git a/src/storm/simulator/DiscreteTimeSparseModelSimulator.cpp b/src/storm/simulator/DiscreteTimeSparseModelSimulator.cpp
deleted file mode 100644
index d45ce861bd..0000000000
--- a/src/storm/simulator/DiscreteTimeSparseModelSimulator.cpp
+++ /dev/null
@@ -1,98 +0,0 @@
-#include "storm/simulator/DiscreteTimeSparseModelSimulator.h"
-#include "storm/models/sparse/Model.h"
-
-namespace storm {
-namespace simulator {
-template<typename ValueType, typename RewardModelType>
-DiscreteTimeSparseModelSimulator<ValueType, RewardModelType>::DiscreteTimeSparseModelSimulator(
-    storm::models::sparse::Model<ValueType, RewardModelType> const& model)
-    : model(model), currentState(*model.getInitialStates().begin()), zeroRewards(model.getNumberOfRewardModels(), storm::utility::zero<ValueType>()) {
-    STORM_LOG_WARN_COND(model.getInitialStates().getNumberOfSetBits() == 1,
-                        "The model has multiple initial states. This simulator assumes it starts from the initial state with the lowest index.");
-    lastRewards = zeroRewards;
-    uint64_t i = 0;
-    for (auto const& rewModPair : model.getRewardModels()) {
-        if (rewModPair.second.hasStateRewards()) {
-            lastRewards[i] += rewModPair.second.getStateReward(currentState);
-        }
-        ++i;
-    }
-}
-
-template<typename ValueType, typename RewardModelType>
-void DiscreteTimeSparseModelSimulator<ValueType, RewardModelType>::setSeed(uint64_t seed) {
-    generator = storm::utility::RandomProbabilityGenerator<ValueType>(seed);
-}
-
-template<typename ValueType, typename RewardModelType>
-bool DiscreteTimeSparseModelSimulator<ValueType, RewardModelType>::randomStep() {
-    // TODO random_uint is slow
-    if (model.getTransitionMatrix().getRowGroupSize(currentState) == 0) {
-        return false;
-    }
-    return step(generator.random_uint(0, model.getTransitionMatrix().getRowGroupSize(currentState) - 1));
-}
-
-template<typename ValueType, typename RewardModelType>
-bool DiscreteTimeSparseModelSimulator<ValueType, RewardModelType>::step(uint64_t action) {
-    // TODO lots of optimization potential.
-    //  E.g., do not sample random numbers if there is only a single transition
-    lastRewards = zeroRewards;
-    ValueType probability = generator.random();
-    STORM_LOG_ASSERT(action < model.getTransitionMatrix().getRowGroupSize(currentState), "Action index higher than number of actions");
-    uint64_t row = model.getTransitionMatrix().getRowGroupIndices()[currentState] + action;
-    uint64_t i = 0;
-    for (auto const& rewModPair : model.getRewardModels()) {
-        if (rewModPair.second.hasStateActionRewards()) {
-            lastRewards[i] += rewModPair.second.getStateActionReward(row);
-        }
-        ++i;
-    }
-    ValueType sum = storm::utility::zero<ValueType>();
-    for (auto const& entry : model.getTransitionMatrix().getRow(row)) {
-        sum += entry.getValue();
-        if (sum >= probability) {
-            currentState = entry.getColumn();
-            i = 0;
-            for (auto const& rewModPair : model.getRewardModels()) {
-                if (rewModPair.second.hasStateRewards()) {
-                    lastRewards[i] += rewModPair.second.getStateReward(currentState);
-                }
-                ++i;
-            }
-            return true;
-        }
-    }
-    // This position should never be reached
-    return false;
-}
-
-template<typename ValueType, typename RewardModelType>
-uint64_t DiscreteTimeSparseModelSimulator<ValueType, RewardModelType>::getCurrentState() const {
-    return currentState;
-}
-
-template<typename ValueType, typename RewardModelType>
-bool DiscreteTimeSparseModelSimulator<ValueType, RewardModelType>::resetToInitial() {
-    currentState = *model.getInitialStates().begin();
-    lastRewards = zeroRewards;
-    uint64_t i = 0;
-    for (auto const& rewModPair : model.getRewardModels()) {
-        if (rewModPair.second.hasStateRewards()) {
-            lastRewards[i] += rewModPair.second.getStateReward(currentState);
-        }
-        ++i;
-    }
-    return true;
-}
-
-template<typename ValueType, typename RewardModelType>
-std::vector<ValueType> const& DiscreteTimeSparseModelSimulator<ValueType, RewardModelType>::getLastRewards() const {
-    return lastRewards;
-}
-
-template class DiscreteTimeSparseModelSimulator<double>;
-template class DiscreteTimeSparseModelSimulator<storm::RationalNumber>;
-
-}  // namespace simulator
-}  // namespace storm
diff --git a/src/storm/simulator/DiscreteTimeSparseModelSimulator.h b/src/storm/simulator/DiscreteTimeSparseModelSimulator.h
deleted file mode 100644
index 88f291e779..0000000000
--- a/src/storm/simulator/DiscreteTimeSparseModelSimulator.h
+++ /dev/null
@@ -1,35 +0,0 @@
-#include <cstdint>
-#include "storm/models/sparse/Model.h"
-#include "storm/utility/random.h"
-
-namespace storm {
-namespace simulator {
-
-/**
- * This class is a low-level interface to quickly sample from Discrete-Time Models
- * stored explicitly as a SparseModel.
- * Additional information about state, actions, should be obtained via the model itself.
- *
- * TODO: It may be nice to write a CPP wrapper that does not require to actually obtain such informations yourself.
- * @tparam ModelType
- */
-template<typename ValueType, typename RewardModelType = storm::models::sparse::StandardRewardModel<ValueType>>
-class DiscreteTimeSparseModelSimulator {
-   public:
-    DiscreteTimeSparseModelSimulator(storm::models::sparse::Model<ValueType, RewardModelType> const& model);
-    void setSeed(uint64_t);
-    bool step(uint64_t action);
-    bool randomStep();
-    std::vector<ValueType> const& getLastRewards() const;
-    uint64_t getCurrentState() const;
-    bool resetToInitial();
-
-   protected:
-    storm::models::sparse::Model<ValueType, RewardModelType> const& model;
-    uint64_t currentState;
-    std::vector<ValueType> lastRewards;
-    std::vector<ValueType> zeroRewards;
-    storm::utility::RandomProbabilityGenerator<ValueType> generator;
-};
-}  // namespace simulator
-}  // namespace storm
\ No newline at end of file
diff --git a/src/storm/simulator/ModelSimulator.cpp b/src/storm/simulator/ModelSimulator.cpp
new file mode 100644
index 0000000000..5c17a54949
--- /dev/null
+++ b/src/storm/simulator/ModelSimulator.cpp
@@ -0,0 +1,61 @@
+#include "storm/simulator/ModelSimulator.h"
+
+#include "storm/utility/constants.h"
+#include "storm/utility/macros.h"
+
+namespace storm {
+namespace simulator {
+
+template<typename ValueType>
+ModelSimulator<ValueType>::ModelSimulator() : currentTime(storm::utility::zero<ValueType>()), randomGenerator() {
+    // Intentionally left empty
+}
+
+template<typename ValueType>
+void ModelSimulator<ValueType>::setSeed(uint64_t newSeed) {
+    randomGenerator = storm::utility::RandomProbabilityGenerator<ValueType>(newSeed);
+}
+
+template<typename ValueType>
+bool ModelSimulator<ValueType>::randomStep() {
+    if (isContinuousTimeModel()) {
+        // First choose time when to leave the state
+        randomTime();
+    }
+
+    if (isCurrentStateDeadlock()) {
+        return false;
+    }
+    if (getCurrentNumberOfChoices() == 1) {
+        return step(0);
+    }
+    // Select action by uniform distribution
+    return step(this->randomGenerator.randomSelect(0, getCurrentNumberOfChoices() - 1));
+}
+
+template<typename ValueType>
+void ModelSimulator<ValueType>::randomTime() {
+    STORM_LOG_ASSERT(isContinuousTimeModel(), "Model must be continuous-time model");
+    ValueType exitRate = getCurrentExitRate();
+    if (!storm::utility::isZero(exitRate)) {
+        STORM_LOG_ASSERT(getCurrentNumberOfChoices() == 1, "Markovian state should have a trivial row group.");
+        ValueType time = this->randomGenerator.randomExponential(exitRate);
+        this->currentTime += time;
+    }
+}
+
+template<typename ValueType>
+ValueType ModelSimulator<ValueType>::getCurrentTime() const {
+    return currentTime;
+}
+
+template<typename ValueType>
+std::vector<ValueType> const& ModelSimulator<ValueType>::getCurrentRewards() const {
+    return currentRewards;
+}
+
+template class ModelSimulator<double>;
+template class ModelSimulator<storm::RationalNumber>;
+
+}  // namespace simulator
+}  // namespace storm
diff --git a/src/storm/simulator/ModelSimulator.h b/src/storm/simulator/ModelSimulator.h
new file mode 100644
index 0000000000..7be755097e
--- /dev/null
+++ b/src/storm/simulator/ModelSimulator.h
@@ -0,0 +1,128 @@
+#pragma once
+
+#include "storm/utility/random.h"
+
+namespace storm {
+namespace simulator {
+
+/*!
+ * Abstract class for simulator of model.
+ */
+template<typename ValueType>
+class ModelSimulator {
+   public:
+    /*!
+     * Constructor.
+     */
+    ModelSimulator();
+
+    /*!
+     * Destructor.
+     */
+    virtual ~ModelSimulator() = default;
+
+    /*!
+     * Set specific seed for random number generator.
+     * Useful to obtain deterministic behaviour of the simulator.
+     *
+     * @param seed Seed.
+     */
+    void setSeed(uint64_t seed);
+
+    /*!
+     * Reset the simulator to the initial state.
+     */
+    virtual void resetToInitial() = 0;
+
+    /*!
+     * Perform a random step from the current state.
+     * For probabilistic states: Randomly select an action and then randomly select a transition.
+     * For Markovian states (only in continuous-time models): Increase the time according to exit rate and then randomly select a transition.
+     *
+     * @return Whether the step was successful or not (e.g. deadlock).
+     */
+    bool randomStep();
+
+    /*
+     * Increase the current time according to the exit rate in the current state (if continuous-time model).
+     * Nothing happens if for discrete-time models.
+     */
+    void randomTime();
+
+    /*!
+     * Perform the given action and then randomly select a transition.
+     *
+     * @param action Action index.
+     * @return True iff step was successful.
+     */
+    virtual bool step(uint64_t action) = 0;
+
+    /*!
+     * Time progressed so far.
+     * Is 0 for discrete-time models.
+     *
+     * @return Progressed time.
+     */
+    ValueType getCurrentTime() const;
+
+    /*!
+     * Return the number of choices for the current state.
+     *
+     * @return Number of current choices.
+     */
+    virtual uint64_t getCurrentNumberOfChoices() const = 0;
+
+    /*!
+     * Return all labels for the current state.
+     * @return State labels.
+     */
+    virtual std::set<std::string> getCurrentStateLabelling() const = 0;
+
+    /*!
+     * Get the rewards for the current state.
+     *
+     * @return The current rewards (in the order from getRewardNames()).
+     */
+    std::vector<ValueType> const& getCurrentRewards() const;
+
+    /*!
+     * The names of the rewards that are returned.
+     * @return Names of reward models.
+     */
+    virtual std::vector<std::string> getRewardNames() const = 0;
+
+    /*!
+     * Get the current exit rate of the state.
+     * @return State exit rate (if Markovian state in a continuous-time model) or 0 (if probabilistic state or a discrete-time model).
+     */
+    virtual ValueType getCurrentExitRate() const = 0;
+
+    /*!
+     * Whether the current state is a deadlock.
+     *
+     * @return True if the current state has no outgoing transitions.
+     */
+    virtual bool isCurrentStateDeadlock() const = 0;
+
+    /*!
+     * Whether the model is a continuous-time model.
+     *
+     * @return True if the model is continuous-time.
+     */
+    virtual bool isContinuousTimeModel() const = 0;
+
+   protected:
+    /// Time which has progressed so far. Only relevant for continuous-time models
+    ValueType currentTime;
+
+    /// The current rewards
+    std::vector<ValueType> currentRewards;
+
+    /// Used to reinitialize the reward vector to all 0.
+    std::vector<ValueType> zeroRewards;
+
+    /// Random number generator
+    storm::utility::RandomProbabilityGenerator<ValueType> randomGenerator;
+};
+}  // namespace simulator
+}  // namespace storm
\ No newline at end of file
diff --git a/src/storm/simulator/PrismProgramSimulator.cpp b/src/storm/simulator/PrismProgramSimulator.cpp
index fa54f1763f..f7f230a445 100644
--- a/src/storm/simulator/PrismProgramSimulator.cpp
+++ b/src/storm/simulator/PrismProgramSimulator.cpp
@@ -1,6 +1,6 @@
 #include "storm/simulator/PrismProgramSimulator.h"
+
 #include "storm/adapters/JsonAdapter.h"
-#include "storm/exceptions/NotSupportedException.h"
 #include "storm/storage/expressions/ExpressionEvaluator.h"
 
 using namespace storm::generator;
@@ -9,29 +9,48 @@ namespace storm {
 namespace simulator {
 
 template<typename ValueType>
-DiscreteTimePrismProgramSimulator<ValueType>::DiscreteTimePrismProgramSimulator(storm::prism::Program const& program,
-                                                                                storm::generator::NextStateGeneratorOptions const& options)
-    : program(program),
+PrismProgramSimulator<ValueType>::PrismProgramSimulator(storm::prism::Program program, storm::generator::NextStateGeneratorOptions const& options)
+    : storm::simulator::ModelSimulator<ValueType>(),
+      program(program),
       currentState(),
       stateGenerator(std::make_shared<storm::generator::PrismNextStateGenerator<ValueType, uint32_t>>(program, options)),
-      zeroRewards(stateGenerator->getNumberOfRewardModels(), storm::utility::zero<ValueType>()),
-      lastActionRewards(zeroRewards),
       stateToId(stateGenerator->getStateSize()),
       idToState() {
+    this->zeroRewards = std::vector<ValueType>(stateGenerator->getNumberOfRewardModels(), storm::utility::zero<ValueType>());
     // Current state needs to be overwritten to actual initial state.
-    resetToInitial();
+    this->resetToInitial();
 }
 
 template<typename ValueType>
-void DiscreteTimePrismProgramSimulator<ValueType>::setSeed(uint64_t newSeed) {
-    generator = storm::utility::RandomProbabilityGenerator<ValueType>(newSeed);
-}
+bool PrismProgramSimulator<ValueType>::step(uint64_t actionNumber) {
+    STORM_LOG_ASSERT(actionNumber < getCurrentNumberOfChoices(), "Action index higher than number of actions");
+    auto choice = getChoices()[actionNumber];
+
+    // Set state-action reward
+    this->currentRewards = this->zeroRewards;
+    if (!isContinuousTimeModel()) {
+        this->currentRewards = choice.getRewards();
+        STORM_LOG_ASSERT(this->currentRewards.size() == stateGenerator->getNumberOfRewardModels(), "Reward vector should have as many rewards as model.");
+    } else {
+        // TODO add support
+        STORM_LOG_WARN("State-action rewards are currently not supported for continuous-time models.");
+    }
+
+    uint64_t nextState;
+    if (choice.size() == 1) {
+        // Select only transition
+        nextState = choice.begin()->first;
+    } else {
+        // Randomly select transition
+        ValueType probability = this->randomGenerator.randomProbability();
+        if (program.getModelType() == storm::prism::Program::ModelType::CTMC ||
+            (program.getModelType() == storm::prism::Program::ModelType::MA && choice.isMarkovian())) {
+            // Scale probability to exit rate
+            probability *= getCurrentExitRate();
+        }
+        nextState = choice.sampleFromDistribution(probability);
+    }
 
-template<typename ValueType>
-bool DiscreteTimePrismProgramSimulator<ValueType>::step(uint64_t actionNumber) {
-    uint32_t nextState = behavior.getChoices()[actionNumber].sampleFromDistribution(generator.random());
-    lastActionRewards = behavior.getChoices()[actionNumber].getRewards();
-    STORM_LOG_ASSERT(lastActionRewards.size() == stateGenerator->getNumberOfRewardModels(), "Reward vector should have as many rewards as model.");
     currentState = idToState[nextState];
     // TODO we do not need to do this in every step!
     clearStateCaches();
@@ -40,24 +59,29 @@ bool DiscreteTimePrismProgramSimulator<ValueType>::step(uint64_t actionNumber) {
 }
 
 template<typename ValueType>
-bool DiscreteTimePrismProgramSimulator<ValueType>::explore() {
+void PrismProgramSimulator<ValueType>::explore() {
     // Load the current state into the next state generator.
     stateGenerator->load(currentState);
     // TODO: This low-level code currently expands all actions, while this is not necessary.
-    // However, using the next state generator ensures compatibliity with the model generator.
+    // However, using the next state generator ensures compatibility with the model generator.
     behavior = stateGenerator->expand(stateToIdCallback);
-    STORM_LOG_ASSERT(behavior.getStateRewards().size() == lastActionRewards.size(), "Reward vectors should have same length.");
+    STORM_LOG_ASSERT(behavior.getStateRewards().size() == this->currentRewards.size(), "Reward vectors should have same length.");
     for (uint64_t i = 0; i < behavior.getStateRewards().size(); i++) {
-        lastActionRewards[i] += behavior.getStateRewards()[i];
+        this->currentRewards[i] += behavior.getStateRewards()[i];
     }
-    return true;
 }
 
 template<typename ValueType>
-bool DiscreteTimePrismProgramSimulator<ValueType>::isSinkState() const {
-    if (behavior.empty()) {
+bool PrismProgramSimulator<ValueType>::isCurrentStateDeadlock() const {
+    return behavior.empty();
+}
+
+template<typename ValueType>
+bool PrismProgramSimulator<ValueType>::isSinkState() const {
+    if (isCurrentStateDeadlock()) {
         return true;
     }
+
     std::set<uint32_t> successorIds;
     for (Choice<ValueType, uint32_t> const& choice : behavior.getChoices()) {
         for (auto it = choice.begin(); it != choice.end(); ++it) {
@@ -74,76 +98,106 @@ bool DiscreteTimePrismProgramSimulator<ValueType>::isSinkState() const {
 }
 
 template<typename ValueType>
-std::vector<std::string> DiscreteTimePrismProgramSimulator<ValueType>::getCurrentStateLabelling() const {
-    std::vector<std::string> labels;
+std::set<std::string> PrismProgramSimulator<ValueType>::getCurrentStateLabelling() const {
+    std::set<std::string> labels;
     for (auto const& label : program.getLabels()) {
         if (stateGenerator->evaluateBooleanExpressionInCurrentState(label.getStatePredicateExpression())) {
-            labels.push_back(label.getName());
+            labels.insert(label.getName());
         }
     }
     return labels;
 }
 
 template<typename ValueType>
-std::vector<generator::Choice<ValueType, uint32_t>> const& DiscreteTimePrismProgramSimulator<ValueType>::getChoices() const {
+std::vector<generator::Choice<ValueType, uint32_t>> const& PrismProgramSimulator<ValueType>::getChoices() const {
     return behavior.getChoices();
 }
 
 template<typename ValueType>
-std::vector<ValueType> const& DiscreteTimePrismProgramSimulator<ValueType>::getLastRewards() const {
-    return lastActionRewards;
+uint64_t PrismProgramSimulator<ValueType>::getCurrentNumberOfChoices() const {
+    return getChoices().size();
 }
 
 template<typename ValueType>
-CompressedState const& DiscreteTimePrismProgramSimulator<ValueType>::getCurrentState() const {
+ValueType PrismProgramSimulator<ValueType>::getCurrentExitRate() const {
+    if (!isContinuousTimeModel()) {
+        // Discrete-time models have no exit rate
+        return storm::utility::zero<ValueType>();
+    }
+
+    if (isCurrentStateDeadlock()) {
+        return storm::utility::zero<ValueType>();
+    }
+    if (getCurrentNumberOfChoices() > 1) {
+        // Probabilistic state
+        return storm::utility::zero<ValueType>();
+    }
+
+    STORM_LOG_ASSERT(getCurrentNumberOfChoices() == 1, "Expected deterministic state.");
+    auto choice = getChoices().front();
+    ValueType totalMass = choice.getTotalMass();
+    // CTMC or (MA and Markovian) => totalMass > 0
+    STORM_LOG_ASSERT(!(program.getModelType() == storm::prism::Program::ModelType::CTMC ||
+                       (program.getModelType() == storm::prism::Program::ModelType::MA && choice.isMarkovian())) ||
+                         storm::utility::isPositive(totalMass),
+                     "Expected positive exit rate.");
+    // MA and not Markovian => totalMass = 1
+    STORM_LOG_ASSERT(!(program.getModelType() == storm::prism::Program::ModelType::MA && !choice.isMarkovian()) || storm::utility::isOne(totalMass),
+                     "Probabilities should sum up to one.");
+    return totalMass;
+}
+
+template<typename ValueType>
+CompressedState const& PrismProgramSimulator<ValueType>::getCurrentState() const {
     return currentState;
 }
 
 template<typename ValueType>
-expressions::SimpleValuation DiscreteTimePrismProgramSimulator<ValueType>::getCurrentStateAsValuation() const {
+expressions::SimpleValuation PrismProgramSimulator<ValueType>::getCurrentStateAsValuation() const {
     return unpackStateIntoValuation(currentState, stateGenerator->getVariableInformation(), program.getManager());
 }
 
 template<typename ValueType>
-std::string DiscreteTimePrismProgramSimulator<ValueType>::getCurrentStateString() const {
+std::string PrismProgramSimulator<ValueType>::getCurrentStateString() const {
     return stateGenerator->stateToString(currentState);
 }
 
 template<typename ValueType>
-storm::json<ValueType> DiscreteTimePrismProgramSimulator<ValueType>::getStateAsJson() const {
+storm::json<ValueType> PrismProgramSimulator<ValueType>::getStateAsJson() const {
     return stateGenerator->currentStateToJson(false);
 }
 
 template<typename ValueType>
-storm::json<ValueType> DiscreteTimePrismProgramSimulator<ValueType>::getObservationAsJson() const {
+storm::json<ValueType> PrismProgramSimulator<ValueType>::getObservationAsJson() const {
     return stateGenerator->currentStateToJson(true);
 }
 
 template<typename ValueType>
-bool DiscreteTimePrismProgramSimulator<ValueType>::resetToInitial() {
-    lastActionRewards = zeroRewards;
+void PrismProgramSimulator<ValueType>::resetToInitial() {
     auto indices = stateGenerator->getInitialStates(stateToIdCallback);
-    STORM_LOG_THROW(indices.size() == 1, storm::exceptions::NotSupportedException, "Program must have a unique initial state");
-    currentState = idToState[indices[0]];
-    return explore();
+    STORM_LOG_WARN_COND(indices.size() == 1,
+                        "The model has multiple initial states. This simulator assumes it starts from the initial state with the lowest index.");
+    resetToState(idToState[indices[0]], storm::utility::zero<ValueType>());
 }
 
 template<typename ValueType>
-bool DiscreteTimePrismProgramSimulator<ValueType>::resetToState(generator::CompressedState const& newState) {
-    lastActionRewards = zeroRewards;
+void PrismProgramSimulator<ValueType>::resetToState(generator::CompressedState const& newState, ValueType time) {
+    this->currentRewards = this->zeroRewards;
+    this->currentTime = time;
     currentState = newState;
-    return explore();
+    explore();
 }
 
 template<typename ValueType>
-bool DiscreteTimePrismProgramSimulator<ValueType>::resetToState(expressions::SimpleValuation const& valuation) {
-    lastActionRewards = zeroRewards;
+void PrismProgramSimulator<ValueType>::resetToState(expressions::SimpleValuation const& valuation, ValueType time) {
+    this->currentRewards = this->zeroRewards;
+    this->currentTime = time;
     currentState = generator::packStateFromValuation(valuation, stateGenerator->getVariableInformation(), true);
-    return explore();
+    explore();
 }
 
 template<typename ValueType>
-std::vector<std::string> DiscreteTimePrismProgramSimulator<ValueType>::getRewardNames() const {
+std::vector<std::string> PrismProgramSimulator<ValueType>::getRewardNames() const {
     std::vector<std::string> names;
     for (uint64_t i = 0; i < stateGenerator->getNumberOfRewardModels(); ++i) {
         names.push_back(stateGenerator->getRewardModelInformation(i).getName());
@@ -152,7 +206,12 @@ std::vector<std::string> DiscreteTimePrismProgramSimulator<ValueType>::getReward
 }
 
 template<typename ValueType>
-uint32_t DiscreteTimePrismProgramSimulator<ValueType>::getOrAddStateIndex(generator::CompressedState const& state) {
+bool PrismProgramSimulator<ValueType>::isContinuousTimeModel() const {
+    return !program.isDiscreteTimeModel();
+}
+
+template<typename ValueType>
+uint32_t PrismProgramSimulator<ValueType>::getOrAddStateIndex(generator::CompressedState const& state) {
     uint32_t newIndex = static_cast<uint32_t>(stateToId.size());
 
     // Check, if the state was already registered.
@@ -166,11 +225,12 @@ uint32_t DiscreteTimePrismProgramSimulator<ValueType>::getOrAddStateIndex(genera
 }
 
 template<typename ValueType>
-void DiscreteTimePrismProgramSimulator<ValueType>::clearStateCaches() {
+void PrismProgramSimulator<ValueType>::clearStateCaches() {
     idToState.clear();
     stateToId = storm::storage::BitVectorHashMap<uint32_t>(stateGenerator->getStateSize());
 }
 
-template class DiscreteTimePrismProgramSimulator<double>;
+template class PrismProgramSimulator<double>;
+template class PrismProgramSimulator<storm::RationalNumber>;
 }  // namespace simulator
 }  // namespace storm
diff --git a/src/storm/simulator/PrismProgramSimulator.h b/src/storm/simulator/PrismProgramSimulator.h
index ef8473a923..d7fad35dc3 100644
--- a/src/storm/simulator/PrismProgramSimulator.h
+++ b/src/storm/simulator/PrismProgramSimulator.h
@@ -1,14 +1,15 @@
 #pragma once
 
+#include "storm/simulator/ModelSimulator.h"
+
 #include "storm/generator/PrismNextStateGenerator.h"
 #include "storm/storage/expressions/SimpleValuation.h"
 #include "storm/storage/prism/Program.h"
-#include "storm/utility/random.h"
 
 namespace storm {
 namespace simulator {
 
-/**
+/*!
  * This class provides a simulator interface on the prism program,
  * and uses the next state generator. While the next state generator has been tuned,
  * it is not targeted for simulation purposes. In particular, we (as of now)
@@ -18,101 +19,126 @@ namespace simulator {
  * On the other hand, this simulator is convenient for stepping through the model
  * as it potentially allows considering the next states.
  * Thus, while a performant alternative would be great, this simulator has its own merits.
- *
- * @tparam ValueType
  */
 template<typename ValueType>
-class DiscreteTimePrismProgramSimulator {
+class PrismProgramSimulator : public ModelSimulator<ValueType> {
    public:
-    /**
+    /*!
      * Initialize the simulator for a given prism program.
      *
      * @param program The prism program. Should have a unique initial state.
      * @param options The generator options that are used to generate successor states.
      */
-    DiscreteTimePrismProgramSimulator(storm::prism::Program const& program, storm::generator::NextStateGeneratorOptions const& options);
-    /**
-     * Set the simulation seed.
+    PrismProgramSimulator(storm::prism::Program program, storm::generator::NextStateGeneratorOptions const& options);
+
+    void resetToInitial() override;
+
+    /*!
+     * Reset the current state to the given state.
+     * @param compressedState State to reset to.
+     * @param time Time to reset to.
      */
-    void setSeed(uint64_t);
-    /**
-     *
+    void resetToState(generator::CompressedState const& compressedState, ValueType time);
+
+    /*!
+     * Reset the current state to the given state.
+     * @param valuationState State to reset to.
+     * @param time Time to reset to.
+     */
+    void resetToState(expressions::SimpleValuation const& valuationState, ValueType time);
+
+    bool step(uint64_t actionNumber) override;
+
+    /*!
+     * Get choices in current state.
      * @return A list of choices that encode the possibilities in the current state.
-     * @note successor states are encoded using state indices that will potentially be invalidated as soon as the internal state of the simulator changes
+     * @note Successor states are encoded using state indices that will potentially be invalidated as soon as the internal state of the simulator changes.
      */
     std::vector<generator::Choice<ValueType, uint32_t>> const& getChoices() const;
 
-    bool isSinkState() const;
+    uint64_t getCurrentNumberOfChoices() const override;
 
-    /**
-     * Make a step and randomly select the successor. The action is given as an argument, the index reflects the index of the getChoices vector that can be
-     * accessed.
-     *
-     * @param actionNumber The action to select.
-     * @return true, if this action can be taken.
+    /*!
+     * Get current state.
+     * @return Current state.
      */
-    bool step(uint64_t actionNumber);
-    /**
-     * Accessor for the last state action reward and the current state reward, added together.
-     * @return A vector with te number of rewards.
-     */
-    std::vector<ValueType> const& getLastRewards() const;
     generator::CompressedState const& getCurrentState() const;
+
+    /*!
+     * Get valuation of current state.
+     * @return Valuation.
+     */
     expressions::SimpleValuation getCurrentStateAsValuation() const;
-    std::vector<std::string> getCurrentStateLabelling() const;
 
+    /*!
+     * Get string representation of current state.
+     * @return String representation of current state.
+     */
+    std::string getCurrentStateString() const;
+
+    /*!
+     * Get json representation of current state.
+     * @return Json representation of current state.
+     */
     storm::json<ValueType> getStateAsJson() const;
 
-    storm::json<ValueType> getObservationAsJson() const;
+    std::set<std::string> getCurrentStateLabelling() const override;
 
-    std::string getCurrentStateString() const;
-    /**
-     * Reset to the (unique) initial state.
-     *
-     * @return
+    /*!
+     * Get json representation of observations.
+     * @return Json representation of observations.
      */
-    bool resetToInitial();
+    storm::json<ValueType> getObservationAsJson() const;
+
+    std::vector<std::string> getRewardNames() const override;
+
+    ValueType getCurrentExitRate() const override;
+
+    bool isCurrentStateDeadlock() const override;
 
-    bool resetToState(generator::CompressedState const& compressedState);
+    /*!
+     * Whether the current state is an absorbing state.
+     * @return True if current state is absorbing.
+     */
+    bool isSinkState() const;
 
-    bool resetToState(expressions::SimpleValuation const& valuationState);
+    bool isContinuousTimeModel() const override;
 
-    /**
-     * The names of the rewards that are returned.
+   private:
+    /*!
+     * Explore the current state and create choices and distribution over successor states.
      */
-    std::vector<std::string> getRewardNames() const;
+    void explore();
 
-   protected:
-    bool explore();
+    /*!
+     * Clear all state caches.
+     */
     void clearStateCaches();
-    /**
+
+    /*!
      * Helper function for (temp) storing states.
+     *
+     * @param state State.
+     * @return Index of state.
      */
-    uint32_t getOrAddStateIndex(generator::CompressedState const&);
+    uint32_t getOrAddStateIndex(generator::CompressedState const& state);
 
     /// The program that we are simulating.
-    storm::prism::Program const& program;
+    storm::prism::Program program;
     /// The current state in the program, in its compressed form.
     generator::CompressedState currentState;
     /// Generator for the next states
     std::shared_ptr<storm::generator::PrismNextStateGenerator<ValueType, uint32_t>> stateGenerator;
     /// Obtained behavior of a state
     generator::StateBehavior<ValueType> behavior;
-    /// Helper for last action reward construction
-    std::vector<ValueType> zeroRewards;
-    /// Stores the action rewards from the last action.
-    std::vector<ValueType> lastActionRewards;
-    /// Random number generator
-    storm::utility::RandomProbabilityGenerator<ValueType> generator;
-    /// Data structure to temp store states.
+    /// Data structure to temporarily store states.
     storm::storage::BitVectorHashMap<uint32_t> stateToId;
-
+    /// Data structure to temporarily store ids to states.
     std::unordered_map<uint32_t, generator::CompressedState> idToState;
 
-   private:
-    // Create a callback for the next-state generator to enable it to request the index of states.
+    /// Create a callback for the next-state generator to enable it to request the index of states.
     std::function<uint32_t(generator::CompressedState const&)> stateToIdCallback =
-        std::bind(&DiscreteTimePrismProgramSimulator<ValueType>::getOrAddStateIndex, this, std::placeholders::_1);
+        std::bind(&PrismProgramSimulator<ValueType>::getOrAddStateIndex, this, std::placeholders::_1);
 };
 }  // namespace simulator
 }  // namespace storm
diff --git a/src/storm/simulator/SparseModelSimulator.cpp b/src/storm/simulator/SparseModelSimulator.cpp
new file mode 100644
index 0000000000..a0783b069a
--- /dev/null
+++ b/src/storm/simulator/SparseModelSimulator.cpp
@@ -0,0 +1,170 @@
+#include "storm/simulator/SparseModelSimulator.h"
+
+#include "storm/models/sparse/Ctmc.h"
+#include "storm/models/sparse/MarkovAutomaton.h"
+
+namespace storm {
+namespace simulator {
+template<typename ValueType, typename RewardModelType>
+SparseModelSimulator<ValueType, RewardModelType>::SparseModelSimulator(std::shared_ptr<storm::models::sparse::Model<ValueType, RewardModelType> const> model)
+    : storm::simulator::ModelSimulator<ValueType>(), model(model) {
+    this->zeroRewards = std::vector<ValueType>(model->getNumberOfRewardModels(), storm::utility::zero<ValueType>());
+    this->resetToInitial();
+
+    if (this->isContinuousTimeModel()) {
+        if (model->getType() == storm::models::ModelType::Ctmc) {
+            exitRates = model->template as<storm::models::sparse::Ctmc<ValueType, RewardModelType>>()->getExitRateVector();
+        } else {
+            STORM_LOG_ASSERT(model->getType() == storm::models::ModelType::MarkovAutomaton, "Model is not a Markov automaton.");
+            auto ma = model->template as<storm::models::sparse::MarkovAutomaton<ValueType, RewardModelType>>();
+            exitRates = ma->getExitRates();
+#ifndef NDEBUG
+            for (size_t i = 0; i < exitRates.size(); ++i) {
+                STORM_LOG_ASSERT((storm::utility::isZero<ValueType>(exitRates[i]) || ma->isMarkovianState(i)) &&
+                                     (!storm::utility::isZero<ValueType>(exitRates[i]) || !ma->isMarkovianState(i)),
+                                 "Exit rate and Markovian state do not match.");
+            }
+#endif
+        }
+    }
+}
+
+template<typename ValueType, typename RewardModelType>
+void SparseModelSimulator<ValueType, RewardModelType>::resetToInitial() {
+    STORM_LOG_WARN_COND(model->getInitialStates().getNumberOfSetBits() == 1,
+                        "The model has multiple initial states. This simulator assumes it starts from the initial state with the lowest index.");
+    currentState = *model->getInitialStates().begin();
+    this->currentTime = storm::utility::zero<ValueType>();
+    // Set state rewards
+    this->currentRewards = this->zeroRewards;
+    uint64_t i = 0;
+    for (auto const& rewModPair : model->getRewardModels()) {
+        if (rewModPair.second.hasStateRewards()) {
+            this->currentRewards[i] += rewModPair.second.getStateReward(currentState);
+        }
+        ++i;
+    }
+}
+
+template<typename ValueType, typename RewardModelType>
+uint64_t SparseModelSimulator<ValueType, RewardModelType>::choice(uint64_t choice) {
+    STORM_LOG_ASSERT(choice < getCurrentNumberOfChoices(), "Action index higher than number of actions");
+    uint64_t row = model->getTransitionMatrix().getRowGroupIndices()[currentState] + choice;
+
+    // Set state-action reward
+    this->currentRewards = this->zeroRewards;
+    uint64_t i = 0;
+    for (auto const& rewModPair : model->getRewardModels()) {
+        if (rewModPair.second.hasStateActionRewards()) {
+            this->currentRewards[i] += rewModPair.second.getStateActionReward(row);
+        }
+        ++i;
+    }
+    return row;
+}
+
+template<typename ValueType, typename RewardModelType>
+void SparseModelSimulator<ValueType, RewardModelType>::transition(uint64_t row, uint64_t column) {
+#ifndef NDEBUG
+    bool columnExists = false;
+    for (auto const& entry : model->getTransitionMatrix().getRow(row)) {
+        if (entry.getColumn() == column) {
+            columnExists = true;
+        }
+    }
+    STORM_LOG_ASSERT(columnExists, "Column " << column << " does not exist");
+#endif
+
+    currentState = column;
+    // Add state reward
+    uint64_t i = 0;
+    for (auto const& rewModPair : model->getRewardModels()) {
+        if (rewModPair.second.hasStateRewards()) {
+            this->currentRewards[i] += rewModPair.second.getStateReward(currentState);
+        }
+        ++i;
+    }
+}
+
+template<typename ValueType, typename RewardModelType>
+bool SparseModelSimulator<ValueType, RewardModelType>::step(uint64_t action) {
+    // Perform action
+    uint64_t row = choice(action);
+
+    if (model->getTransitionMatrix().getRow(row).getNumberOfEntries() == 1) {
+        // Select the only transition
+        uint64_t column = model->getTransitionMatrix().getRow(row).begin()->getColumn();
+        transition(row, column);
+        return true;
+    }
+
+    // Randomly select transition
+    ValueType probability = this->randomGenerator.randomProbability();
+    if (model->getType() == storm::models::ModelType::Ctmc) {
+        // Scale probability to exit rate
+        probability *= getCurrentExitRate();
+    }
+    ValueType sum = storm::utility::zero<ValueType>();
+    for (auto const& entry : model->getTransitionMatrix().getRow(row)) {
+        sum += entry.getValue();
+        if (sum >= probability) {
+            transition(row, entry.getColumn());
+            return true;
+        }
+    }
+
+    STORM_LOG_ASSERT(false, "This line should not be reached.");
+    return false;
+}
+
+template<typename ValueType, typename RewardModelType>
+bool SparseModelSimulator<ValueType, RewardModelType>::step(uint64_t action, uint64_t column) {
+    uint64_t row = choice(action);
+    transition(row, column);
+    return true;
+}
+
+template<typename ValueType, typename RewardModelType>
+std::vector<std::string> SparseModelSimulator<ValueType, RewardModelType>::getRewardNames() const {
+    std::vector<std::string> names;
+    for (auto name : model->getRewardModels()) {
+        names.push_back(name.first);
+    }
+    return names;
+}
+
+template<typename ValueType, typename RewardModelType>
+ValueType SparseModelSimulator<ValueType, RewardModelType>::getCurrentExitRate() const {
+    return exitRates[currentState];
+}
+
+template<typename ValueType, typename RewardModelType>
+uint64_t SparseModelSimulator<ValueType, RewardModelType>::getCurrentState() const {
+    return currentState;
+}
+
+template<typename ValueType, typename RewardModelType>
+uint64_t SparseModelSimulator<ValueType, RewardModelType>::getCurrentNumberOfChoices() const {
+    return model->getTransitionMatrix().getRowGroupSize(currentState);
+}
+
+template<typename ValueType, typename RewardModelType>
+std::set<std::string> SparseModelSimulator<ValueType, RewardModelType>::getCurrentStateLabelling() const {
+    return model->getStateLabeling().getLabelsOfState(currentState);
+}
+
+template<typename ValueType, typename RewardModelType>
+bool SparseModelSimulator<ValueType, RewardModelType>::isCurrentStateDeadlock() const {
+    return model->getTransitionMatrix().getRowGroupSize(currentState) == 0;
+}
+
+template<typename ValueType, typename RewardModelType>
+bool SparseModelSimulator<ValueType, RewardModelType>::isContinuousTimeModel() const {
+    return !model->isDiscreteTimeModel();
+}
+
+template class SparseModelSimulator<double>;
+template class SparseModelSimulator<storm::RationalNumber>;
+
+}  // namespace simulator
+}  // namespace storm
diff --git a/src/storm/simulator/SparseModelSimulator.h b/src/storm/simulator/SparseModelSimulator.h
new file mode 100644
index 0000000000..71bbb685c6
--- /dev/null
+++ b/src/storm/simulator/SparseModelSimulator.h
@@ -0,0 +1,77 @@
+#pragma once
+
+#include "storm/models/sparse/Model.h"
+#include "storm/simulator/ModelSimulator.h"
+
+namespace storm {
+namespace simulator {
+
+/*!
+ * This class is a low-level interface to quickly sample from sparse models.
+ * Additional information about state, actions, should be obtained via the model itself.
+ *
+ * TODO: It may be nice to write a CPP wrapper that does not require to actually obtain such information yourself.
+ */
+template<typename ValueType, typename RewardModelType = storm::models::sparse::StandardRewardModel<ValueType>>
+class SparseModelSimulator : public storm::simulator::ModelSimulator<ValueType> {
+   public:
+    /*!
+     * Constructor.
+     * @param model Sparse model.
+     */
+    SparseModelSimulator(std::shared_ptr<storm::models::sparse::Model<ValueType, RewardModelType> const> model);
+
+    void resetToInitial() override;
+
+    bool step(uint64_t action) override;
+
+    /*!
+     * Perform action and then select transition corresponding to column.
+     */
+    bool step(uint64_t action, uint64_t column);
+
+    /*!
+     * Get current state id.
+     * @return Current state id.
+     */
+    uint64_t getCurrentState() const;
+
+    uint64_t getCurrentNumberOfChoices() const override;
+
+    std::set<std::string> getCurrentStateLabelling() const override;
+
+    std::vector<std::string> getRewardNames() const override;
+
+    ValueType getCurrentExitRate() const override;
+
+    bool isCurrentStateDeadlock() const override;
+
+    bool isContinuousTimeModel() const override;
+
+   private:
+    /*!
+     * Select choice, update state-action rewards and return corresponding row.
+     *
+     * @param choice Choice index.
+     * @return Row in sparse matrix corresponding to chosen state-action.
+     */
+    uint64_t choice(uint64_t choice);
+
+    /*!
+     * Perform transition from selected row to given column.
+     * @param row Row in sparse matrix corresponding to chosen state-action.
+     * @param column Successor state.
+     */
+    void transition(uint64_t row, uint64_t column);
+
+    /// The underlying sparse model
+    std::shared_ptr<storm::models::sparse::Model<ValueType, RewardModelType> const> model;
+
+    /// The current state id
+    uint64_t currentState;
+
+    /// Exit rates for continuous-time models
+    std::vector<ValueType> exitRates;
+};
+}  // namespace simulator
+}  // namespace storm
diff --git a/src/storm/simulator/TraceSimulator.cpp b/src/storm/simulator/TraceSimulator.cpp
new file mode 100644
index 0000000000..39a197df1a
--- /dev/null
+++ b/src/storm/simulator/TraceSimulator.cpp
@@ -0,0 +1,100 @@
+#include "storm/simulator/TraceSimulator.h"
+
+#include "storm/exceptions/UnsupportedModelException.h"
+#include "storm/utility/constants.h"
+#include "storm/utility/macros.h"
+
+namespace storm {
+namespace simulator {
+
+template<typename ValueType>
+TraceSimulator<ValueType>::TraceSimulator(std::shared_ptr<storm::simulator::ModelSimulator<ValueType>> simulator) : simulator(simulator) {
+    // Intentionally left empty.
+}
+
+template<typename ValueType>
+SimulationTraceResult TraceSimulator<ValueType>::simulateStepBoundedTrace(size_t stepBound, std::optional<std::string> const& goalLabel) {
+    simulator->resetToInitial();
+    size_t steps = 0;
+    while (steps <= stepBound) {
+        if (reachedLabel(goalLabel)) {
+            return SimulationTraceResult::GOAL_REACHED;
+        }
+        if (!simulator->randomStep()) {
+            return SimulationTraceResult::DEADLOCK;
+        }
+        steps++;
+    }
+    return SimulationTraceResult::BOUND_REACHED;
+}
+
+template<typename ValueType>
+SimulationTraceResult TraceSimulator<ValueType>::simulateTimeBoundedTrace(ValueType timeBound, std::optional<std::string> const& goalLabel) {
+    simulator->resetToInitial();
+
+    while (simulator->getCurrentTime() <= timeBound) {
+        if (reachedLabel(goalLabel)) {
+            return SimulationTraceResult::GOAL_REACHED;
+        }
+        if (!simulator->randomStep()) {
+            return SimulationTraceResult::DEADLOCK;
+        }
+    }
+    return SimulationTraceResult::BOUND_REACHED;
+}
+
+template<typename ValueType>
+ValueType TraceSimulator<ValueType>::simulateStepBoundedReachability(std::string const& label, uint64_t stepBound, uint64_t numberRuns) {
+    STORM_LOG_WARN_COND(!simulator->isContinuousTimeModel(), "Simulating step-bounded reachability on a continuous-time model.");
+    uint64_t numberGoal = 0;
+    storm::simulator::SimulationTraceResult res;
+    for (size_t i = 0; i < numberRuns; ++i) {
+        res = simulateStepBoundedTrace(stepBound, label);
+        if (res == storm::simulator::SimulationTraceResult::GOAL_REACHED) {
+            ++numberGoal;
+        } else {
+            if (res == storm::simulator::SimulationTraceResult::DEADLOCK) {
+                STORM_LOG_WARN("Simulation reached a deadlock state.");
+            }
+        }
+    }
+    return storm::utility::convertNumber<ValueType>(numberGoal) / storm::utility::convertNumber<ValueType>(numberRuns);
+}
+
+template<typename ValueType>
+ValueType TraceSimulator<ValueType>::simulateTimeBoundedReachability(std::string const& label, ValueType timebound, uint64_t numberRuns) {
+    STORM_LOG_THROW(simulator->isContinuousTimeModel(), storm::exceptions::UnsupportedModelException,
+                    "Time-bounded reachability can only be simulated on continuous-time models.");
+    uint64_t numberGoal = 0;
+    storm::simulator::SimulationTraceResult res;
+    for (size_t i = 0; i < numberRuns; ++i) {
+        res = simulateTimeBoundedTrace(timebound, label);
+        if (res == storm::simulator::SimulationTraceResult::GOAL_REACHED) {
+            ++numberGoal;
+        } else {
+            if (res == storm::simulator::SimulationTraceResult::DEADLOCK) {
+                STORM_LOG_WARN("Simulation reached a deadlock state.");
+            }
+        }
+    }
+    return storm::utility::convertNumber<ValueType>(numberGoal) / storm::utility::convertNumber<ValueType>(numberRuns);
+}
+
+template<typename ValueType>
+bool TraceSimulator<ValueType>::reachedLabel(std::optional<std::string> const& goalLabel) const {
+    if (!goalLabel.has_value()) {
+        return false;
+    }
+    return simulator->getCurrentStateLabelling().contains(goalLabel.value());
+}
+
+template<typename ValueType>
+std::shared_ptr<storm::simulator::ModelSimulator<ValueType> const> TraceSimulator<ValueType>::getSimulator() const {
+    return simulator;
+}
+
+template class TraceSimulator<double>;
+template class TraceSimulator<storm::RationalNumber>;
+
+}  // namespace simulator
+}  // namespace storm
diff --git a/src/storm/simulator/TraceSimulator.h b/src/storm/simulator/TraceSimulator.h
new file mode 100644
index 0000000000..60cb22d2f5
--- /dev/null
+++ b/src/storm/simulator/TraceSimulator.h
@@ -0,0 +1,82 @@
+#pragma once
+
+#include "storm/simulator/ModelSimulator.h"
+
+namespace storm {
+namespace simulator {
+
+/*!
+ * The result of simulating a trace can be either that
+ * - the goal label was reached
+ * - the step/time bound was reached
+ * - or a deadlock was encountered.
+ */
+enum class SimulationTraceResult { GOAL_REACHED, BOUND_REACHED, DEADLOCK };
+
+/*!
+ * Perform simulation runs of a model.
+ */
+template<typename ValueType>
+class TraceSimulator {
+   public:
+    /*!
+     * Constructor.
+     * @param simulator Underlying simulator.
+     */
+    TraceSimulator(std::shared_ptr<storm::simulator::ModelSimulator<ValueType>> simulator);
+
+    /*!
+     * Simulate a single (untimed) trace until either the goal label or the step bound is reached.
+     * @param stepBound Step bound.
+     * @param goalLabel Label to reach.
+     * @return Result of trace (step bound reached, goal reached or deadlock).
+     */
+    SimulationTraceResult simulateStepBoundedTrace(size_t stepBound, std::optional<std::string> const& goalLabel);
+
+    /*!
+     * Simulate a single (timed) trace until either the goal label or the time bound is reached.
+     * @param timeBound Time bound.
+     * @param goalLabel Label to reach.
+     * @return Result of trace (time bound reached, goal reached or deadlock).
+     */
+    SimulationTraceResult simulateTimeBoundedTrace(ValueType timeBound, std::optional<std::string> const& goalLabel);
+
+    /*!
+     * Perform Monte Carlo simulation of a step-bounded reachability property.
+     * Approximates P=? [F<=stepBound label].
+     * @param label Goal label to reach.
+     * @param stepBound Step bound.
+     * @param numberRuns Number of simulation runs to perform.
+     * @return Probability of satisfying the property. Calculated as the ratio of runs reaching the label within stepBound / numberRuns.
+     */
+    ValueType simulateStepBoundedReachability(std::string const& label, uint64_t stepBound, uint64_t numberRuns);
+
+    /*!
+     * Perform Monte Carlo simulation of a time-bounded reachability property (for continuous-time models).
+     * Approximates P=? [F<=timeBound label].
+     * @param label Goal label to reach.
+     * @param timebound Time bound.
+     * @param numberRuns Number of simulation runs to perform.
+     * @return Probability of satisfying the property. Calculated as the ratio of runs reaching the label within timeBound / numberRuns.
+     */
+    ValueType simulateTimeBoundedReachability(std::string const& label, ValueType timebound, uint64_t numberRuns);
+
+    /*!
+     * Return the underlying simulator.
+     * @return Simulator.
+     */
+    std::shared_ptr<storm::simulator::ModelSimulator<ValueType> const> getSimulator() const;
+
+   private:
+    /*!
+     * Check whether the current state of the simulator contains the given label.
+     * @param goalLabel Label to reach.
+     * @return True if current state contains goal label. Always returns false if no goal label was provided.
+     */
+    bool reachedLabel(std::optional<std::string> const& goalLabel) const;
+
+    /// The underlying simulator
+    std::shared_ptr<storm::simulator::ModelSimulator<ValueType>> simulator;
+};
+}  // namespace simulator
+}  // namespace storm
\ No newline at end of file
diff --git a/src/storm/utility/random.cpp b/src/storm/utility/random.cpp
index fbd0f9e2ff..44d3e93673 100644
--- a/src/storm/utility/random.cpp
+++ b/src/storm/utility/random.cpp
@@ -2,48 +2,59 @@
 
 #include <limits>
 
+#include "storm/exceptions/NotSupportedException.h"
+#include "storm/utility/macros.h"
+
 namespace storm {
 namespace utility {
-RandomProbabilityGenerator<double>::RandomProbabilityGenerator() : distribution(0.0, 1.0) {
+
+RandomProbabilityGenerator<double>::RandomProbabilityGenerator() : probabilityDistribution(0.0, 1.0) {
     std::random_device rd;
     engine = std::mt19937(rd());
 }
 
-RandomProbabilityGenerator<double>::RandomProbabilityGenerator(uint64_t seed) : distribution(0.0, 1.0), engine(seed) {}
+RandomProbabilityGenerator<double>::RandomProbabilityGenerator(uint64_t seed) : engine(seed), probabilityDistribution(0.0, 1.0) {}
+
+double RandomProbabilityGenerator<double>::randomProbability() {
+    return probabilityDistribution(engine);
+}
 
-double RandomProbabilityGenerator<double>::random() {
-    return distribution(engine);
+uint64_t RandomProbabilityGenerator<double>::randomSelect(uint64_t min, uint64_t max) {
+    decltype(uniformDistribution.param()) newInterval(min, max);
+    uniformDistribution.param(newInterval);
+    return uniformDistribution(engine);
+    // return std::uniform_int_distribution<uint64_t>(min, max)(engine);
 }
 
-uint64_t RandomProbabilityGenerator<double>::random_uint(uint64_t min, uint64_t max) {
-    return std::uniform_int_distribution<uint64_t>(min, max)(engine);
+double RandomProbabilityGenerator<double>::randomExponential(double rate) {
+    decltype(exponentialDistribution.param()) newRate(rate);
+    exponentialDistribution.param(newRate);
+    return exponentialDistribution(engine);
 }
 
-RandomProbabilityGenerator<RationalNumber>::RandomProbabilityGenerator() : distribution(0, std::numeric_limits<uint64_t>::max()) {
+RandomProbabilityGenerator<storm::RationalNumber>::RandomProbabilityGenerator() : probabilityDistribution(0, std::numeric_limits<uint64_t>::max()) {
     std::random_device rd;
     engine = std::mt19937(rd());
 }
 
-RandomProbabilityGenerator<RationalNumber>::RandomProbabilityGenerator(uint64_t seed) : distribution(0, std::numeric_limits<uint64_t>::max()), engine(seed) {}
+RandomProbabilityGenerator<storm::RationalNumber>::RandomProbabilityGenerator(uint64_t seed)
+    : engine(seed), probabilityDistribution(0, std::numeric_limits<uint64_t>::max()) {}
 
-RationalNumber RandomProbabilityGenerator<RationalNumber>::random() {
-    return carl::rationalize<RationalNumber>(distribution(engine)) / carl::rationalize<RationalNumber>(std::numeric_limits<uint64_t>::max());
+storm::RationalNumber RandomProbabilityGenerator<storm::RationalNumber>::randomProbability() {
+    return carl::rationalize<storm::RationalNumber>(probabilityDistribution(engine)) /
+           carl::rationalize<storm::RationalNumber>(std::numeric_limits<uint64_t>::max());
 }
 
-uint64_t RandomProbabilityGenerator<RationalNumber>::random_uint(uint64_t min, uint64_t max) {
-    return std::uniform_int_distribution<uint64_t>(min, max)(engine);
+uint64_t RandomProbabilityGenerator<storm::RationalNumber>::randomSelect(uint64_t min, uint64_t max) {
+    decltype(uniformDistribution.param()) newInterval(min, max);
+    uniformDistribution.param(newInterval);
+    return uniformDistribution(engine);
+    // return std::uniform_int_distribution<uint64_t>(min, max)(engine);
 }
 
-BernoulliDistributionGenerator::BernoulliDistributionGenerator(double prob) : distribution(prob) {}
-
-bool BernoulliDistributionGenerator::random(boost::mt19937& engine) {
-    return distribution(engine);
+storm::RationalNumber RandomProbabilityGenerator<storm::RationalNumber>::randomExponential(storm::RationalNumber rate) {
+    STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Exponential distribution is not supported for rational numbers");
 }
 
-ExponentialDistributionGenerator::ExponentialDistributionGenerator(double rate) : distribution(rate) {}
-
-double ExponentialDistributionGenerator::random(boost::mt19937& engine) {
-    return distribution(engine);
-}
 }  // namespace utility
 }  // namespace storm
diff --git a/src/storm/utility/random.h b/src/storm/utility/random.h
index fb4f16e774..498caf5731 100644
--- a/src/storm/utility/random.h
+++ b/src/storm/utility/random.h
@@ -1,18 +1,24 @@
 #pragma once
 
-#include <boost/random.hpp>
 #include <random>
+
 #include "storm/adapters/RationalNumberAdapter.h"
 
 namespace storm {
 namespace utility {
+
 template<typename ValueType>
 class RandomProbabilityGenerator {
    public:
     RandomProbabilityGenerator();
     RandomProbabilityGenerator(uint64_t seed);
-    ValueType random() const;
-    uint64_t random_uint(uint64_t min, uint64_t max);
+    ValueType randomProbability();
+    uint64_t randomSelect(uint64_t min, uint64_t max);
+    ValueType randomExponential(ValueType rate);
+
+   protected:
+    std::mt19937 engine;
+    std::uniform_int_distribution<uint64_t> uniformDistribution;
 };
 
 template<>
@@ -20,12 +26,15 @@ class RandomProbabilityGenerator<double> {
    public:
     RandomProbabilityGenerator();
     RandomProbabilityGenerator(uint64_t seed);
-    double random();
-    uint64_t random_uint(uint64_t min, uint64_t max);
+    double randomProbability();
+    uint64_t randomSelect(uint64_t min, uint64_t max);
+    double randomExponential(double rate);
 
    private:
-    std::uniform_real_distribution<double> distribution;
     std::mt19937 engine;
+    std::uniform_int_distribution<uint64_t> uniformDistribution;
+    std::uniform_real_distribution<double> probabilityDistribution;
+    std::exponential_distribution<double> exponentialDistribution;
 };
 
 template<>
@@ -33,31 +42,14 @@ class RandomProbabilityGenerator<storm::RationalNumber> {
    public:
     RandomProbabilityGenerator();
     RandomProbabilityGenerator(uint64_t seed);
-    RationalNumber random();
-    uint64_t random_uint(uint64_t min, uint64_t max);
+    storm::RationalNumber randomProbability();
+    uint64_t randomSelect(uint64_t min, uint64_t max);
+    storm::RationalNumber randomExponential(storm::RationalNumber rate);
 
    private:
-    std::uniform_int_distribution<uint64_t> distribution;
     std::mt19937 engine;
+    std::uniform_int_distribution<uint64_t> uniformDistribution;
+    std::uniform_int_distribution<uint64_t> probabilityDistribution;
 };
-
-class BernoulliDistributionGenerator {
-   public:
-    BernoulliDistributionGenerator(double prob);
-    bool random(boost::mt19937& engine);
-
-   private:
-    boost::random::bernoulli_distribution<> distribution;
-};
-
-class ExponentialDistributionGenerator {
-   public:
-    ExponentialDistributionGenerator(double rate);
-    double random(boost::mt19937& engine);
-
-   private:
-    boost::random::exponential_distribution<> distribution;
-};
-
 }  // namespace utility
 }  // namespace storm
\ No newline at end of file
diff --git a/src/test/storm-dft/bdd/TestBddModularizer.cpp b/src/test/storm-dft/bdd/BddModularizerTest.cpp
similarity index 100%
rename from src/test/storm-dft/bdd/TestBddModularizer.cpp
rename to src/test/storm-dft/bdd/BddModularizerTest.cpp
diff --git a/src/test/storm-dft/bdd/TestBddVarOrdering.cpp b/src/test/storm-dft/bdd/BddVarOrderingTest.cpp
similarity index 97%
rename from src/test/storm-dft/bdd/TestBddVarOrdering.cpp
rename to src/test/storm-dft/bdd/BddVarOrderingTest.cpp
index aa7720430a..f6e9ee634d 100644
--- a/src/test/storm-dft/bdd/TestBddVarOrdering.cpp
+++ b/src/test/storm-dft/bdd/BddVarOrderingTest.cpp
@@ -7,7 +7,7 @@
 
 namespace {
 
-TEST(TestBddVarOrdering, VariableOrdering) {
+TEST(BddVarOrderingTest, VariableOrdering) {
 #ifdef STORM_HAVE_SYLVAN
     auto dft = storm::dft::api::loadDFTGalileoFile<double>(STORM_TEST_RESOURCES_DIR "/dft/bdd/AndOrTest.dft");
     auto manager{std::make_shared<storm::dft::storage::SylvanBddManager>()};
@@ -42,7 +42,7 @@ TEST(TestBddVarOrdering, VariableOrdering) {
 #endif
 }
 
-TEST(TestBddVarOrdering, OrderParser) {
+TEST(BddVarOrderingTest, OrderParser) {
 #ifdef STORM_HAVE_SYLVAN
     auto dft = storm::dft::api::loadDFTGalileoFile<double>(STORM_TEST_RESOURCES_DIR "/dft/bdd/AndOrTest.dft");
 
diff --git a/src/test/storm-dft/bdd/TestBdd.cpp b/src/test/storm-dft/bdd/SftBddTest.cpp
similarity index 100%
rename from src/test/storm-dft/bdd/TestBdd.cpp
rename to src/test/storm-dft/bdd/SftBddTest.cpp
diff --git a/src/test/storm-dft/simulator/DftSimulatorTest.cpp b/src/test/storm-dft/simulator/DftSimulatorTest.cpp
index 93beba164c..f8562c3b99 100644
--- a/src/test/storm-dft/simulator/DftSimulatorTest.cpp
+++ b/src/test/storm-dft/simulator/DftSimulatorTest.cpp
@@ -24,8 +24,7 @@ std::pair<double, double> simulateDft(std::string const& file, double timebound,
     storm::dft::storage::DFTStateGenerationInfo stateGenerationInfo(dft->buildStateGenerationInfo(symmetries));
 
     // Init random number generator
-    // storm::utility::setLogLevel(l3pp::LogLevel::TRACE);
-    boost::mt19937 gen(5u);
+    storm::utility::RandomProbabilityGenerator<double> gen(5u);
     storm::dft::simulator::DFTTraceSimulator<double> simulator(*dft, stateGenerationInfo, gen);
 
     size_t count = 0;
@@ -69,12 +68,7 @@ TEST(DftSimulatorTest, VotingUnreliability) {
     result = simulateDftProb(STORM_TEST_RESOURCES_DIR "/dft/voting3.dft", 1, 10000);
     EXPECT_NEAR(result, 0.3496529873, 0.01);
     result = simulateDftProb(STORM_TEST_RESOURCES_DIR "/dft/voting4.dft", 1, 10000);
-#if BOOST_VERSION > 106400
     EXPECT_NEAR(result, 0.693568287, 0.01);
-#else
-    // Older Boost versions yield different value
-    EXPECT_NEAR(result, 0.693568287, 0.015);
-#endif
 }
 
 TEST(DftSimulatorTest, PandUnreliability) {
diff --git a/src/test/storm-dft/simulator/DftTraceGeneratorTest.cpp b/src/test/storm-dft/simulator/DftTraceGeneratorTest.cpp
index a0feb68aa4..fc4a49a7dc 100644
--- a/src/test/storm-dft/simulator/DftTraceGeneratorTest.cpp
+++ b/src/test/storm-dft/simulator/DftTraceGeneratorTest.cpp
@@ -151,7 +151,7 @@ TYPED_TEST(DftTraceGeneratorTest, RandomStepsAnd) {
     EXPECT_EQ(this->getConfig().useSR && this->getConfig().useDC, pair.second.hasSymmetries());
 
     // Init random number generator
-    boost::mt19937 gen(5u);
+    storm::utility::RandomProbabilityGenerator<double> gen(5u);
     storm::dft::simulator::DFTTraceSimulator<double> simulator(*dft, pair.second, gen);
 
     auto state = simulator.getCurrentState();
@@ -160,19 +160,13 @@ TYPED_TEST(DftTraceGeneratorTest, RandomStepsAnd) {
     // First random step
     storm::dft::simulator::SimulationStepResult res = simulator.randomStep();
     EXPECT_EQ(res, storm::dft::simulator::SimulationStepResult::SUCCESSFUL);
-#if BOOST_VERSION > 106400
-    // Older Boost versions yield different value
-    EXPECT_NEAR(simulator.getCurrentTime(), 0.522079, 1e-6);
-#endif
+    EXPECT_NEAR(simulator.getCurrentTime(), 0.113522, 1e-6);
     state = simulator.getCurrentState();
     EXPECT_FALSE(state->hasFailed(dft->getTopLevelIndex()));
 
     res = simulator.randomStep();
     EXPECT_EQ(res, storm::dft::simulator::SimulationStepResult::SUCCESSFUL);
-#if BOOST_VERSION > 106400
-    // Older Boost versions yield different value
-    EXPECT_NEAR(simulator.getCurrentTime(), 0.522079 + 0.9497214, 1e-6);
-#endif
+    EXPECT_NEAR(simulator.getCurrentTime(), 0.113522 + 0.904286, 1e-6);
     state = simulator.getCurrentState();
     EXPECT_TRUE(state->hasFailed(dft->getTopLevelIndex()));
 }
@@ -183,7 +177,7 @@ TYPED_TEST(DftTraceGeneratorTest, Reset) {
     EXPECT_EQ(this->getConfig().useSR && this->getConfig().useDC, pair.second.hasSymmetries());
 
     // Init random number generator
-    boost::mt19937 gen(5u);
+    storm::utility::RandomProbabilityGenerator<double> gen(5u);
     storm::dft::simulator::DFTTraceSimulator<double> simulator(*dft, pair.second, gen);
 
     auto state = simulator.getCurrentState();
@@ -192,19 +186,13 @@ TYPED_TEST(DftTraceGeneratorTest, Reset) {
     // First random step
     storm::dft::simulator::SimulationStepResult res = simulator.randomStep();
     EXPECT_EQ(res, storm::dft::simulator::SimulationStepResult::SUCCESSFUL);
-#if BOOST_VERSION > 106400
-    // Older Boost versions yield different value
-    EXPECT_NEAR(simulator.getCurrentTime(), 0.522079, 1e-6);
-#endif
+    EXPECT_NEAR(simulator.getCurrentTime(), 0.113522, 1e-6);
     auto stateStep1 = simulator.getCurrentState();
     EXPECT_FALSE(stateStep1->hasFailed(dft->getTopLevelIndex()));
 
     res = simulator.randomStep();
     EXPECT_EQ(res, storm::dft::simulator::SimulationStepResult::SUCCESSFUL);
-#if BOOST_VERSION > 106400
-    // Older Boost versions yield different value
-    EXPECT_NEAR(simulator.getCurrentTime(), 0.522079 + 0.9497214, 1e-6);
-#endif
+    EXPECT_NEAR(simulator.getCurrentTime(), 0.113522 + 0.904286, 1e-6);
     state = simulator.getCurrentState();
     EXPECT_TRUE(state->hasFailed(dft->getTopLevelIndex()));
 
@@ -215,10 +203,7 @@ TYPED_TEST(DftTraceGeneratorTest, Reset) {
 
     res = simulator.randomStep();
     EXPECT_EQ(res, storm::dft::simulator::SimulationStepResult::SUCCESSFUL);
-#if BOOST_VERSION > 106400
-    // Older Boost versions yield different value
-    EXPECT_NEAR(simulator.getCurrentTime(), 0.522079 + 0.9497214 + 2.4686932, 1e-6);
-#endif
+    EXPECT_NEAR(simulator.getCurrentTime(), 0.113522 + 0.904286 + 7.769302, 1e-6);
     state = simulator.getCurrentState();
     EXPECT_TRUE(state->hasFailed(dft->getTopLevelIndex()));
 }
@@ -228,7 +213,7 @@ TYPED_TEST(DftTraceGeneratorTest, Fdep) {
     auto dft = pair.first;
 
     // Init random number generator. Will not be important as we are choosing the steps deterministically.
-    boost::mt19937 gen(5u);
+    storm::utility::RandomProbabilityGenerator<double> gen(5u);
     storm::dft::simulator::DFTTraceSimulator<double> simulator(*dft, pair.second, gen);
 
     // Start with initial state
diff --git a/src/test/storm-dft/simulator/ImportanceFunction.cpp b/src/test/storm-dft/simulator/ImportanceFunctionTest.cpp
similarity index 87%
rename from src/test/storm-dft/simulator/ImportanceFunction.cpp
rename to src/test/storm-dft/simulator/ImportanceFunctionTest.cpp
index 8a3ee078cb..3889324a3f 100644
--- a/src/test/storm-dft/simulator/ImportanceFunction.cpp
+++ b/src/test/storm-dft/simulator/ImportanceFunctionTest.cpp
@@ -28,7 +28,7 @@ TEST(ImportanceFunctionTest, RandomStepsAnd) {
     auto dft = pair.first;
 
     // Init random number generator
-    boost::mt19937 gen(5u);
+    storm::utility::RandomProbabilityGenerator<double> gen(5u);
     storm::dft::simulator::DFTTraceSimulator<double> simulator(*dft, pair.second, gen);
 
     // Init importance function
@@ -44,20 +44,14 @@ TEST(ImportanceFunctionTest, RandomStepsAnd) {
     // First random step
     storm::dft::simulator::SimulationStepResult res = simulator.randomStep();
     EXPECT_EQ(res, storm::dft::simulator::SimulationStepResult::SUCCESSFUL);
-#if BOOST_VERSION > 106400
-    // Older Boost versions yield different value
-    EXPECT_NEAR(simulator.getCurrentTime(), 0.522079, 1e-6);
-#endif
+    EXPECT_NEAR(simulator.getCurrentTime(), 0.113522, 1e-6);
     state = simulator.getCurrentState();
     EXPECT_FALSE(state->hasFailed(dft->getTopLevelIndex()));
     EXPECT_EQ(imp.getImportance(state), 1);
 
     res = simulator.randomStep();
     EXPECT_EQ(res, storm::dft::simulator::SimulationStepResult::SUCCESSFUL);
-#if BOOST_VERSION > 106400
-    // Older Boost versions yield different value
-    EXPECT_NEAR(simulator.getCurrentTime(), 0.522079 + 0.9497214, 1e-6);
-#endif
+    EXPECT_NEAR(simulator.getCurrentTime(), 0.113522 + 0.904286, 1e-6);
     state = simulator.getCurrentState();
     EXPECT_TRUE(state->hasFailed(dft->getTopLevelIndex()));
     EXPECT_EQ(imp.getImportance(state), 2);
diff --git a/src/test/storm-dft/simulator/SamplingTest.cpp b/src/test/storm-dft/simulator/SamplingTest.cpp
index 1fb4003fe4..dfe5c0cdaa 100644
--- a/src/test/storm-dft/simulator/SamplingTest.cpp
+++ b/src/test/storm-dft/simulator/SamplingTest.cpp
@@ -6,20 +6,13 @@
 namespace {
 
 TEST(SamplingTest, SampleExponential) {
-#if BOOST_VERSION < 106400
-    // Boost changed implementation of exponential distribution
-    // -> different values are returned
-    GTEST_SKIP();
-    return;
-#endif
-    boost::mt19937 gen(5u);
-    storm::utility::ExponentialDistributionGenerator dist(5);
+    storm::utility::RandomProbabilityGenerator<double> dist(5u);
 
     // Ensure that pseudo random numbers are the same on all machines
-    double reference[] = {0.18241937154,   0.0522078772595, 0.0949721368604, 0.246869315378, 0.765000791199,
-                          0.0177096648877, 0.225167598601,  0.23538530391,   1.01605360643,  0.138846355094};
+    double reference[] = {0.011352194735198223, 0.35595966617002761,  0.0904286229501646,  0.77693019063043178, 0.018822806624073801,
+                          0.10108027454097809,  0.087433897364154994, 0.13335480534208866, 0.93816592681955202, 0.33034562021306652};
     for (int i = 0; i < 10; ++i) {
-        EXPECT_NEAR(dist.random(gen), reference[i], 1e-10);
+        EXPECT_NEAR(dist.randomExponential(5), reference[i], 1e-10);
     }
 }
 
diff --git a/src/test/storm/simulator/PrismProgramSimulator.cpp b/src/test/storm/simulator/PrismProgramSimulator.cpp
deleted file mode 100644
index 50e0b0929f..0000000000
--- a/src/test/storm/simulator/PrismProgramSimulator.cpp
+++ /dev/null
@@ -1,55 +0,0 @@
-#include "storm-config.h"
-#include "test/storm_gtest.h"
-
-#include "storm-parsers/parser/PrismParser.h"
-#include "storm/environment/Environment.h"
-#include "storm/simulator/PrismProgramSimulator.h"
-
-TEST(PrismProgramSimulatorTest, KnuthYaoDieTest) {
-#ifndef STORM_HAVE_Z3
-    GTEST_SKIP() << "Z3 not available.";
-#endif
-    storm::Environment env;
-    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_TEST_RESOURCES_DIR "/mdp/die_c1.nm");
-    storm::builder::BuilderOptions options;
-    options.setBuildAllRewardModels();
-
-    storm::simulator::DiscreteTimePrismProgramSimulator<double> sim(program, options);
-    sim.setSeed(42);
-    EXPECT_EQ("coin_flips", sim.getRewardNames()[0]);
-    auto rew = sim.getLastRewards();
-    rew = sim.getLastRewards();
-    EXPECT_EQ(1ul, rew.size());
-    EXPECT_EQ(0.0, rew[0]);
-    auto labels = sim.getCurrentStateLabelling();
-    EXPECT_EQ(0ul, labels.size());
-    EXPECT_EQ(2ul, sim.getChoices().size());
-    sim.step(0);
-    rew = sim.getLastRewards();
-    EXPECT_EQ(1ul, rew.size());
-    EXPECT_EQ(0.0, rew[0]);
-    labels = sim.getCurrentStateLabelling();
-    EXPECT_EQ(0ul, labels.size());
-    EXPECT_EQ(1ul, sim.getChoices().size());
-    sim.step(0);
-    rew = sim.getLastRewards();
-    EXPECT_EQ(1ul, rew.size());
-    EXPECT_EQ(1.0, rew[0]);
-    labels = sim.getCurrentStateLabelling();
-    EXPECT_EQ(0ul, labels.size());
-    sim.step(0);
-    rew = sim.getLastRewards();
-    EXPECT_EQ(1ul, rew.size());
-    EXPECT_EQ(1.0, rew[0]);
-    labels = sim.getCurrentStateLabelling();
-    EXPECT_EQ(2ul, labels.size());
-    EXPECT_TRUE(std::count(labels.begin(), labels.end(), "done") == 1);
-    EXPECT_TRUE(std::count(labels.begin(), labels.end(), "five") == 1);
-    sim.step(0);
-    rew = sim.getLastRewards();
-    EXPECT_EQ(1ul, rew.size());
-    EXPECT_EQ(0.0, rew[0]);
-    labels = sim.getCurrentStateLabelling();
-    EXPECT_TRUE(std::count(labels.begin(), labels.end(), "done") == 1);
-    EXPECT_TRUE(std::count(labels.begin(), labels.end(), "five") == 1);
-}
diff --git a/src/test/storm/simulator/PrismProgramSimulatorTest.cpp b/src/test/storm/simulator/PrismProgramSimulatorTest.cpp
new file mode 100644
index 0000000000..c16a6b0f42
--- /dev/null
+++ b/src/test/storm/simulator/PrismProgramSimulatorTest.cpp
@@ -0,0 +1,206 @@
+#include "storm-config.h"
+#include "test/storm_gtest.h"
+
+#include "storm-parsers/parser/PrismParser.h"
+#include "storm/adapters/JsonAdapter.h"
+#include "storm/simulator/PrismProgramSimulator.h"
+
+TEST(PrismProgramSimulatorTest, KnuthYaoDieDtmc) {
+#ifndef STORM_HAVE_Z3
+    GTEST_SKIP() << "Z3 not available.";
+#endif
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_TEST_RESOURCES_DIR "/dtmc/die.pm");
+    storm::builder::BuilderOptions options;
+    options.setBuildAllRewardModels();
+    options.setBuildAllLabels();
+
+    storm::simulator::PrismProgramSimulator<double> sim(program, options);
+    sim.setSeed(42);
+    EXPECT_FALSE(sim.isContinuousTimeModel());
+    EXPECT_EQ("coin_flips", sim.getRewardNames()[0]);
+
+    EXPECT_EQ(0ul, sim.getStateAsJson()["s"]);
+    EXPECT_EQ(0ul, sim.getStateAsJson()["d"]);
+    auto rew = sim.getCurrentRewards();
+    rew = sim.getCurrentRewards();
+    EXPECT_EQ(1ul, rew.size());
+    EXPECT_EQ(0.0, rew[0]);
+    auto labels = sim.getCurrentStateLabelling();
+    EXPECT_EQ(0ul, labels.size());
+    EXPECT_EQ(1ul, sim.getCurrentNumberOfChoices());
+    EXPECT_EQ(0, sim.getCurrentTime());
+}
+
+TEST(PrismProgramSimulatorTest, KnuthYaoDieMdp) {
+#ifndef STORM_HAVE_Z3
+    GTEST_SKIP() << "Z3 not available.";
+#endif
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_TEST_RESOURCES_DIR "/mdp/die_c1.nm");
+    storm::builder::BuilderOptions options;
+    options.setBuildAllRewardModels();
+    options.setBuildAllLabels();
+
+    storm::simulator::PrismProgramSimulator<double> sim(program, options);
+    sim.setSeed(42);
+    EXPECT_FALSE(sim.isContinuousTimeModel());
+    EXPECT_EQ("coin_flips", sim.getRewardNames()[0]);
+
+    EXPECT_EQ(0ul, sim.getStateAsJson()["s"]);
+    EXPECT_EQ(0ul, sim.getStateAsJson()["d"]);
+    auto rew = sim.getCurrentRewards();
+    rew = sim.getCurrentRewards();
+    EXPECT_EQ(1ul, rew.size());
+    EXPECT_EQ(0.0, rew[0]);
+    auto labels = sim.getCurrentStateLabelling();
+    EXPECT_EQ(0ul, labels.size());
+    EXPECT_EQ(2ul, sim.getCurrentNumberOfChoices());
+
+    sim.step(0);
+    EXPECT_EQ(02ul, sim.getStateAsJson()["s"]);
+    EXPECT_EQ(0ul, sim.getStateAsJson()["d"]);
+    rew = sim.getCurrentRewards();
+    EXPECT_EQ(1ul, rew.size());
+    EXPECT_EQ(0.0, rew[0]);
+    EXPECT_EQ(0ul, sim.getCurrentStateLabelling().size());
+    EXPECT_EQ(1ul, sim.getCurrentNumberOfChoices());
+
+    sim.randomStep();
+    EXPECT_EQ(5ul, sim.getStateAsJson()["s"]);
+    EXPECT_EQ(0ul, sim.getStateAsJson()["d"]);
+    rew = sim.getCurrentRewards();
+    EXPECT_EQ(1ul, rew.size());
+    EXPECT_EQ(1.0, rew[0]);
+    EXPECT_EQ(0ul, sim.getCurrentStateLabelling().size());
+    EXPECT_EQ(1ul, sim.getCurrentNumberOfChoices());
+
+    sim.randomStep();
+    EXPECT_EQ(7ul, sim.getStateAsJson()["s"]);
+    EXPECT_EQ(5ul, sim.getStateAsJson()["d"]);
+    rew = sim.getCurrentRewards();
+    EXPECT_EQ(1ul, rew.size());
+    EXPECT_EQ(1.0, rew[0]);
+    labels = sim.getCurrentStateLabelling();
+    EXPECT_EQ(2ul, labels.size());
+    EXPECT_TRUE(labels.contains("done"));
+    EXPECT_TRUE(labels.contains("five"));
+    EXPECT_EQ(1ul, sim.getCurrentNumberOfChoices());
+
+    sim.randomStep();
+    EXPECT_EQ(7ul, sim.getStateAsJson()["s"]);
+    EXPECT_EQ(5ul, sim.getStateAsJson()["d"]);
+    rew = sim.getCurrentRewards();
+    EXPECT_EQ(1ul, rew.size());
+    EXPECT_EQ(0.0, rew[0]);
+    labels = sim.getCurrentStateLabelling();
+    EXPECT_EQ(2ul, labels.size());
+    EXPECT_TRUE(labels.contains("done"));
+    EXPECT_TRUE(labels.contains("five"));
+
+    EXPECT_EQ(0, sim.getCurrentTime());
+}
+
+TEST(PrismProgramSimulatorTest, SimpleCtmc) {
+#ifndef STORM_HAVE_Z3
+    GTEST_SKIP() << "Z3 not available.";
+#endif
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_TEST_RESOURCES_DIR "/ctmc/simple2.sm");
+    storm::builder::BuilderOptions options;
+    options.setBuildAllRewardModels();
+    options.setBuildAllLabels();
+
+    storm::simulator::PrismProgramSimulator<double> sim(program, options);
+    sim.setSeed(5);
+    EXPECT_TRUE(sim.isContinuousTimeModel());
+    EXPECT_EQ(2ul, sim.getRewardNames().size());
+
+    EXPECT_EQ(0ul, sim.getStateAsJson()["s"]);
+    EXPECT_EQ(2ul, sim.getCurrentRewards().size());
+    EXPECT_EQ(0ul, sim.getCurrentStateLabelling().size());
+    EXPECT_EQ(1ul, sim.getCurrentNumberOfChoices());
+    EXPECT_EQ(0, sim.getCurrentTime());
+
+    sim.randomStep();
+    EXPECT_EQ(2ul, sim.getStateAsJson()["s"]);
+    EXPECT_EQ(0ul, sim.getCurrentStateLabelling().size());
+    EXPECT_EQ(1ul, sim.getCurrentNumberOfChoices());
+    EXPECT_NEAR(0.007095, sim.getCurrentTime(), 1e-6);
+
+    sim.randomStep();
+    EXPECT_EQ(3ul, sim.getStateAsJson()["s"]);
+    EXPECT_EQ(0ul, sim.getCurrentStateLabelling().size());
+    EXPECT_EQ(1ul, sim.getCurrentNumberOfChoices());
+    EXPECT_NEAR(0.007095 + 0.452143, sim.getCurrentTime(), 1e-6);
+
+    sim.randomStep();
+    EXPECT_EQ(4ul, sim.getStateAsJson()["s"]);
+    EXPECT_EQ(0ul, sim.getCurrentStateLabelling().size());
+    EXPECT_EQ(1ul, sim.getCurrentNumberOfChoices());
+    EXPECT_NEAR(0.007095 + 0.452143 + 0.094114, sim.getCurrentTime(), 1e-6);
+}
+
+TEST(PrismProgramSimulatorTest, SimpleMA) {
+#ifndef STORM_HAVE_Z3
+    GTEST_SKIP() << "Z3 not available.";
+#endif
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_TEST_RESOURCES_DIR "/ma/simple.ma");
+    storm::builder::BuilderOptions options;
+    options.setBuildAllRewardModels();
+    options.setBuildAllLabels();
+
+    storm::simulator::PrismProgramSimulator<double> sim(program, options);
+    sim.setSeed(5);
+    EXPECT_TRUE(sim.isContinuousTimeModel());
+    EXPECT_EQ(0ul, sim.getRewardNames().size());
+
+    // 1st run
+    EXPECT_EQ(0ul, sim.getStateAsJson()["s"]);
+    EXPECT_EQ(0ul, sim.getCurrentRewards().size());
+    auto labels = sim.getCurrentStateLabelling();
+    EXPECT_EQ(0ul, labels.size());
+    EXPECT_EQ(2ul, sim.getCurrentNumberOfChoices());
+    EXPECT_EQ(0, sim.getCurrentTime());
+
+    sim.step(1);
+    EXPECT_EQ(0ul, sim.getStateAsJson()["s"]);
+    labels = sim.getCurrentStateLabelling();
+    EXPECT_EQ(0ul, labels.size());
+    EXPECT_EQ(2ul, sim.getCurrentNumberOfChoices());
+    EXPECT_EQ(0, sim.getCurrentTime());
+
+    sim.step(1);
+    EXPECT_EQ(2ul, sim.getStateAsJson()["s"]);
+    EXPECT_EQ(0ul, sim.getCurrentStateLabelling().size());
+    EXPECT_EQ(1ul, sim.getCurrentNumberOfChoices());
+    EXPECT_EQ(0, sim.getCurrentTime());
+
+    sim.randomStep();
+    EXPECT_EQ(4ul, sim.getStateAsJson()["s"]);
+    EXPECT_EQ(0ul, sim.getCurrentStateLabelling().size());
+    EXPECT_EQ(1ul, sim.getCurrentNumberOfChoices());
+    EXPECT_NEAR(0.037679, sim.getCurrentTime(), 1e-6);
+
+    // 2nd run
+    sim.resetToInitial();
+    EXPECT_EQ(0ul, sim.getStateAsJson()["s"]);
+    EXPECT_EQ(0, sim.getCurrentTime());
+    sim.step(0);
+    EXPECT_EQ(1ul, sim.getStateAsJson()["s"]);
+    EXPECT_EQ(0, sim.getCurrentTime());
+    sim.randomStep();
+    EXPECT_EQ(0ul, sim.getStateAsJson()["s"]);
+    EXPECT_NEAR(0.388465, sim.getCurrentTime(), 1e-6);
+
+    sim.step(0);
+    EXPECT_EQ(1ul, sim.getStateAsJson()["s"]);
+    EXPECT_NEAR(0.388465, sim.getCurrentTime(), 1e-6);
+    sim.randomStep();
+    EXPECT_EQ(0ul, sim.getStateAsJson()["s"]);
+    EXPECT_NEAR(0.388465 + 0.050540, sim.getCurrentTime(), 1e-6);
+
+    sim.step(0);
+    EXPECT_EQ(1ul, sim.getStateAsJson()["s"]);
+    EXPECT_NEAR(0.388465 + 0.050540, sim.getCurrentTime(), 1e-6);
+    sim.randomStep();
+    EXPECT_EQ(3ul, sim.getStateAsJson()["s"]);
+    EXPECT_NEAR(0.388465 + 0.050540 + 0.066677, sim.getCurrentTime(), 1e-6);
+}
diff --git a/src/test/storm/simulator/SparseModelSimulatorTest.cpp b/src/test/storm/simulator/SparseModelSimulatorTest.cpp
new file mode 100644
index 0000000000..651bd6eaef
--- /dev/null
+++ b/src/test/storm/simulator/SparseModelSimulatorTest.cpp
@@ -0,0 +1,209 @@
+#include "storm-config.h"
+#include "test/storm_gtest.h"
+
+#include "storm-parsers/parser/PrismParser.h"
+#include "storm/api/builder.h"
+#include "storm/simulator/SparseModelSimulator.h"
+
+TEST(SparseModelSimulatorTest, KnuthYaoDieDtmc) {
+#ifndef STORM_HAVE_Z3
+    GTEST_SKIP() << "Z3 not available.";
+#endif
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_TEST_RESOURCES_DIR "/dtmc/die.pm");
+    storm::builder::BuilderOptions options;
+    options.setBuildAllRewardModels();
+    options.setBuildAllLabels();
+    auto model = storm::api::buildSparseModel<double>(program, options)->template as<storm::models::sparse::Dtmc<double>>();
+
+    storm::simulator::SparseModelSimulator<double> sim(model);
+    sim.setSeed(42);
+    EXPECT_FALSE(sim.isContinuousTimeModel());
+    EXPECT_EQ("coin_flips", sim.getRewardNames()[0]);
+
+    EXPECT_EQ(0ul, sim.getCurrentState());
+    auto rew = sim.getCurrentRewards();
+    rew = sim.getCurrentRewards();
+    EXPECT_EQ(1ul, rew.size());
+    EXPECT_EQ(0.0, rew[0]);
+    auto labels = sim.getCurrentStateLabelling();
+    EXPECT_EQ(1ul, labels.size());
+    EXPECT_EQ("init", *labels.begin());
+    EXPECT_EQ(1ul, sim.getCurrentNumberOfChoices());
+    EXPECT_EQ(0, sim.getCurrentTime());
+}
+
+TEST(SparseModelSimulatorTest, KnuthYaoDieMdp) {
+#ifndef STORM_HAVE_Z3
+    GTEST_SKIP() << "Z3 not available.";
+#endif
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_TEST_RESOURCES_DIR "/mdp/die_c1.nm");
+    storm::builder::BuilderOptions options;
+    options.setBuildAllRewardModels();
+    options.setBuildAllLabels();
+    auto model = storm::api::buildSparseModel<double>(program, options)->template as<storm::models::sparse::Mdp<double>>();
+
+    storm::simulator::SparseModelSimulator<double> sim(model);
+    sim.setSeed(42);
+    EXPECT_FALSE(sim.isContinuousTimeModel());
+    EXPECT_EQ("coin_flips", sim.getRewardNames()[0]);
+    EXPECT_EQ(0ul, sim.getCurrentState());
+    auto rew = sim.getCurrentRewards();
+    rew = sim.getCurrentRewards();
+    EXPECT_EQ(1ul, rew.size());
+    EXPECT_EQ(0.0, rew[0]);
+    auto labels = sim.getCurrentStateLabelling();
+    EXPECT_EQ(1ul, labels.size());
+    EXPECT_EQ("init", *labels.begin());
+    EXPECT_EQ(2ul, sim.getCurrentNumberOfChoices());
+
+    sim.step(0);
+    EXPECT_EQ(2ul, sim.getCurrentState());
+    rew = sim.getCurrentRewards();
+    EXPECT_EQ(1ul, rew.size());
+    EXPECT_EQ(0.0, rew[0]);
+    EXPECT_EQ(0ul, sim.getCurrentStateLabelling().size());
+    EXPECT_EQ(1ul, sim.getCurrentNumberOfChoices());
+
+    sim.randomStep();
+    EXPECT_EQ(5ul, sim.getCurrentState());
+    rew = sim.getCurrentRewards();
+    EXPECT_EQ(1ul, rew.size());
+    EXPECT_EQ(1.0, rew[0]);
+    EXPECT_EQ(0ul, sim.getCurrentStateLabelling().size());
+    EXPECT_EQ(1ul, sim.getCurrentNumberOfChoices());
+
+    sim.randomStep();
+    EXPECT_EQ(11ul, sim.getCurrentState());
+    rew = sim.getCurrentRewards();
+    EXPECT_EQ(1ul, rew.size());
+    EXPECT_EQ(1.0, rew[0]);
+    labels = sim.getCurrentStateLabelling();
+    EXPECT_EQ(2ul, labels.size());
+    EXPECT_TRUE(labels.contains("done"));
+    EXPECT_TRUE(labels.contains("five"));
+    EXPECT_EQ(1ul, sim.getCurrentNumberOfChoices());
+
+    sim.randomStep();
+    EXPECT_EQ(11ul, sim.getCurrentState());
+    rew = sim.getCurrentRewards();
+    EXPECT_EQ(1ul, rew.size());
+    EXPECT_EQ(0.0, rew[0]);
+    labels = sim.getCurrentStateLabelling();
+    EXPECT_EQ(2ul, labels.size());
+    EXPECT_TRUE(labels.contains("done"));
+    EXPECT_TRUE(labels.contains("five"));
+
+    EXPECT_EQ(0, sim.getCurrentTime());
+}
+
+TEST(SparseModelSimulatorTest, SimpleCtmc) {
+#ifndef STORM_HAVE_Z3
+    GTEST_SKIP() << "Z3 not available.";
+#endif
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_TEST_RESOURCES_DIR "/ctmc/simple2.sm");
+    storm::builder::BuilderOptions options;
+    options.setBuildAllRewardModels();
+    options.setBuildAllLabels();
+    auto model = storm::api::buildSparseModel<double>(program, options)->template as<storm::models::sparse::Ctmc<double>>();
+
+    storm::simulator::SparseModelSimulator<double> sim(model);
+    sim.setSeed(5);
+    EXPECT_TRUE(sim.isContinuousTimeModel());
+    EXPECT_EQ(2ul, sim.getRewardNames().size());
+
+    EXPECT_EQ(0ul, sim.getCurrentState());
+    EXPECT_EQ(2ul, sim.getCurrentRewards().size());
+    auto labels = sim.getCurrentStateLabelling();
+    EXPECT_EQ(1ul, labels.size());
+    EXPECT_EQ("init", *labels.begin());
+    EXPECT_EQ(1ul, sim.getCurrentNumberOfChoices());
+    EXPECT_EQ(0, sim.getCurrentTime());
+
+    sim.randomStep();
+    EXPECT_EQ(2ul, sim.getCurrentState());
+    EXPECT_EQ(0ul, sim.getCurrentStateLabelling().size());
+    EXPECT_EQ(1ul, sim.getCurrentNumberOfChoices());
+    EXPECT_NEAR(0.007095, sim.getCurrentTime(), 1e-6);
+
+    sim.randomStep();
+    EXPECT_EQ(3ul, sim.getCurrentState());
+    EXPECT_EQ(0ul, sim.getCurrentStateLabelling().size());
+    EXPECT_EQ(1ul, sim.getCurrentNumberOfChoices());
+    EXPECT_NEAR(0.007095 + 0.452143, sim.getCurrentTime(), 1e-6);
+
+    sim.randomStep();
+    EXPECT_EQ(4ul, sim.getCurrentState());
+    EXPECT_EQ(0ul, sim.getCurrentStateLabelling().size());
+    EXPECT_EQ(1ul, sim.getCurrentNumberOfChoices());
+    EXPECT_NEAR(0.007095 + 0.452143 + 0.094114, sim.getCurrentTime(), 1e-6);
+}
+
+TEST(SparseModelSimulatorTest, SimpleMA) {
+#ifndef STORM_HAVE_Z3
+    GTEST_SKIP() << "Z3 not available.";
+#endif
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_TEST_RESOURCES_DIR "/ma/simple.ma");
+    storm::builder::BuilderOptions options;
+    options.setBuildAllRewardModels();
+    options.setBuildAllLabels();
+    auto model = storm::api::buildSparseModel<double>(program, options)->template as<storm::models::sparse::MarkovAutomaton<double>>();
+
+    storm::simulator::SparseModelSimulator<double> sim(model);
+    sim.setSeed(5);
+    EXPECT_TRUE(sim.isContinuousTimeModel());
+    EXPECT_EQ(0ul, sim.getRewardNames().size());
+
+    // 1st run
+    EXPECT_EQ(0ul, sim.getCurrentState());
+    EXPECT_EQ(0ul, sim.getCurrentRewards().size());
+    auto labels = sim.getCurrentStateLabelling();
+    EXPECT_EQ(1ul, labels.size());
+    EXPECT_EQ("init", *labels.begin());
+    EXPECT_EQ(2ul, sim.getCurrentNumberOfChoices());
+    EXPECT_EQ(0, sim.getCurrentTime());
+
+    sim.step(1);
+    EXPECT_EQ(0ul, sim.getCurrentState());
+    labels = sim.getCurrentStateLabelling();
+    EXPECT_EQ(1ul, labels.size());
+    EXPECT_EQ("init", *labels.begin());
+    EXPECT_EQ(2ul, sim.getCurrentNumberOfChoices());
+    EXPECT_EQ(0, sim.getCurrentTime());
+
+    sim.step(1);
+    EXPECT_EQ(2ul, sim.getCurrentState());
+    EXPECT_EQ(0ul, sim.getCurrentStateLabelling().size());
+    EXPECT_EQ(1ul, sim.getCurrentNumberOfChoices());
+    EXPECT_EQ(0, sim.getCurrentTime());
+
+    sim.randomStep();
+    EXPECT_EQ(4ul, sim.getCurrentState());
+    EXPECT_EQ(0ul, sim.getCurrentStateLabelling().size());
+    EXPECT_EQ(1ul, sim.getCurrentNumberOfChoices());
+    EXPECT_NEAR(0.037679, sim.getCurrentTime(), 1e-6);
+
+    // 2nd run
+    sim.resetToInitial();
+    EXPECT_EQ(0ul, sim.getCurrentState());
+    EXPECT_EQ(0, sim.getCurrentTime());
+    sim.step(0);
+    EXPECT_EQ(1ul, sim.getCurrentState());
+    EXPECT_EQ(0, sim.getCurrentTime());
+    sim.randomStep();
+    EXPECT_EQ(0ul, sim.getCurrentState());
+    EXPECT_NEAR(0.388465, sim.getCurrentTime(), 1e-6);
+
+    sim.step(0);
+    EXPECT_EQ(1ul, sim.getCurrentState());
+    EXPECT_NEAR(0.388465, sim.getCurrentTime(), 1e-6);
+    sim.randomStep();
+    EXPECT_EQ(0ul, sim.getCurrentState());
+    EXPECT_NEAR(0.388465 + 0.050540, sim.getCurrentTime(), 1e-6);
+
+    sim.step(0);
+    EXPECT_EQ(1ul, sim.getCurrentState());
+    EXPECT_NEAR(0.388465 + 0.050540, sim.getCurrentTime(), 1e-6);
+    sim.randomStep();
+    EXPECT_EQ(3ul, sim.getCurrentState());
+    EXPECT_NEAR(0.388465 + 0.050540 + 0.066677, sim.getCurrentTime(), 1e-6);
+}
diff --git a/src/test/storm/simulator/TraceSimulatorTest.cpp b/src/test/storm/simulator/TraceSimulatorTest.cpp
new file mode 100644
index 0000000000..10566f7bd8
--- /dev/null
+++ b/src/test/storm/simulator/TraceSimulatorTest.cpp
@@ -0,0 +1,119 @@
+#include "storm-config.h"
+#include "test/storm_gtest.h"
+
+#include "storm-parsers/api/storm-parsers.h"
+#include "storm/api/builder.h"
+#include "storm/simulator/PrismProgramSimulator.h"
+#include "storm/simulator/SparseModelSimulator.h"
+#include "storm/simulator/TraceSimulator.h"
+
+template<typename CValueType>
+class PrismProgramSimulator {
+   public:
+    typedef CValueType ValueType;
+
+    static std::shared_ptr<storm::simulator::TraceSimulator<ValueType>> createSimulator(std::string const& file, uint64_t seed) {
+        storm::prism::Program program = storm::api::parseProgram(file, true);
+        program = storm::utility::prism::preprocess(program, "");  // Needed for Tandem for some reason
+        storm::builder::BuilderOptions options;
+        options.setBuildAllRewardModels();
+        options.setBuildAllLabels();
+
+        auto prismProgramSimulator = std::make_shared<storm::simulator::PrismProgramSimulator<ValueType>>(program, options);
+        prismProgramSimulator->setSeed(seed);
+        return std::make_shared<storm::simulator::TraceSimulator<ValueType>>(prismProgramSimulator);
+    }
+};
+
+template<typename CValueType>
+class SparseModelSimulator {
+   public:
+    typedef CValueType ValueType;
+
+    static std::shared_ptr<storm::simulator::TraceSimulator<ValueType>> createSimulator(std::string const& file, uint64_t seed) {
+        storm::prism::Program program = storm::api::parseProgram(file, true);
+        storm::builder::BuilderOptions options;
+        options.setBuildAllLabels();
+        auto model = storm::api::buildSparseModel<ValueType>(program, options);
+
+        auto sparseModelSimulator = std::make_shared<storm::simulator::SparseModelSimulator<ValueType>>(model);
+        sparseModelSimulator->setSeed(seed);
+        return std::make_shared<storm::simulator::TraceSimulator<ValueType>>(sparseModelSimulator);
+    }
+};
+
+template<typename TestType>
+class TraceSimulatorTest : public ::testing::Test {
+   public:
+    typedef typename TestType::ValueType ValueType;
+
+    TraceSimulatorTest() {}
+
+    void SetUp() override {
+#ifndef STORM_HAVE_Z3
+        GTEST_SKIP() << "Z3 not available.";
+#endif
+    }
+
+    std::shared_ptr<storm::simulator::TraceSimulator<ValueType>> createSimulator(std::string const& file) {
+        return TestType::createSimulator(file, 5);
+    }
+
+    ValueType parseNumber(std::string const& input) {
+        return storm::utility::convertNumber<ValueType>(input);
+    }
+};
+
+typedef ::testing::Types<PrismProgramSimulator<double>, SparseModelSimulator<double>> TestingTypes;
+
+TYPED_TEST_SUITE(TraceSimulatorTest, TestingTypes, );
+
+TYPED_TEST(TraceSimulatorTest, KnuthYaoDieDtmc) {
+    auto simulator = this->createSimulator(STORM_TEST_RESOURCES_DIR "/dtmc/die.pm");
+
+    double result = simulator->simulateStepBoundedReachability("one", 100, 10000);
+    EXPECT_NEAR(this->parseNumber("1/6"), result, 0.01);
+
+    result = simulator->simulateStepBoundedReachability("two", 100, 10000);
+    EXPECT_NEAR(this->parseNumber("1/6"), result, 0.01);
+}
+
+TYPED_TEST(TraceSimulatorTest, KnuthYaoDiceMdp) {
+    auto simulator = this->createSimulator(STORM_TEST_RESOURCES_DIR "/mdp/two_dice.nm");
+
+    double result = simulator->simulateStepBoundedReachability("two", 100, 10000);
+    EXPECT_NEAR(this->parseNumber("1/32"), result, 0.01);
+
+    result = simulator->simulateStepBoundedReachability("three", 100, 10000);
+    EXPECT_NEAR(this->parseNumber("2/32"), result, 0.01);
+
+    result = simulator->simulateStepBoundedReachability("four", 100, 2 * 10000);
+    EXPECT_NEAR(this->parseNumber("3/32"), result, 0.01);
+}
+
+TYPED_TEST(TraceSimulatorTest, EmbeddedCtmc) {
+    auto simulator = this->createSimulator(STORM_TEST_RESOURCES_DIR "/ctmc/embedded2.sm");
+
+    double result = simulator->simulateTimeBoundedReachability("network_full", 10000, 10000);
+    EXPECT_NEAR(0.0019216435246119591, result, 0.01);
+
+    result = simulator->simulateTimeBoundedReachability("fail_io", 10000, 10000);
+    EXPECT_NEAR(0.001911229643, result, 0.01);
+}
+
+TYPED_TEST(TraceSimulatorTest, TandemCtmc) {
+    auto simulator = this->createSimulator(STORM_TEST_RESOURCES_DIR "/ctmc/tandem5.sm");
+
+    double result = simulator->simulateTimeBoundedReachability("network_full", 10, 10000);
+    EXPECT_NEAR(0.015446370562428037, result, 0.01);
+
+    result = simulator->simulateTimeBoundedReachability("first_queue_full", 10, 10000);
+    EXPECT_NEAR(0.999999837225515, result, 0.01);
+}
+
+TYPED_TEST(TraceSimulatorTest, ServerMA) {
+    auto simulator = this->createSimulator(STORM_TEST_RESOURCES_DIR "/ma/server.ma");
+
+    double result = simulator->simulateTimeBoundedReachability("error", 1, 10000);
+    EXPECT_TRUE(storm::utility::isBetween(0.382120908, result, 0.455504));
+}
diff --git a/src/test/storm/storage/JaniLocalEliminatorTests.cpp b/src/test/storm/storage/JaniLocalEliminatorTest.cpp
similarity index 100%
rename from src/test/storm/storage/JaniLocalEliminatorTests.cpp
rename to src/test/storm/storage/JaniLocalEliminatorTest.cpp