diff --git a/Project.toml b/Project.toml
index a77faaae..28066a39 100644
--- a/Project.toml
+++ b/Project.toml
@@ -17,6 +17,7 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsAPI = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 
 [weakdeps]
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
diff --git a/docs/make.jl b/docs/make.jl
index 6c5d8c16..52adc5ac 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -33,6 +33,7 @@ makedocs(
             "Muti-attribute attentional drift diffusion Model" => "maaDDM.md",
             "Poisson Race" => "poisson_race.md",
             "Racing Diffusion Model (RDM)" => "rdm.md",
+            "Starting-time Drift Diffusion Model (stDDM)" => "stDDM.md",
             "Wald Model" => "wald.md",
             "Wald Mixture Model" => "wald_mixture.md",
         ],
diff --git a/docs/src/stDDM.md b/docs/src/stDDM.md
new file mode 100644
index 00000000..6ca7dd81
--- /dev/null
+++ b/docs/src/stDDM.md
@@ -0,0 +1,100 @@
+# Starting-time Drift Diffusion Model (stDDM)
+
+The relative starting time drift diffusion model (stDDM) characterizes the contributions of multiple unique attributes to the rate of evidence accumulation. Compared to the DDM, which assumes a constant evidence accumulation rate within each trial, the stDDM allows different attributes to enter the evidence accumulation process at various time points relative to one another. By doing so, the stDDM quantifies both the weights given to each attribute and their onset times (Amasino et al., 2019; Barakchian et al., 2021; Chen et al., 2022; Maier et al., 2020; Sullivan and Huettel, 2021).
+
+# Example
+In this example, we will demonstrate how to use the stDDM in a generic two-alternative forced-choice task with two arbitrary attributes.
+
+## Load Packages
+The first step is to load the required packages.
+
+```@example stDDM
+using SequentialSamplingModels
+using Plots
+using Random
+
+Random.seed!(8741)
+```
+
+## Create Model Object
+In the code below, we will define parameters for the stDDM and create a model object to store the parameter values. 
+
+### Drift Rate
+The drift rate controls the speed and direction in which information accumulates. Here, each drift coefficient indicates the weighting strengths given to the first and second attributes (e.g., taste and health, payoff and delay, self and other), respectively, to the total drift rate in a given trial, where the drift rate accumulates relative evidence in favor of an option.
+```@example stDDM
+ν = [2.5,2.0]
+```
+### Diffusion Noise
+Diffusion noise is the amount of within trial noise in the evidence accumulation process. 
+```@example stDDM 
+σ = 1.0
+```
+
+### starting time
+The starting time parameter $s$ denotes how much earlier one attribute begins to affect the evidence accumulation process relative to the other(s). If $s$ is negative, attribute 1 evidence is accumulated before attribute 2 evidence; if $s$ is positive, attribute 1 evidence is accumulated after attribute 2 evidence. The absolute value of $s$ indicates the difference in starting times for the two attributes.
+```@example stDDM 
+s = 0.10 
+```
+
+### Starting Point
+An indicator of an an initial bias towards a decision. The z parameter is relative to a (i.e. it ranges from 0 to 1).
+```@example stDDM 
+z = 0.50
+```
+
+### Drift Rates Dispersion
+Dispersion parameters of the drift rate are drawn from a multivariate normal distribution, with the mean vector ν describing the distribution of actual drift rates from specific trials. The standard deviation or across-trial variability is captured by the η vector, and the corresponding correlation between the two attributes is denoted by ρ.
+```@example stDDM
+η = [1.0,1.0]
+ρ = 0.3
+```
+
+### Threshold
+The threshold α represents the amount of evidence required to make a decision.
+```@example stDDM 
+α = 1.5
+```
+
+### Non-Decision Time
+Non-decision time is an additive constant representing encoding and motor response time. 
+```@example stDDM 
+τ = 0.30
+```
+
+### stDDM Constructor 
+Now that values have been asigned to the parameters, we will pass them to `stDDM` to generate the model object.
+```@example stDDM 
+dist = stDDM(;ν, σ, s, z, η, ρ, α, τ,)
+```
+
+## Simulate Model
+Now that the model is defined, we will generate $10,000$ choices and reaction times using `rand`. 
+ ```@example stDDM 
+ choices,rts = rand(dist, 10_000)
+```
+
+## Compute Choice Probability
+The choice probability $\Pr(C=c)$ is computed by passing the model and choice index to `cdf`.
+ ```@example stDDM 
+cdf(dist, 1)
+```
+To compute the joint probability of choosing $c$ within $t$ seconds, i.e., $\Pr(T \leq t \wedge C=c)$, pass a third argument for $t$.
+
+## Plot Simulation
+The code below overlays the PDF on reaction time histograms for each option.
+ ```@example stDDM 
+histogram(dist)
+plot!(dist; t_range=range(.301, 3.0, length=100))
+```
+
+# References
+
+Amasino, D.R., Sullivan, N.J., Kranton, R.E. et al. Amount and time exert independent influences on intertemporal choice. Nat Hum Behav 3, 383–392 (2019). https://doi.org/10.1038/s41562-019-0537-2
+
+Barakchian, Z., Beharelle, A.R. & Hare, T.A. Healthy decisions in the cued-attribute food choice paradigm have high test-retest reliability. Sci Rep, (2021). https://doi.org/10.1038/s41598-021-91933-6
+
+Chen, HY., Lombardi, G., Li, SC. et al. Older adults process the probability of winning sooner but weigh it less during lottery decisions. Sci Rep, (2022). https://doi.org/10.1038/s41598-022-15432-y
+
+Maier, S.U., Raja Beharelle, A., Polanía, R. et al. Dissociable mechanisms govern when and how strongly reward attributes affect decisions. Nat Hum Behav 4, 949–963 (2020). https://doi.org/10.1038/s41562-020-0893-y
+
+Sullivan, N.J., Huettel, S.A. Healthful choices depend on the latency and rate of information accumulation. Nat Hum Behav 5, 1698–1706 (2021). https://doi.org/10.1038/s41562-021-01154-0
\ No newline at end of file
diff --git a/src/SequentialSamplingModels.jl b/src/SequentialSamplingModels.jl
index 1fe262bb..e109c7a8 100644
--- a/src/SequentialSamplingModels.jl
+++ b/src/SequentialSamplingModels.jl
@@ -28,6 +28,7 @@ module SequentialSamplingModels
     import Distributions: rand
     import Distributions: std
     import StatsAPI: params
