diff --git a/core/MarkovProcesses.jl b/core/MarkovProcesses.jl
new file mode 100644
index 0000000000000000000000000000000000000000..e1f47fea94f4bc535978564b820ba4be79d6e39c
--- /dev/null
+++ b/core/MarkovProcesses.jl
@@ -0,0 +1,14 @@
+module MarkovProcesses
+
+import Base: +, -, getfield, getindex
+
+export Model, ContinuousTimeModel, DiscreteTimeModel
+export simulate, set_param!, get_param
+include("model.jl")
+
+export Observations, AbstractTrajectory
+export +,-,δ,get_obs_variables
+include("observations.jl")
+
+end
+
diff --git a/core/model.jl b/core/model.jl
new file mode 100644
index 0000000000000000000000000000000000000000..4c093a8e6a86923c76add5c46ed9d94e863b7b8b
--- /dev/null
+++ b/core/model.jl
@@ -0,0 +1,10 @@
+
+abstract type Model end 
+abstract type ContinuousTimeModel <: Model end 
+abstract type DiscreteTimeModel <: Model end 
+
+function check_consistency(m::Model) end
+function simulate(m::Model, n::Int; bound::Real = Inf)::AbstractObservations end
+function set_param!(m::Model, p::AbstractVector{Real})::Nothing end
+function get_param(m::Model)::AbstractVector{Real} end
+
diff --git a/core/observations.jl b/core/observations.jl
new file mode 100644
index 0000000000000000000000000000000000000000..931f6bbd21f9900f5098e609e53db7fb88dbdfd4
--- /dev/null
+++ b/core/observations.jl
@@ -0,0 +1,26 @@
+
+abstract type AbstractTrajectory end
+ContinuousObservations = AbstractVector{AbstractTrajectory}
+
+struct Trajectory <: AbstractTrajectory
+    m::ContinuousTimeModel
+    values::AbstractMatrix{Real}
+    times::AbstractMatrix{Real}
+    transitions::AbstractVector{Union{String,Missing,Nothing}}
+end
+
+struct ObservedTrajectory <: AbstractTrajectory
+    m::ContinuousTimeModel
+    values::AbstractMatrix{Real}
+    times::AbstractMatrix{Real}
+end
+
+function +(σ1::AbstractTrajectory,σ2::AbstractTrajectory) end
+function -(σ1::AbstractTrajectory,σ2::AbstractTrajectory) end
+function δ(σ1::AbstractTrajectory,t::Real) end
+function get_obs_variables(σ::Trajectory) end
+function get_obs_variables(σ::ObservedTrajectory) end
+function get_values(σ::AbstractTrajectory, variable::String) end
+function get_times(σ::AbstractTrajectory, variable::String) end
+function getindex(σ::AbstractTrajectory, idx::String) end
+
diff --git a/tests/load_module.jl b/tests/load_module.jl
new file mode 100644
index 0000000000000000000000000000000000000000..548f0bcef52722fcd973676ef4597913dfef2547
--- /dev/null
+++ b/tests/load_module.jl
@@ -0,0 +1,3 @@
+
+using MarkovProcesses
+
diff --git a/tests/models/sir.jl b/tests/models/sir.jl
new file mode 100644
index 0000000000000000000000000000000000000000..c4db3e0345ed961c5df5ebb3747056b6070eb2bc
--- /dev/null
+++ b/tests/models/sir.jl
@@ -0,0 +1,52 @@
+
+import StaticArrays: SVector, SMatrix, @SMatrix
+
+State = SVector{3, Int}
+Parameters = SVector{2, Real}
+
+d=3
+dobs=1
+k=2
+dict=Dict("S" => 1, "I" => 2, "R" => 3)
+l_name_param = ["ki", "kr"]
+p = Parameters(0.0012, 0.05)
+x0 = State(95, 5, 0)
+
+function f(xn::State, p::Parameters, tn::Real)
+    a1 = p[1] * xn[1] * xn[2]
+    a2 = p[2] * xn[2]
+    l_a = SVector(a1, a2)
+    asum = sum(l_a)
+    # column-major order
+    l_nu = @SMatrix [-1.0 0.0;
+                     1.0 -1.0;
+                     0.0 1.0]
+    
+    u1, u2 = rand(), rand()
+    tau = - log(u1) / asum
+    b_inf = 0.0
+    b_sup = a1
+    reaction = 0
+    for i = 1:2
+        if b_inf < asum*u2 < b_sup
+            reaction = i
+            break
+        end
+        b_inf += l_a[i]
+        b_sup += l_a[i+1]
+    end
+ 
+    nu = l_nu[:,reaction]
+    xnplus1 = State(xn[1]+nu[1], xn[2]+nu[2], xn[3]+nu[3])
+    tnplus1 = tn + tau
+    transition = "R$(reaction)"
+
+    return xnplus1
+end
+
+g = SVector(1)
+
+# Gamma should be constructed automatically in the case of 
+
+#m = Model(d,dobs,k,dict,l_name,param,p,x0,f,g)
+
diff --git a/tests/sim_sir.jl b/tests/sim_sir.jl
new file mode 100644
index 0000000000000000000000000000000000000000..2feb9003229886f4e3ff4271f9232a1a6d1dfd1f
--- /dev/null
+++ b/tests/sim_sir.jl
@@ -0,0 +1,12 @@
+
+using MarkovProcesses 
+using PyPlot
+
+include("models/sir.jl")
+
+σ = simulate(SIR)
+plt.figure()
+plot(σ["S,times"], σ["S,values"])
+plt.savefig("sim_sir.png")
+plt.close()
+