diff --git a/src/SpinShuttling.jl b/src/SpinShuttling.jl
index 00b4382..699e526 100644
--- a/src/SpinShuttling.jl
+++ b/src/SpinShuttling.jl
@@ -215,7 +215,7 @@ function sampling(samplingfunction::Function, M::Int)::Union{Tuple{Real,Real},Tu
     N = length(samplingfunction(1))
     A = N > 1 ? zeros(N) : 0
     Q = copy(A)
-    for k in 1:M
+    Threads.@threads for k in 1:M
         x = samplingfunction(k)::Union{Real,Vector{<:Real}}
         Q = Q +(k-1)/k*(x-A).^ 2
         A = A + (x-A)/k
diff --git a/src/stochastics.jl b/src/stochastics.jl
index 0183e7c..2b232b2 100644
--- a/src/stochastics.jl
+++ b/src/stochastics.jl
@@ -135,7 +135,7 @@ function covariancematrix(P₁::Matrix{<:Real}, P₂::Matrix{<:Real}, process::G
     N = size(P₁, 1)
     P₁ = P₁'; P₂ = P₂';
     A = Matrix{Real}(undef, N, N)
-    for i in 1:N
+    Threads.@threads for i in 1:N
         for j in 1:N
             A[i, j] = covariance(P₁[:, i], P₂[:, j], process)
         end
@@ -150,7 +150,7 @@ function covariancematrix(P::Matrix{<:Real}, process::GaussianRandomField)::Symm
     N = size(P, 1)
     P = P'
     A = Matrix{Real}(undef, N, N)
-    for i in 1:N
+    Threads.@threads for i in 1:N
         for j in i:N
             A[i, j] = covariance(P[:, i], P[:, j], process)
         end
@@ -164,6 +164,7 @@ numerical quadrature of the covariance matrix.
 Using Simpson's rule by default.
 """
 function characteristicfunction(R::RandomFunction)::Tuple{Vector{<:Real},Vector{<:Number}}
+    # need further optimization
     dt=R.P[2,1]-R.P[1,1]
     N=size(R.Σ,1)
     @assert N%2==1