+    import StatsBase:cor2cov
 
     export AbstractDDM
     export AbstractaDDM
@@ -36,7 +37,8 @@ module SequentialSamplingModels
     export AbstractLCA
     export AbstractLNR
     export AbstractPoissonRace
-    export AbstractRDM
+    export AbstractRDM 
+    export AbstractstDDM 
     export AbstractWald
     export aDDM
     export CDDM
@@ -50,6 +52,7 @@ module SequentialSamplingModels
     export PoissonRace
     export SSM1D
     export SSM2D
+    export stDDM
     export ContinuousMultivariateSSM
     export Wald
     export WaldMixture
@@ -92,4 +95,5 @@ module SequentialSamplingModels
     include("ext_functions.jl")
     include("ex_gaussian.jl")
     include("poisson_race.jl")
+    include("stDDM.jl")
 end
diff --git a/src/stDDM.jl b/src/stDDM.jl
new file mode 100644
index 00000000..9c596db5
--- /dev/null
+++ b/src/stDDM.jl
@@ -0,0 +1,201 @@
+"""
+    stDDM{T<:Real} <: AbstractstDDM
+
+An object for the starting-time diffusion decision model.
+
+# Parameters 
+
+- `ν`:  vector of drift rate weights for attribute one and two
+- `σ`:  diffusion noise
+- `s`:  initial latency bias (positive for attribute two, negative for attribute one)
+- `z`:  initial evidence 
+- `η`:  vector of variability in drift rate for attribute one and two
+- `ρ`:  correlation between drift rate for attributes
+- `α`:  evidence threshold
+- `τ`:  non-decision time 
+
+# Constructors 
+
+    stDDM(ν, σ, s, z, η, ρ, α, τ)
+
+    stDDM(;ν = [0.5,0.6],σ = 1,s = 0.50, z = 0.50, η = [1.0,1.0], ρ = 0.00, α = 1.0, τ = .300)
+        
+# Example 
+
+```julia 
+using SequentialSamplingModels
+
+ν = [0.5, 0.6]
+σ = 1
+s = 0.50
+z = 0.50
+η = [1.0, 1.0]
+ρ = 0.00
+α = 1.0
+τ = 0.300
+
+# Create stDDM model instance
+dist = stDDM(;ν, σ, s, z, η, ρ, α, τ)
+
+choices,rts = rand(dist, 500)
+```
+
+# References
+
+Amasino, D.R., Sullivan, N.J., Kranton, R.E. et al. Amount and time exert independent influences on intertemporal choice. Nat Hum Behav 3, 383–392 (2019). https://doi.org/10.1038/s41562-019-0537-2
+
+Barakchian, Z., Beharelle, A.R. & Hare, T.A. Healthy decisions in the cued-attribute food choice paradigm have high test-retest reliability. Sci Rep, (2021). https://doi.org/10.1038/s41598-021-91933-6
+
+Chen, HY., Lombardi, G., Li, SC. et al. Older adults process the probability of winning sooner but weigh it less during lottery decisions. Sci Rep, (2022). https://doi.org/10.1038/s41598-022-15432-y
+
+Lombardi, G., & Hare, T. Piecewise constant averaging methods allow for fast and accurate hierarchical Bayesian estimation of drift diffusion models with time-varying evidence accumulation rates. PsyArXiv, (2021). https://doi.org/10.31234/osf.io/5azyx
+
+Sullivan, N.J., Huettel, S.A. Healthful choices depend on the latency and rate of information accumulation. Nat Hum Behav 5, 1698–1706 (2021). https://doi.org/10.1038/s41562-021-01154-0
+"""
+
+mutable struct stDDM{T<:Real} <: AbstractstDDM
+    ν::Vector{T}
+    σ::T
+    s::T
+    z::T
+    η::Vector{T}
+    ρ::T
+    α::T
+    τ::T
+end
+
+function stDDM(ν, σ, s, z, η, ρ, α, τ)
+    _, σ ,s ,z ,_ ,ρ , α, τ = promote(ν[1],σ, s, z, η[1], ρ, α, τ)
+    ν = convert(Vector{typeof(τ)}, ν)
+    η = convert(Vector{typeof(τ)}, η)
+    return stDDM(ν, σ, s, z, η, ρ, α, τ)
+end
+
+function stDDM(;ν = [0.5,0.6], 
+    σ = 1,
+    s = 0.50, 
+    z = 0.50, 
+    η = fill(1.0, length(ν)), 
+    ρ = 0.0,
+    α = 1.0, 
+    τ = .300
+    )
+    return stDDM(ν, σ, s, z, η, ρ, α, τ)
+end
+
+function params(d::AbstractstDDM)
+    (d.ν, d.σ, d.s, d.z, d.η, d.ρ, d.α, d.τ)    
+end
+
+get_pdf_type(d::AbstractstDDM) = Approximate
+
+"""
+    rand(dist::AbstractstDDM)
+
+Generate a random choice-rt pair for starting-time diffusion decision model.
+
+# Arguments
+- `rng`: a random number generator
+- `dist`: model object for the starting-time diffusion decision model. 
+- `Δt`: time-step for simulation
+"""
+function rand(rng::AbstractRNG, dist::AbstractstDDM; kwargs...)
+    return simulate_trial(rng, dist; kwargs...)
+end
+
+"""
+    simulate_trial(rng::AbstractRNG, dist::AbstractstDDM;  Δt, max_steps)
+
+Generate a single simulated trial from the starting-time diffusion decision model.
+
+# Arguments
+
+- `rng`: a random number generator
+- `model::AbstractstDDM`: a starting-time diffusion decision model object
+- `Δt`: time-step for simulation
+- `max_steps`: total/max time for simulation
+"""
+function simulate_trial(rng::AbstractRNG, d::AbstractstDDM; Δt = .001, max_steps=6)
+    (;ν, σ, s, z, η, ρ, α, τ) = d
+
+    lt = Int(max_steps / Δt)
+    start_step = abs(Int(s / Δt))
+    
+    t = τ
+    choice = 0  # Initialize choice with a default value
+
+    X = z * α
+    deciding = true
+    cont = 1
+    Ρ = [1.0 ρ; ρ 1.0]    
+    Σ = cor2cov(Ρ,η)
+   
+    evidence = rand(MvNormal([ν[1], ν[2]], Σ))
+
+    while deciding && cont <= lt
+        δ1 = cont ≤ start_step && s > 0 ? 0.0 : 1.0  
+        δ2 = cont ≤ start_step && s < 0 ? 0.0 : 1.0 
+
+        noise = rand(rng, Normal(0, σ)) * √(Δt)
+
+        X += (evidence[1] * δ1 + evidence[2] * δ2) * Δt + noise
+        # X += (evidence[1] * a1Δ * δ1 + evidence[2] * a1Δ * δ2) * Δt + noise #note something to consider and fix is attribute difference
+        if X > α
+            choice = 1
+            deciding = false
+        elseif X < 0
+            choice = 2
+            deciding = false
+        end
+        t += Δt
+
+        cont += 1
+    end
+    
+    return (;choice,rt=t)
+
+end
+
+
+"""
+    simulate(rng::AbstractRNG, model::AbstractstDDM; Δt)
+
+Returns a matrix containing evidence samples of the stDDM decision process. In the matrix, rows 
+represent samples of evidence per time step and columns represent different accumulators.
+
+# Arguments
+
+- `rng`: a random number generator
+- `model::AbstractstDDM`: a starting-time diffusion decision model diffusion model object
+- `Δt`: time-step for simulation
+"""
+function simulate(rng::AbstractRNG, model::AbstractstDDM; Δt = .001)
+    (;ν, σ, s, z, η, ρ, α, τ) = model
+
+    x = α * z
+    t = 0.0
+    evidence = [x]
+    time_steps = [t]
+    cont = 1
+    start_step = abs(Int(s / Δt))
+
+    Ρ = [1.0 ρ; ρ 1.0]    
+    Σ = cor2cov(Ρ,η) 
+    while (x < α) && (x > 0)
+        t += Δt
+
+        δ1 = cont ≤ start_step && s > 0 ? 0.0 : 1.0  
+        δ2 = cont ≤ start_step && s < 0 ? 0.0 : 1.0 
+
+        increment = rand(rng, MvNormal([ν[1], ν[2]], Σ))
+        noise = rand(rng, Normal(0, σ)) * √(Δt)
+
+        x += (increment[1] * δ1 + increment[2] * δ2) * Δt + noise
+
+        push!(evidence, x)
+        push!(time_steps, t)
+        cont += 1
+    end
+
+    return time_steps,evidence
+end
\ No newline at end of file
diff --git a/src/type_system.jl b/src/type_system.jl
index 1c61846d..553a237c 100644
--- a/src/type_system.jl
+++ b/src/type_system.jl
@@ -70,7 +70,14 @@ abstract type AbstractLCA <: SSM2D end
 
 An abstract type for the Poisson race model.
 """
-abstract type AbstractPoissonRace <: SSM2D end
+abstract type AbstractPoissonRace <:SSM2D end
+
+"""
+    AbstractstDDM <: SSM2D
+
+An abstract type for the starting-time diffusion decision model.
+"""
+abstract type AbstractstDDM <: SSM2D end
 
 """
     AbstractRDM <: SSM2D
