diff --git a/emper/Runtime.cpp b/emper/Runtime.cpp
index 4b6b00800bbbd8bb16e61f15b47c35b8c5a1bb86..c9d73304bb1b305d9aeaed4433d2145e23e099ce 100644
--- a/emper/Runtime.cpp
+++ b/emper/Runtime.cpp
@@ -12,6 +12,7 @@
 
 #include <cstdlib>	// for rand, srand, abort
 #include <cstring>
+#include <fstream>	// IWYU pragma: keep
 #include <iostream>
 #include <memory>	 // for __shared_ptr_access, shared_ptr
 #include <string>	 // for string
@@ -199,14 +200,13 @@ Runtime::~Runtime() {
 	// Print the stats after all worker threads have terminated and before we delete
 	// objects bound to the runtime's lifetime
 	if constexpr (emper::STATS) {
-		printLastRuntimeStats();
-
-		// TODO: The following should probably be part of printLastRuntimeStats?
-		if constexpr (emper::IO) {
-			emper::io::Stats::printStats(globalIo, ioContexts);
+		const char* envValueStatsFile = std::getenv("EMPER_STATS_FILE");
+		if (envValueStatsFile) {
+			std::ofstream out(envValueStatsFile, std::ofstream::out);
+			printStats(out);
+		} else {
+			printStats(std::cout);
 		}
-
-		workerSleepStrategy.printStats();
 	}
 
 	for (unsigned int i = 0; i < workerCount; ++i) {
@@ -385,9 +385,9 @@ void Runtime::initiateAndWaitUntilTermination() {
 	waitUntilFinished();
 }
 
-void Runtime::printStats(bool detailed) {
+void Runtime::printStats(std::ostream& out, bool detailed) {
 	if constexpr (!emper::STATS) {
-		std::cout << "printStats(): EMPER stats disabled by static configuration" << std::endl;
+		out << "printStats(): EMPER stats disabled by static configuration" << std::endl;
 		return;
 	}
 
@@ -410,7 +410,7 @@ void Runtime::printStats(bool detailed) {
 		stats = strategy->getStats();
 	}
 
-	stats->print();
+	stats->print(out);
 
 	emper::stats::Worker globalWorkerStats(workerCount);
 	auto workerStatsSnapshot = workerStats->getSnapshot();
@@ -418,24 +418,30 @@ void Runtime::printStats(bool detailed) {
 		globalWorkerStats += stats;
 	}
 
-	globalWorkerStats.print();
+	globalWorkerStats.print(out);
 	if (detailed) {
 		for (auto& stats : workerStatsSnapshot) {
 			stats.print();
 		}
 	}
 
-	fromAnywhereStats->print();
+	fromAnywhereStats->print(out);
+
+	if constexpr (emper::IO) {
+		emper::io::Stats::printStats(globalIo, ioContexts, out);
+	}
+
+	workerSleepStrategy.printStats(out);
 }
 
-void Runtime::printLastRuntimeStats() {
+void Runtime::printLastRuntimeStats(std::ostream& out) {
 	std::lock_guard<std::mutex> lock(currentRuntimeMutex);
 	if (currentRuntime == nullptr) {
 		LOGE("Could not print current runtime stats");
 		return;
 	}
 
-	currentRuntime->printStats();
+	currentRuntime->printStats(out);
 }
 
 void Runtime::executeAndWait(std::function<void()> f) {
diff --git a/emper/Runtime.hpp b/emper/Runtime.hpp
index f80e7a43509c4e2a13357318af42c9c2bf0c5853..18a16649759adb38e21071bea2f504cae6eda21f 100644
--- a/emper/Runtime.hpp
+++ b/emper/Runtime.hpp
@@ -9,6 +9,7 @@
 #include <cstdint>		 // for intptr_t
 #include <cstdlib>		 // for abort
 #include <functional>	 // for function
+#include <iostream>
 #include <memory>
 #include <mutex>
 #include <optional>
@@ -101,7 +102,7 @@ class Runtime : public Logger<LogSubsystem::RUNTI> {
 
 	static RuntimeStrategyFactory& DEFAULT_STRATEGY;
 
-	static void printLastRuntimeStats();
+	static void printLastRuntimeStats(std::ostream& out = std::cout);
 
 	static auto getDefaultWorkerCount() -> workerid_t;
 
@@ -206,7 +207,7 @@ class Runtime : public Logger<LogSubsystem::RUNTI> {
 
 	void initiateAndWaitUntilTermination();
 
-	void printStats(bool detailed = false);
+	void printStats(std::ostream& out = std::cout, bool detailed = false);
 
 	static auto inRuntime() -> bool { return Worker::isWorkerThread(); }
 
diff --git a/emper/RuntimeStrategyStats.hpp b/emper/RuntimeStrategyStats.hpp
index fcd66652f42168d0f10a0505e2e58308db1a183a..784332c8f2feab9b686d4b35f61ed6d801535545 100644
--- a/emper/RuntimeStrategyStats.hpp
+++ b/emper/RuntimeStrategyStats.hpp
@@ -5,5 +5,5 @@
 class RuntimeStrategyStats {
  public:
 	virtual ~RuntimeStrategyStats() = default;
-	virtual void print() = 0;
+	virtual void print(std::ostream& out) = 0;
 };
diff --git a/emper/io/Stats.cpp b/emper/io/Stats.cpp
index d1f76561ab59cda60c96c5b02905b165a9e8bef6..627074af94bd8fec39cde22b756205d95011f305 100644
--- a/emper/io/Stats.cpp
+++ b/emper/io/Stats.cpp
@@ -117,10 +117,10 @@ auto operator<<(std::ostream& os, const Stats& s) -> std::ostream& {
 	return os;
 }
 
-void Stats::printStats(IoContext* globalIoContext,
-											 const std::vector<IoContext*>& workerIoContexts) {
+void Stats::printStats(IoContext* globalIoContext, const std::vector<IoContext*>& workerIoContexts,
+											 std::ostream& out) {
 	if (globalIoContext) {
-		std::cout << globalIoContext->getStats();
+		out << globalIoContext->getStats();
 	}
 
 	// Use a stats object to calculate the averages
@@ -156,7 +156,7 @@ void Stats::printStats(IoContext* globalIoContext,
 	}
 
 	// print avgs
-	std::cout << avgs << std::endl;
+	out << avgs << std::endl;
 }
 
 }	 // namespace emper::io
diff --git a/emper/io/Stats.hpp b/emper/io/Stats.hpp
index bf303b68028ee2666fd34983604c5a55a1b2be83..8a104ce3775657dd2aaaa20f0b4b82f45c7d9a2b 100644
--- a/emper/io/Stats.hpp
+++ b/emper/io/Stats.hpp
@@ -128,7 +128,8 @@ class Stats {
 	boost::circular_buffer<unsigned> reReapCount_last;
 
 	static void printStats(IoContext* globalIoContext,
-												 const std::vector<IoContext*>& workerIoContexts);
+												 const std::vector<IoContext*>& workerIoContexts,
+												 std::ostream& out = std::cout);
 
 	Stats() : reReapCount_last(10) {}
 
diff --git a/emper/sleep_strategy/AbstractWorkerSleepStrategy.hpp b/emper/sleep_strategy/AbstractWorkerSleepStrategy.hpp
index 33b7a7ae5774b5a8096b1a5be8029f10c50e7c4d..47a8e6e69002ab467dd2ac28dc0fe7ccf91f943b 100644
--- a/emper/sleep_strategy/AbstractWorkerSleepStrategy.hpp
+++ b/emper/sleep_strategy/AbstractWorkerSleepStrategy.hpp
@@ -45,6 +45,6 @@ class AbstractWorkerSleepStrategy {
 
 	void sleep() { static_cast<T*>(this)->sleep(); }
 
-	inline void printStats() { static_cast<T*>(this)->printStats(); }
+	inline void printStats(std::ostream& out = std::cout) { static_cast<T*>(this)->printStats(out); }
 };
 }	 // namespace emper::sleep_strategy
diff --git a/emper/strategies/AbstractWorkStealingStats.cpp b/emper/strategies/AbstractWorkStealingStats.cpp
index c6e6e79620c715b91139fe944ce98a5d56345591..6b53c066cf5f730d488edc7dc231236429e59999 100644
--- a/emper/strategies/AbstractWorkStealingStats.cpp
+++ b/emper/strategies/AbstractWorkStealingStats.cpp
@@ -16,28 +16,28 @@ AbstractWorkStealingStats::AbstractWorkStealingStats(AbstractWorkStealingStrateg
 	}
 }
 
-void AbstractWorkStealingStats::print() {
-	std::cout << "total-scheduled-fibers-to-local: "
-						<< std::to_string(comulatedWorkerStats.scheduledFibersToLocal) << std::endl
-						<< "total-scheduled-fibers-to-overflow-queue: "
-						<< std::to_string(comulatedWorkerStats.scheduledFibersToOverflowQueue) << std::endl
-						<< "global-max-queue-length: " << std::to_string(comulatedWorkerStats.maxQueueLength)
-						<< std::endl
-						<< "total-next-fiber-from-local: "
-						<< std::to_string(comulatedWorkerStats.nextFiberFromLocal) << std::endl
-						<< "total-next-fiber-hint-local: "
-						<< std::to_string(comulatedWorkerStats.nextFiberFromHintLocal) << std::endl
-						<< "total-next-fiber-hint-anywherequeue: "
-						<< std::to_string(comulatedWorkerStats.nextFiberFromHintAnywhere) << std::endl
-						<< "total-next-fiber-stolen: " << std::to_string(comulatedWorkerStats.nextFiberStolen)
-						<< std::endl
-						<< "total-next-fiber-from-anywhere-queue: "
-						<< std::to_string(comulatedWorkerStats.nextFiberFromAnywhereQueue) << std::endl
-						<< "total-fibers-lifted-from-anywhere-queue: "
-						<< std::to_string(comulatedWorkerStats.totalFibersLiftedFromAnywhereQueue) << std::endl
-						<< "max-fibers-lifted-from-anywhere-queue: "
-						<< std::to_string(comulatedWorkerStats.maxFibersLiftedFromAnywhereQueue) << std::endl
-						<< "avg-fibers-lifted-from-anywhere-queue: "
-						<< std::to_string(comulatedWorkerStats.avgFibersLiftedFromAnywhereQueue.getAverage())
-						<< std::endl;
+void AbstractWorkStealingStats::print(std::ostream& out) {
+	out << "total-scheduled-fibers-to-local: "
+			<< std::to_string(comulatedWorkerStats.scheduledFibersToLocal) << std::endl
+			<< "total-scheduled-fibers-to-overflow-queue: "
+			<< std::to_string(comulatedWorkerStats.scheduledFibersToOverflowQueue) << std::endl
+			<< "global-max-queue-length: " << std::to_string(comulatedWorkerStats.maxQueueLength)
+			<< std::endl
+			<< "total-next-fiber-from-local: " << std::to_string(comulatedWorkerStats.nextFiberFromLocal)
+			<< std::endl
+			<< "total-next-fiber-hint-local: "
+			<< std::to_string(comulatedWorkerStats.nextFiberFromHintLocal) << std::endl
+			<< "total-next-fiber-hint-anywherequeue: "
+			<< std::to_string(comulatedWorkerStats.nextFiberFromHintAnywhere) << std::endl
+			<< "total-next-fiber-stolen: " << std::to_string(comulatedWorkerStats.nextFiberStolen)
+			<< std::endl
+			<< "total-next-fiber-from-anywhere-queue: "
+			<< std::to_string(comulatedWorkerStats.nextFiberFromAnywhereQueue) << std::endl
+			<< "total-fibers-lifted-from-anywhere-queue: "
+			<< std::to_string(comulatedWorkerStats.totalFibersLiftedFromAnywhereQueue) << std::endl
+			<< "max-fibers-lifted-from-anywhere-queue: "
+			<< std::to_string(comulatedWorkerStats.maxFibersLiftedFromAnywhereQueue) << std::endl
+			<< "avg-fibers-lifted-from-anywhere-queue: "
+			<< std::to_string(comulatedWorkerStats.avgFibersLiftedFromAnywhereQueue.getAverage())
+			<< std::endl;
 }
diff --git a/emper/strategies/AbstractWorkStealingStats.hpp b/emper/strategies/AbstractWorkStealingStats.hpp
index f8d19002993b517860908451ea91db41838e2418..fee4c6f7d814eac3da8cef334115301e9b370c9e 100644
--- a/emper/strategies/AbstractWorkStealingStats.hpp
+++ b/emper/strategies/AbstractWorkStealingStats.hpp
@@ -2,6 +2,7 @@
 // Copyright © 2021 Florian Schmaus
 #pragma once
 
+#include <iostream>
 #include <vector>
 
 #include "RuntimeStrategyStats.hpp"
@@ -14,7 +15,7 @@ class AbstractWorkStealingStats : public RuntimeStrategyStats {
 	std::vector<AbstractWorkStealingWorkerStats> workerStats;
 	AbstractWorkStealingWorkerStats comulatedWorkerStats;
 
-	AbstractWorkStealingStats(AbstractWorkStealingStrategy &strategy);
+	AbstractWorkStealingStats(AbstractWorkStealingStrategy& strategy);
 
-	void print() override;
+	void print(std::ostream& out = std::cout) override;
 };
diff --git a/emper/strategies/laws/LawsStrategyStats.cpp b/emper/strategies/laws/LawsStrategyStats.cpp
index 7d718508c7edc42b6e79d1b21596c86a60fc4da9..cf9e534d856e5d910fa81c473e908ca448a8f500 100644
--- a/emper/strategies/laws/LawsStrategyStats.cpp
+++ b/emper/strategies/laws/LawsStrategyStats.cpp
@@ -16,21 +16,21 @@ LawsStrategyStats::LawsStrategyStats(LawsStrategy& lawsStrategy)
 	}
 }
 
-void LawsStrategyStats::print() {
-	AbstractWorkStealingStats::print();
+void LawsStrategyStats::print(std::ostream& out) {
+	AbstractWorkStealingStats::print(out);
 
-	std::cout << "total-scheduled-fibers-to-priority: "
-						<< std::to_string(comulatedWorkerStats.scheduledFibersToPriority) << std::endl
-						<< "total-dispatched-fibers-from-priority: "
-						<< std::to_string(comulatedWorkerStats.dispatchedFibersFromPriority) << std::endl
-						<< "total-dispatched-fibers-from-local: "
-						<< std::to_string(comulatedWorkerStats.dispatchedFibersFromLocal) << std::endl
-						<< "total-dispatched-fibers-from-hint-local: "
-						<< std::to_string(comulatedWorkerStats.dispatchedFibersFromHintLocal) << std::endl
-						<< "total-dispatched-fibers-from-hint-anywhere: "
-						<< std::to_string(comulatedWorkerStats.dispatchedFibersFromHintAnywhere) << std::endl
-						<< "total-dispatched-fibers-stolen: "
-						<< std::to_string(comulatedWorkerStats.dispatchedFibersStolen) << std::endl
-						<< "total-dispatched-fibers-from-anywhere-queue: "
-						<< std::to_string(comulatedWorkerStats.dispatchedFibersFromAnywhereQueue) << std::endl;
+	out << "total-scheduled-fibers-to-priority: "
+			<< std::to_string(comulatedWorkerStats.scheduledFibersToPriority) << std::endl
+			<< "total-dispatched-fibers-from-priority: "
+			<< std::to_string(comulatedWorkerStats.dispatchedFibersFromPriority) << std::endl
+			<< "total-dispatched-fibers-from-local: "
+			<< std::to_string(comulatedWorkerStats.dispatchedFibersFromLocal) << std::endl
+			<< "total-dispatched-fibers-from-hint-local: "
+			<< std::to_string(comulatedWorkerStats.dispatchedFibersFromHintLocal) << std::endl
+			<< "total-dispatched-fibers-from-hint-anywhere: "
+			<< std::to_string(comulatedWorkerStats.dispatchedFibersFromHintAnywhere) << std::endl
+			<< "total-dispatched-fibers-stolen: "
+			<< std::to_string(comulatedWorkerStats.dispatchedFibersStolen) << std::endl
+			<< "total-dispatched-fibers-from-anywhere-queue: "
+			<< std::to_string(comulatedWorkerStats.dispatchedFibersFromAnywhereQueue) << std::endl;
 }
diff --git a/emper/strategies/laws/LawsStrategyStats.hpp b/emper/strategies/laws/LawsStrategyStats.hpp
index 35c30d9feb9297a5229efd1e3ed40ac069f2704a..3baedc2a484c6c804c8f5817d2029d2631f1a569 100644
--- a/emper/strategies/laws/LawsStrategyStats.hpp
+++ b/emper/strategies/laws/LawsStrategyStats.hpp
@@ -2,6 +2,7 @@
 // Copyright © 2020-2021 Florian Schmaus
 #pragma once
 
+#include <iostream>
 #include <vector>
 
 #include "strategies/AbstractWorkStealingStats.hpp"
@@ -16,5 +17,5 @@ class LawsStrategyStats : public AbstractWorkStealingStats {
 
 	LawsStrategyStats(LawsStrategy& lawsStrategy);
 
-	void print() override;
+	void print(std::ostream& out = std::cout) override;
 };