diff --git a/core/MarkovProcesses.jl b/core/MarkovProcesses.jl
index fe716b1e52d1118077504a5cb9705a42552744ea..01657933bb049cef75f4f27e419d6ff64dc049b3 100644
--- a/core/MarkovProcesses.jl
+++ b/core/MarkovProcesses.jl
@@ -29,7 +29,7 @@ export InvariantPredicateFunction, CheckConstraintsFunction, UpdateStateFunction
 # Trajectory related methods
 export +, -, δ, dist_lp, euclidean_distance
 export get_obs_var, length_states, length_obs_var
-export get_state_from_time, get_var_from_time, vectorize
+export get_state_from_time, get_var_from_time, vectorize, trajectory_from_csv
 export isbounded, states, times, transitions
 export check_consistency, issteadystate, isaccepted
 
diff --git a/core/trajectory.jl b/core/trajectory.jl
index da97f76eb0adce7c817a6565db09997f0a4b658a..371241530a42c23e10ec6a8277a465ecab0e3e88 100644
--- a/core/trajectory.jl
+++ b/core/trajectory.jl
@@ -308,3 +308,23 @@ function +(σ1::AbstractTrajectory,σ2::AbstractTrajectory) end
 function -(σ1::AbstractTrajectory,σ2::AbstractTrajectory) end
 δ(σ::AbstractTrajectory,idx::Int) = times(σ)[i+1] - times(σ)[i]
 
+function trajectory_from_csv(csv_file, model::ContinuousTimeModel)
+    csv_mat_values, header = readdlm(csv_file, ',', header = true)
+    nbr_states = size(csv_mat_values, 1)
+    times = zeros(nbr_states)
+    values = Vector{Vector{Int}}(undef, length(model.g))
+    transitions = fill(nothing, nbr_states)
+    for i = eachindex(header)
+        model_var = header[i]
+        if model_var == "time"
+            times = csv_mat_values[:,i]
+        elseif model_var == "transitions"
+            transitions = csv_mat_values[:,i]
+        else
+            @assert Symbol(model_var) in model.g "Variable is not observed in the model"
+            values[model._map_obs_var_idx[Symbol(model_var)]] = csv_mat_values[:,i]
+        end
+    end
+    return Trajectory(model, values, times, transitions)
+end
+