@@ -116,7 +123,7 @@ Base.length(d::SSM2D) = 2
 rand(d::SSM2D; kwargs...) = rand(Random.default_rng(), d; kwargs...)
 rand(d::ContinuousMultivariateSSM; kwargs...) = rand(Random.default_rng(), d; kwargs...)
 rand(d::ContinuousMultivariateSSM, n::Int; kwargs...) =
-    rand(Random.default_rng(), d, n; kwargs...)
+rand(Random.default_rng(), d, n; kwargs...)
 
 """
     rand(rng::AbstractRNG, d::SSM2D, N::Int; kwargs...)
@@ -127,20 +134,21 @@ with more than one choice option.
 # Arguments
 
 - `d::SSM2D`: a 2D sequential sampling model.
-- `n_sim::Int`: the number of simulated choices and rts  
+- `N::Int`: the number of simulated choices and rts  
 
 # Keywords
 
 - `kwargs...`: optional keyword arguments 
 """
-function rand(rng::AbstractRNG, d::SSM2D, n_sim::Int)
-    choice = fill(0, n_sim)
-    rt = fill(0.0, n_sim)
-    for i = 1:n_sim
-        choice[i], rt[i] = rand(rng, d)
+function rand(rng::AbstractRNG, d::SSM2D, N::Int; kwargs...)
+    choice = fill(0, N)
+    rt = fill(0.0, N)
+    for i = 1:N
+        choice[i], rt[i] = rand(rng, d; kwargs...)
     end
     return (; choice, rt)
 end
+rand(d::SSM2D, N::Int; kwargs...) = rand(Random.default_rng(), d, N; kwargs...)
 
 """
     logpdf(d::SSM2D, data::NamedTuple) 
@@ -183,8 +191,7 @@ Computes the probability density for a 2D sequential sampling model.
 pdf(d::SSM2D, data::NamedTuple, args...; kwargs...) =
     pdf.(d, data.choice, data.rt, args...; kwargs...)
 
-pdf(d::SSM2D, data::AbstractArray{Real,2}, args...; kwargs...) =
-    pdf(d, Int(data[1]), data[2], args...; kwargs...)
+pdf(d::SSM2D, data::AbstractArray{Real,2}) = pdf(d, Int(data[1]), data[2])
 
 """
     cdf(d::SSM2D, choice::Int, ub=10)
diff --git a/test/stddm_tests.jl b/test/stddm_tests.jl
new file mode 100644
index 00000000..8e7ad2fb
--- /dev/null
+++ b/test/stddm_tests.jl
@@ -0,0 +1,192 @@
+@safetestset "stDDM Tests" begin
+    @safetestset "stDDM predictions" begin
+        using SequentialSamplingModels
+        using Test
+        using Random
+        Random.seed!(8414)
+
+        parms = (ν = [1.0, 1.2], α = 0.8, τ = 0.30 , s = 0.0, z = 0.50)
+        model = stDDM(; parms...)
+        choice, rt = rand(model, 10_000)
+        # use validated simulator
+        #Lombardi, G., & Hare, T. Piecewise constant averaging methods allow for fast and 
+        # accurate hierarchical Bayesian estimation of drift diffusion models with 
+        # time-varying evidence accumulation rates. PsyArXiv, (2021). https://doi.org/10.31234/osf.io/5azyx
+        #code to validate
+        #https://github.com/galombardi/method_HtSSM_aDDM/tree/master/RecoveryFitting/stDDM
+
+        @test mean(choice .== 1) ≈  0.8192  atol = 1e-1
+        @test mean(rt[choice.==1]) ≈ 0.4332 atol = 1e-1
+        @test mean(rt[choice.==2]) ≈ 0.4555 atol = 1e-1
+        @test std(rt[choice.==1]) ≈ 0.1106 atol = 1e-1
+        @test std(rt[choice.==2]) ≈ 0.1318 atol = 1e-1
+
+        parms = (ν = [2, 1], α =1.5, τ = 0.30 , s = 0.20, z = 0.50)
+        model = stDDM(; parms...)
+        choice, rt = rand(model, 10_000)
+
+        #R:      0.8883 0.6639 0.6303 0.2377 0.3146
+        #Python: 0.8918 0.6694 0.6285 0.2478 0.3156
+        #Julia:  0.8897 0.6635 0.6298 0.2464 0.3037
+
+        @test mean(choice .== 1)   ≈  0.8883 atol = 1e-1
+        @test mean(rt[choice.==1]) ≈  0.6639  atol = 1e-1
+        @test mean(rt[choice.==2]) ≈  0.6303  atol = 1e-1
+        @test std(rt[choice.==1])  ≈  0.2377 atol = 1e-1
+        @test std(rt[choice.==2])  ≈  0.3146 atol = 1e-1
+
+        parms = (ν = [.7, .9], α =2, τ = 0.40 , s = 0.60, z = 0.40)
+        model = stDDM(; parms...)
+        choice, rt = rand(model, 10_000)
+
+        #R:      0.7417 1.1632 1.0788 0.5329 0.6273
+        #Python: 0.7345 1.1686 1.0444 0.5515 0.6125
+        #Julia:  0.7335 1.1642 1.0546 0.533 0.6285
+
+        @test mean(choice .== 1)   ≈  0.7417 atol = 1e-1
+        @test mean(rt[choice.==1]) ≈  1.1632  atol = 1e-1
+        @test mean(rt[choice.==2]) ≈  1.0788  atol = 1e-1
+        @test std(rt[choice.==1])  ≈  0.5329 atol = 1e-1
+        @test std(rt[choice.==2])  ≈  0.6273 atol = 1e-1
+
+    end
+    
+#   R:
+#   Rstddm <- function(d_v, d_h, thres, nDT, tIn, bias, sd_n, N, seed = NULL) {
+#   if (!is.null(seed)) set.seed(seed)
+  
+#   T <- 6 # Total time
+#   dt <- 0.001 # Time step
+#   lt <- as.integer(T / dt) # Number of time steps
+#   vec_tHealth <- rep(1, lt)
+#   vec_tVal <- rep(1, lt)
+#   aux <- abs(as.integer(tIn / dt))
+  
+#   # Adjusting health and value vectors based on tIn
+#   if (tIn > 0) {
+#     vec_tVal[1:aux] <- 0
+#   } else if (tIn < 0) {
+#     vec_tHealth[1:aux] <- 0
+#   }
+  
+#   simulate_ddm <- function() {
+#     X <- bias * thres # Initial accumulation value
+#     flag <- FALSE
+#     cont <- 0
+#     vecOut <- T
+    
+#     Sigma = matrix(c(1,0,0,1),2,2)
+#     d_vht <- MASS::mvrnorm(1, c(d_v, d_h), Sigma)
+#     while (!flag && cont < lt) {
+#       noise <- rnorm(1, mean = 0, sd = sd_n) * sqrt(dt)
+#       X <- X + (d_vht[[1]] *1* vec_tVal[cont + 1] + d_vht[[2]] *1* vec_tHealth[cont + 1]) * dt + noise #note 1's are VD and HD
+      
+#       if (X > thres) {
+#         flag <- TRUE
+#         vecOut <- nDT + cont * dt
+#       } else if (X < 0) {
+#         flag <- TRUE
+#         vecOut <- -nDT - cont * dt
+#       }
+#       cont <- cont + 1
+#     }
+    
+#     return(vecOut)
+#   }
+  
+#   results <- numeric(N)
+#   for (i in 1:N) {
+#     results[i] <- simulate_ddm()
+#   }
+  
+#   return(results)
+# }
+#res = Rstddm(d_v = 1, d_h = 1.2, thres = 0.8, nDT = 0.30, tIn = 0.0, bias = .5, sd_n = 1, N = 10000, seed = 8414)
+#res = Rstddm(d_v = 2, d_h = 1, thres = 1.5, nDT = 0.30, tIn = 0.2, bias = .5, sd_n = 1, N = 10000, seed = 8414)
+#res = Rstddm(d_v = .7, d_h = .9, thres = 2, nDT = 0.40, tIn = 0.6, bias = .5, sd_n = 1, N = 10000, seed = 8414)
+#
+#   Python:
+# import numpy as np
+# def pystddm(d_v, d_h, thres, nDT, tIn, bias, sd_n, N, seed=None):
+#     # Set a seed for reproducibility if provided
+#     if seed is not None:
+#         np.random.seed(seed)
+    
+#     # Constants and Preparations
+#     T = 6  # Total time
+#     dt = 0.001  # Time step
+#     lt = int(T / dt)  # Number of time steps
+#     vec_tHealth = np.ones(lt)
+#     vec_tVal = np.ones(lt)
+#     aux = abs(int(tIn / dt))
+    
+#     # Adjusting health and value vectors based on tIn
+#     if tIn > 0:
+#         vec_tVal[:aux] = 0
+#     elif tIn < 0:
+#         vec_tHealth[:aux] = 0
+    
+#     def simulate_ddm():
+#         X = bias * thres  # Initial accumulation value
+#         flag = False
+#         cont = 0
+#         vecOut = T
+        
+#         Sigma = np.array([[1, 0], [0, 1]])
+#         d_vht = np.random.multivariate_normal([d_v, d_h], Sigma)
+        
+#         while not flag and cont < lt:
+#             noise = np.random.normal(0, sd_n) * np.sqrt(dt)
+#             X += (d_vht[0] * vec_tVal[cont] + d_vht[1] * vec_tHealth[cont]) * dt + noise
+            
+#             if X > thres:
+#                 flag = True
+#                 vecOut = nDT + cont * dt
+#             elif X < 0:
+#                 flag = True
+#                 vecOut = -nDT - cont * dt
+#             cont += 1
+        
+#         return vecOut
+    
+#     # Apply the simulation for each output element
+#     results = [simulate_ddm() for _ in range(N)]
+    
+#     return results
+
+#res = pystddm(d_v = 1, d_h = 1.2, thres = 0.8, nDT = 0.30, tIn = 0.0, bias = .5, sd_n = 1, N = 10000, seed = 8414)
+#res = pystddm(d_v = 2, d_h = 1, thres = 1.5, nDT = 0.30, tIn = 0.2, bias = .5, sd_n = 1, N = 10000, seed = 8414)
+#res = pystddm(
+#     d_v=.7,    # Weight of the drift for the first attribute
+#     d_h=.9,    # Weight of the drift for the second attribute
+#     thres=2,  # Decision threshold
+#     nDT=0.4,  # Additional decision time
+#     tIn=0.60,  # Time influence of the second attribute
+#     bias=.4,   # Initial value (starting point) of the accumulation
+#     sd_n=1,   # Standard deviation of the noise
+#     N=10000,    # Number of simulations
+#     seed=8414  # Seed for reproducibility
+# )
+
+    @safetestset "simulate" begin
+        using SequentialSamplingModels
+        using Test
+        using Random
+
+        Random.seed!(8477)
+        α = 0.80
+        Δt = 0.0001
+        dist = stDDM(; α, ν = [3, 3])
+
+        time_steps, evidence = simulate(dist; Δt)
+
+        @test time_steps[1] ≈ 0
+        @test length(time_steps) == length(evidence)
+        @test evidence[end] ≈ α atol = 0.05
+
+        dist = stDDM(; α, ν = [-3, -3])
+        time_steps, evidence = simulate(dist; Δt)
+        @test evidence[end] ≈ 0.0 atol = 0.05
+    end
+
+end