diff --git a/src/DiffFusion.jl b/src/DiffFusion.jl
index 1c8d837c..4f8887aa 100644
--- a/src/DiffFusion.jl
+++ b/src/DiffFusion.jl
@@ -17,6 +17,7 @@ using Random
 using Roots
 using Sobol
 using SparseArrays
+using StatsBase
 using Zygote
 
 import Base.length
@@ -74,6 +75,7 @@ include("payoffs/RatesPayoffs.jl")
 include("payoffs/RatesOptions.jl")
 include("payoffs/AmcPayoffs.jl")
 include("payoffs/AssetOptions.jl")
+include("payoffs/BarrierOptions.jl")
 
 include("products/Cashflows.jl")
 include("products/AssetOptionFlows.jl")
@@ -88,7 +90,9 @@ include("products/CashAndAssetLegs.jl")
 include("products/BermudanSwaptionLeg.jl")
 
 include("utils/Bachelier.jl")
+include("utils/Barriers.jl")
 include("utils/Black.jl")
+include("utils/BrownianBridge.jl")
 include("utils/Gradients.jl")
 include("utils/Integrations.jl")
 include("utils/InterpolationMethods.jl")
diff --git a/src/paths/AbstractPath.jl b/src/paths/AbstractPath.jl
index b284eef8..b6dc9a1a 100644
--- a/src/paths/AbstractPath.jl
+++ b/src/paths/AbstractPath.jl
@@ -94,6 +94,15 @@ function forward_asset(p::AbstractPath, t::ModelTime, T::ModelTime, key::String)
     error("AbstractPath needs to implement forward_asset method.")
 end
 
+"""
+    forward_asset_and_zero_bonds(p::AbstractPath, t::ModelTime, T::ModelTime, key::String)
+
+Calculate asset and zero bond components for forward asset price calculation.
+"""
+function forward_asset_and_zero_bonds(p::AbstractPath, t::ModelTime, T::ModelTime, key::String)
+    error("AbstractPath needs to implement forward_asset_and_zero_bonds method.")
+end
+
 
 """
     fixing(p::AbstractPath, t::ModelTime, key::String)
diff --git a/src/paths/PathMethods.jl b/src/paths/PathMethods.jl
index 1dbdb7e9..3cd08863 100644
--- a/src/paths/PathMethods.jl
+++ b/src/paths/PathMethods.jl
@@ -350,6 +350,67 @@ function forward_asset(p::Path, t::ModelTime, T::ModelTime, key::String)
 end
 
 
+
+"""
+    forward_asset_zero_bonds(p::Path, t::ModelTime, T::ModelTime, key::String)
+
+Calculate asset (plus deterministic jumps) as well as domestic and foreign
+zero bond price associated with an `Asset` key.
+
+This function implements methodology redundant to `forward_asset(...)`. But it
+returns asset and zero bonds separately.
+
+This function is used for barrier option pricing.
+"""
+function forward_asset_and_zero_bonds(p::Path, t::ModelTime, T::ModelTime, key::String)
+    @assert t <= T
+    (context_key, ts_key_for, ts_key_dom, op) = context_keys(key)
+    @assert op in (_empty_context_key, "-")
+    entry = p.context.assets[context_key]
+    #
+    spot_alias = entry.asset_spot_alias
+    spot = p.ts_dict[spot_alias](T, TermstructureScalar)  # capture any discrete jumps until T to be consistent with forward_asset
+    #
+    ts_alias_dom = entry.domestic_termstructure_dict[ts_key_dom]
+    df_dom_t = discount(t, p.ts_dict, ts_alias_dom)
+    df_dom_T = discount(T, p.ts_dict, ts_alias_dom)
+    #
+    ts_alias_for = entry.foreign_termstructure_dict[ts_key_for]
+    df_for_t = discount(t, p.ts_dict, ts_alias_for)
+    df_for_T = discount(T, p.ts_dict, ts_alias_for)
+    #
+    if isnothing(entry.asset_model_alias) &&
+        isnothing(entry.foreign_model_alias) &&
+        isnothing(entry.domestic_model_alias)
+        # we can take a short-cut for fully deterministic models
+        e = ones(length(p))
+        return ((spot*df_for_t/df_dom_t) .*e, (df_dom_T/df_dom_t).*e, (df_for_T/df_for_t).*e)
+    end
+    #
+    X = state_variable(p.sim, t, p.interpolation)
+    SX = model_state(X, p.state_alias_dict)
+    #
+    y_ast = zeros(length(p))
+    y_dom = zeros(length(p))
+    y_for = zeros(length(p))
+    if !isnothing(entry.asset_model_alias)
+        y_ast .+= log_asset(p.sim.model, entry.asset_model_alias, t, SX)
+    end
+    if !isnothing(entry.domestic_model_alias)
+        y_ast .+= log_bank_account(p.sim.model, entry.domestic_model_alias, t, SX)
+        y_dom .-= log_zero_bond(p.sim.model, entry.domestic_model_alias, t, T, SX)
+    end
+    if !isnothing(entry.foreign_model_alias)
+        y_ast .-= log_bank_account(p.sim.model, entry.foreign_model_alias, t, SX)
+        y_for .-= log_zero_bond(p.sim.model, entry.foreign_model_alias, t, T, SX)
+    end
+    asset = (spot * df_for_t / df_dom_t) .* exp.(y_ast)
+    zb_dom = (df_dom_T / df_dom_t) .* exp.(y_dom)
+    zb_for = (df_for_T / df_for_t) .* exp.(y_for)
+    return (asset, zb_dom, zb_for)
+end
+
+
 """
     fixing(p::Path, t::ModelTime, key::String)
 
diff --git a/src/payoffs/AssetOptions.jl b/src/payoffs/AssetOptions.jl
index c472d31f..fb2b48b4 100644
--- a/src/payoffs/AssetOptions.jl
+++ b/src/payoffs/AssetOptions.jl
@@ -13,6 +13,9 @@ at `expiry_time`.
 
 Option forward price is calculated as expectation in T-forward measure where T corresponds
 to the expiry time. Conditioning (for time-t price) is on information at `obs_time`.
+
+Strike price `strike_price` must be time-t (`obs_time`) measurable. Otherwise, we *look into
+the future*.
 """
 struct VanillaAssetOption <: Payoff
     obs_time::ModelTime
diff --git a/src/payoffs/BarrierOptions.jl b/src/payoffs/BarrierOptions.jl
new file mode 100644
index 00000000..ad64b3da
--- /dev/null
+++ b/src/payoffs/BarrierOptions.jl
@@ -0,0 +1,231 @@
+
+"""
+    struct BarrierAssetOption <: Payoff
+        obs_time::ModelTime
+        expiry_time::ModelTime
+        forward_price::ForwardAsset
+        strike_price::Payoff
+        call_put::ModelValue
+        barrier_level::Payoff
+        barrier_direction::ModelValue
+        barrier_type::ModelValue
+        rebate_price::ModelValue
+        no_hit_times::AbstractVector
+    end
+
+The time-t forward price of an option paying [ϕ(F-K)]^+. Forward asset price *F* is determined
+at `expiry_time`.
+
+Option forward price is calculated as expectation in T-forward measure where T corresponds
+to the expiry time. Conditioning (for time-t price) is on information at `obs_time`. This
+requires particular care when using Black-Scholes pricing functions.
+
+Strike price `strike_price` and barrier level `barrier_level` must be time-t (`obs_time`)
+measurable. Otherwise, we *look into the future*.
+
+`barrier_direction` is -1 for up-barrier and +1 for down-barrier. `barrier_type` is -1 for
+in-barrier and +1 for out-barrier. See also `black_scholes_barrier_price`.
+
+`no_hit_times` is a list of times where past hit events are observed and with which
+no-hit probability is estimated. First time is zero and last time is `obs_time`. Must be
+of length 2 or more.
+"""
+struct BarrierAssetOption <: Payoff
+    obs_time::ModelTime
+    expiry_time::ModelTime
+    forward_price::ForwardAsset
+    strike_price::Payoff
+    call_put::ModelValue
+    barrier_level::Payoff
+    barrier_direction::ModelValue
+    barrier_type::ModelValue
+    rebate_price::ModelValue
+    no_hit_times::AbstractVector
+end
+
+
+"""
+    BarrierAssetOption(
+        forward_price::ForwardAsset,
+        strike_price::Payoff,
+        barrier_level::Payoff,
+        option_type::String,
+        rebate_price::ModelValue,
+        number_of_no_hit_times::Integer,
+        )
+
+Create a `BarrierAssetOption` payoff.
+
+String `option_type` is of the form [D|U][O|I][C|P].
+
+No-hit times are calculates as equally spaced times of
+length `number_of_no_hit_times`.
+"""
+function BarrierAssetOption(
+    forward_price::ForwardAsset,
+    strike_price::Payoff,
+    barrier_level::Payoff,
+    option_type::String,
+    rebate_price::ModelValue,
+    number_of_no_hit_times::Integer,
+    )
+    #
+    @assert obs_time(strike_price) ≤ forward_price.obs_time  # payoff must be time-t measurable
+    @assert length(option_type) == 3
+    @assert number_of_no_hit_times ≥ 2
+    # derive option properties
+    o_type = uppercase(option_type)
+    @assert o_type[1] in ('D', 'U')
+    @assert o_type[2] in ('O', 'I')
+    @assert o_type[3] in ('C', 'P')
+    if o_type[1] == 'D'
+        η = 1
+    else
+        η = -1
+    end
+    if o_type[2] == 'O'
+        χ = 1
+    else
+        χ = -1
+    end
+    if o_type[3] == 'C'
+        ϕ = 1
+    else
+        ϕ = -1
+    end
+    #
+    no_hit_times = LinRange(0.0, forward_price.obs_time, number_of_no_hit_times)
+    #
+    return BarrierAssetOption(
+        forward_price.obs_time,
+        forward_price.maturity_time,
+        forward_price,
+        strike_price,
+        ϕ,
+        barrier_level,
+        η,
+        χ,
+        rebate_price,
+        no_hit_times,
+    )
+end
+
+
+"""
+    obs_time(p::BarrierAssetOption)
+
+Return BarrierAssetOption observation time.
+"""
+function obs_time(p::BarrierAssetOption)
+    return p.obs_time
+end
+
+
+"""
+    obs_times(p::BarrierAssetOption)
+
+Return all BarrierAssetOption observation times. 
+"""
+function obs_times(p::BarrierAssetOption)
+    times = Set(obs_time(p))
+    times = union(times, obs_times(p.strike_price))
+    times = union(times, p.no_hit_times)
+    return times
+end
+
+
+"""
+    at(p::BarrierAssetOption, path::AbstractPath)
+
+Evaluate a `BarrierAssetOption` at a given `path`, *X(omega)*.
+"""
+function at(p::BarrierAssetOption, path::AbstractPath)
+    # F = at(p.forward_price, path)
+    (S, Pd, Pf) = forward_asset_and_zero_bonds(
+        path,
+        p.forward_price.obs_time,
+        p.forward_price.maturity_time,
+        p.forward_price.key
+    )
+    F = S .* Pf ./ Pd
+    #
+    K = p.strike_price(path)
+    H = p.barrier_level(path)
+    ν² = asset_variance(path, p.obs_time, p.expiry_time, p.forward_price.key)
+    T = p.expiry_time - p.obs_time
+    σ = 0.0
+    if T > 0.0
+        σ = sqrt.(ν² / T)
+    end
+    # calculate hit-price
+    if p.barrier_type == 1  # out-Barrier
+        hit_price = p.rebate_price
+    else
+        # Vanilla option price
+        if all(ν² .≤ 0.0)
+            # calculate intrinsic value
+            hit_price = max.(p.call_put*(F - K), 0.0)
+        else
+            hit_price = black_price(K, F, sqrt.(ν²), p.call_put)
+        end
+    end
+    # calculate no-hit price
+    #
+    σ = max.(σ, 0.001)  # we need a positive volatility
+    T = max(T, 1.0/365/24) # we approximate intrinsic value by option value 1h prior to exercise
+
+    no_hit_price = black_scholes_barrier_price(
+        K,
+        H,
+        p.rebate_price,
+        p.barrier_direction,
+        p.barrier_type,
+        p.call_put,
+        S, Pd, Pd ./ Pf, σ, T
+    ) ./ Pd  # we want the forward price here
+    # calculate no hit on path
+    #
+    S = hcat(
+        [ asset(path, t, p.forward_price.key) for t in p.no_hit_times]...
+    )
+    logS = log.(S)
+    logS_returns = logS[:,2:end] .- logS[:,1:end-1]
+    variances = std(logS_returns, dims=1).^2
+    no_hit_prob = barrier_no_hit_probability(
+        log.(H),
+        p.barrier_direction,
+        logS,
+        variances,
+    )
+    #
+    return no_hit_prob .* no_hit_price .+ (1.0 .- no_hit_prob) .* hit_price
+end
+
+
+"""
+    string(p::VanillaAssetOption)
+
+Formatted (and shortened) output for VanillaAssetOption payoff.
+"""
+string(p::BarrierAssetOption) = begin
+    option_type = ""
+    if p.barrier_direction == -1
+        option_type = "U"
+    end
+    if p.barrier_direction == +1
+        option_type = "D"
+    end
+    if p.barrier_type == -1
+        option_type = option_type * "I"
+    end
+    if p.barrier_type == +1
+        option_type = option_type * "O"
+    end
+    if p.call_put == 1.0
+        option_type = option_type * "Call"
+    end
+    if p.call_put == -1.0
+        option_type = option_type * "Put"
+    end
+    @sprintf("%s(%s, X = %s, H = %s)", option_type, string(p.forward_price), string(p.strike_price), string(p.barrier_level))
+end
diff --git a/src/products/AssetOptionFlows.jl b/src/products/AssetOptionFlows.jl
index 50cd80dc..ddc3eaf9 100644
--- a/src/products/AssetOptionFlows.jl
+++ b/src/products/AssetOptionFlows.jl
@@ -35,7 +35,7 @@ end
 
 Return the payoff representing the simulated expected amount of the `VanillaAssetOptionFlow`.
 
-This implementation is an approximation and does not capture convexity adjustments.
+This implementation is an approximation and does not capture payment delay convexity adjustments.
 """
 function expected_amount(cf::VanillaAssetOptionFlow, obs_time::ModelTime)
     if obs_time ≥ cf.expiry_time
@@ -46,3 +46,68 @@ function expected_amount(cf::VanillaAssetOptionFlow, obs_time::ModelTime)
     return VanillaAssetOption(F, K, cf.call_put)
 end
 
+
+"""
+    struct BarrierAssetOptionFlow <: CashFlow
+        expiry_time::ModelTime
+        pay_time::ModelTime
+        strike_price::ModelValue
+        barrier_level::ModelValue
+        rebate_price::ModelValue
+        option_type::String
+        asset_key::String
+        hit_obs_step_size::ModelTime
+    end
+
+A `CashFlow` representing a Single Barrier option on an `Asset`.
+
+`option_type` is of the form [D|U][O|I][C|P].
+
+`hit_obs_step_size` represents a modelling parameter to build
+the grid of *past* observation times to monitor barrier hit
+along the path.
+"""
+struct BarrierAssetOptionFlow <: CashFlow
+    expiry_time::ModelTime
+    pay_time::ModelTime
+    strike_price::ModelValue
+    barrier_level::ModelValue
+    rebate_price::ModelValue
+    option_type::String
+    asset_key::String
+    hit_obs_step_size::ModelTime
+end
+
+
+"""
+    amount(cf::BarrierAssetOptionFlow)
+
+Return the payoff of the `BarrierAssetOptionFlow`.
+"""
+function amount(cf::BarrierAssetOptionFlow)
+    F = ForwardAsset(cf.expiry_time, cf.expiry_time, cf.asset_key)
+    X = Fixed(cf.strike_price)
+    H = Fixed(cf.barrier_level)
+    nnht = max(2, Int(round(cf.expiry_time / cf.hit_obs_step_size)))
+    return BarrierAssetOption(F, X, H, cf.option_type, cf.rebate_price, nnht)
+end
+
+
+"""
+    expected_amount(cf::BarrierAssetOptionFlow, obs_time::ModelTime)
+
+Return the payoff representing the simulated expected amount of the `BarrierAssetOptionFlow`.
+
+This implementation is an approximation and does not capture payment delay convexity adjustments.
+"""
+function expected_amount(cf::BarrierAssetOptionFlow, obs_time::ModelTime)
+    if obs_time ≥ cf.expiry_time
+        return amount(cf)
+    end
+    F = ForwardAsset(obs_time, cf.expiry_time, cf.asset_key)
+    X = Fixed(cf.strike_price)
+    H = Fixed(cf.barrier_level)
+    nnht = max(2, Int(round(obs_time / cf.hit_obs_step_size)))
+    return BarrierAssetOption(F, X, H, cf.option_type, cf.rebate_price, nnht)
+end
+
diff --git a/src/utils/Barriers.jl b/src/utils/Barriers.jl
new file mode 100644
index 00000000..1f82ac74
--- /dev/null
+++ b/src/utils/Barriers.jl
@@ -0,0 +1,201 @@
+
+"""
+    black_scholes_vanilla_price(X, ϕ, S, DF_r, DF_b, ν)
+
+Vanilla option pricing formula for calls and puts.
+
+Notation:
+  - X option strike
+  - ϕ = 1 (call), -1 (put); option type
+  - S underlying spot
+  - DF_r = exp(-r T), r is (domestic) risk-free rate
+  - DF_b = exp(-b T), b is *cost-of-carry*, i.e. b = r - q (see Haug, sec. 1.2.1)
+  - ν = σ √T, σ option volatility, T time to option expiry.
+"""
+function black_scholes_vanilla_price(X, ϕ, S, DF_r, DF_b, ν)
+    N(x) = cdf.(Normal(), x)
+    F = S ./ DF_b
+    d1 = log.(F ./ X) ./ ν .+ ν ./ 2
+    d2 = d1 .- ν
+    return DF_r .* ϕ .* (F .* N(ϕ .* d1) .- X .* N(ϕ .* d2))
+end
+
+
+"""
+    black_scholes_barrier_price(X, H, K, η, χ, ϕ, S, DF_r, DF_b, σ, T)
+
+Barrier option pricing formulas following Haug, "The Complete Guide to Option Pricing Formulas".
+
+Notation:
+  - r is (domestic) risk-free rate
+  - b is *cost-of-carry*, see Haug, sec. 1.2.1:
+      - b = r - q, where q is a continuous dividend yield,
+      - b = r - rf, where rf is a foreign currency rate,
+      - b = 0 for options on futures.
+  - DF_r = exp(-r T)
+  - DF_b = exp(-b T)
+
+  - σ option volatility
+  - T time to option expiry
+  - ν = σ √T
+
+  - S underlying spot
+  - X option strike
+  - H option barrier
+  - K option rebate
+
+  - η = 1 (down), -1 (up); barrier direction
+  - χ = 1 (out), -1 (in); barrier type
+  - ϕ = 1 (call), -1 (put); option type
+
+"""
+function black_scholes_barrier_price(X, H, K, η, χ, ϕ, S, DF_r, DF_b, σ, T)
+
+    ν = σ * sqrt(T)
+    if any(ν .<= 0.0)
+        error("Intrinsic value not yet implemented.")
+    end
+
+    down = Int(η == 1)
+    up = 1 - down
+    out = Int(χ == 1)
+    in = 1 - out
+    
+    barrier_hit = down .* (S .< H) + up .* (S .> H)  # OR of disjoint events
+    barrier_hit_price = in .* black_scholes_vanilla_price(X, ϕ, S, DF_r, DF_b, ν) # + out * DF_r * K, for rebate at expiry
+    
+    N(x) = cdf.(Normal(), x)
+    # scalar auxiliary variables
+    r = -log.(DF_r) ./ T
+    b = -log.(DF_b) ./ T
+    μ = (b .- (σ.^2)/2)./(σ.^2)
+    λ = sqrt.(max.(μ.^2 .+ 2 .* r ./ σ.^2, 0.0))
+    α = (1 .+ μ) .* ν
+
+    # common vector-valued input terms
+    log_S_X = log.(S ./ X)
+    log_S_H = log.(S ./ H)
+    log_H_X = log.(H ./ X)
+    x1 = log_S_X ./ ν .+ α                # -> A
+    x2 = log_S_H ./ ν .+ α                # -> B, E
+    y1 = (log_H_X .- log_S_H) ./ ν .+ α   # -> C
+    y2 = -log_S_H ./ ν .+ α               # -> D, E
+    z  = -log_S_H ./ ν .+ (λ .* ν)        # -> F
+
+    # option component terms; TODO: move to where it is really needed and avoid unnecessary calculations
+    A = ϕ .* S .* DF_r./DF_b                       .* N(ϕ .* x1) .- ϕ .* X .* DF_r                  .* N(ϕ .* x1 .- ϕ .* ν)
+    B = ϕ .* S .* DF_r./DF_b                       .* N(ϕ .* x2) .- ϕ .* X .* DF_r                  .* N(ϕ .* x2 .- ϕ .* ν)
+    C = ϕ .* S .* DF_r./DF_b .* (H./S).^(2*(μ.+1)) .* N(η .* y1) .- ϕ .* X .* DF_r .* (H./S).^(2*μ) .* N(η .* y1 .- η .* ν)
+    D = ϕ .* S .* DF_r./DF_b .* (H./S).^(2*(μ.+1)) .* N(η .* y2) .- ϕ .* X .* DF_r .* (H./S).^(2*μ) .* N(η .* y2 .- η .* ν)
+    #
+    E = K .* DF_r .* (                  N(η .* x2 .- η .* ν) .- (H./S).^(2 .* μ) .* N(η .* y2           .- η .* ν))
+    F = K         .* ((H./S).^(μ.+λ) .* N(η .* z)            .+ (H./S).^(μ .- λ) .* N(η .* z  .- 2 .* λ .* η .* ν))
+
+    # if KO rebate is paid at expiry: λ -> μ
+    z_T_e = -log_S_H ./ ν .+ (μ .* ν)  # -> F_T_e
+    F_T_e = K .* ((H./S).^(2*μ)     .* N(η * z)              .+                     N(η .* z  .- 2 .* μ .* η .* ν))
+    
+    # For the no-hit case we need to distinguish 16 different cases.
+    # We follow the ordering in Haug
+    if χ == -1  # in-options
+        if ϕ == 1  # call
+            if η == 1  # DIC
+                if X > H
+                    no_hit_price = C
+                else
+                    no_hit_price = A .- B .+ D
+                end
+            else  # UIC
+                if X > H
+                    no_hit_price = A
+                else
+                    no_hit_price = B .- C .+ D
+                end
+            end
+        else  # put
+            if η == 1  # DIP
+                if X > H
+                    no_hit_price = B .- C .+ D
+                else
+                    no_hit_price = A
+                end
+            else  # UIP
+                if X > H
+                    no_hit_price = A .- B .+ D
+                else
+                    no_hit_price = C
+                end
+            end
+        end
+        # add rebate for in-options (if needed)
+        no_hit_price = no_hit_price .+ E
+    else  # out-options
+        if ϕ == 1  # call
+            if η == 1  # DOC
+                if X > H
+                    no_hit_price = A .- C
+                else
+                    no_hit_price = B .- D
+                end
+            else  # UOC
+                if X > H
+                    no_hit_price = 0 .* S  # beware dimensions
+                else
+                    no_hit_price = A - B + C - D
+                end
+            end
+        else  # put
+            if η == 1  # DOP
+                if X > H
+                    no_hit_price = A .- B .+ C .- D
+                else
+                    no_hit_price = 0 .* S  # beware dimensions
+                end
+            else  # UOP
+                if X > H
+                    no_hit_price = B .- D
+                else
+                    no_hit_price = A .- C
+                end
+            end
+        end
+        # add rebate (at KO-time) for out-options (if needed)
+        no_hit_price = no_hit_price .+ F
+    end
+    
+    return barrier_hit .* barrier_hit_price + (1.0 .- barrier_hit) .* no_hit_price
+end
+
+"""
+    black_scholes_barrier_price(X, H, K, option_type::String, S, DF_r, DF_b, σ, T)
+
+Barrier option pricing formulas following Haug.
+
+String `option_type` is of the form
+
+    [D|U][O|I][C|P]
+
+"""
+function black_scholes_barrier_price(X, H, K, option_type::String, S, DF_r, DF_b, σ, T)
+    @assert length(option_type) == 3
+    o_type = uppercase(option_type)
+    @assert o_type[1] in ('D', 'U')
+    @assert o_type[2] in ('O', 'I')
+    @assert o_type[3] in ('C', 'P')
+    if o_type[1] == 'D'
+        η = 1
+    else
+        η = -1
+    end
+    if o_type[2] == 'O'
+        χ = 1
+    else
+        χ = -1
+    end
+    if o_type[3] == 'C'
+        ϕ = 1
+    else
+        ϕ = -1
+    end
+    return black_scholes_barrier_price(X, H, K, η, χ, ϕ, S, DF_r, DF_b, σ, T)
+end
diff --git a/src/utils/BrownianBridge.jl b/src/utils/BrownianBridge.jl
new file mode 100644
index 00000000..b770d6d3
--- /dev/null
+++ b/src/utils/BrownianBridge.jl
@@ -0,0 +1,89 @@
+
+
+
+"""
+    up_hit_probability(u, w0, w1, ν²)
+
+Return the probablity that a Brownian motion W(t) hits an upper barrer level `u`
+given that W(t0) = w0, w(t1) = w1 and variance of W(t1) conditional on W(T0) is
+ν².
+
+Formula is derived from S. Shreve, Stochastic Calculus for Finance II,
+Corollary 3.7.4.
+"""
+function up_hit_probability(u, w0, w1, ν²)
+    return exp.(-2.0 .* (u .- min.(w0, u)) .* (u .- min.(w1, u)) ./ ν²)
+end
+
+
+"""
+    down_hit_probability(d, w0, w1, ν²)
+
+Return the probablity that a Brownian motion W(t) hits a lower barrer level `d`
+given that W(t0) = w0, w(t1) = w1 and variance of W(t1) conditional on W(T0) is
+ν².
+
+Formula is derived from S. Shreve, Stochastic Calculus for Finance II,
+Corollary 3.7.4. and applied to -W(t).
+"""
+function down_hit_probability(d, w0, w1, ν²)
+    return exp.(-2.0 .* (max.(w0, d) .- d) .* (max.(w1, d) .- d) ./ ν²)
+end
+
+
+"""
+    barrier_no_hit_probability(
+        barrier_level,
+        barrier_direction,
+        brownian_levels::AbstractMatrix,
+        brownian_variances::AbstractMatrix,
+        )
+
+Calculate the (path-wise) probability that a Brownian motion W(t) does
+*not* hit an upper or lower barrier level given a list of observations.
+
+`barrier_level` is a scalar value.
+
+`barrier_direction` is +1 for down-barrier and -1 for up-barrier.
+
+`brownian_levels` is of size (p, n). Here, p is the number of
+simulated Brownian motion paths and n is the number of discrete
+observations.
+
+`brownian_variances` is of size (p, n-1). That is, we allow path-
+dependent variance. For non-path-dependence variance, use a
+matrix with p=1.
+"""
+function barrier_no_hit_probability(
+    barrier_level,
+    barrier_direction,
+    brownian_levels::AbstractMatrix,
+    brownian_variances::AbstractMatrix,
+    )
+    # ensure broadcasting
+    p = max(size(brownian_levels, 1), size(brownian_variances, 1))
+    @assert size(brownian_levels, 1) == p || size(brownian_levels, 1) == 1
+    @assert size(brownian_variances, 1) == p || size(brownian_variances, 1) == 1
+    #
+    @assert size(brownian_levels, 2) ≥ 2
+    @assert size(brownian_levels, 2) == size(brownian_variances, 2) + 1
+    #
+    @assert barrier_direction in (-1, 1)
+    #
+    if barrier_direction == -1
+        hit_probability_func = up_hit_probability
+    else
+        hit_probability_func = down_hit_probability
+    end
+    no_hit = ones(p)
+    for k in 1:(size(brownian_levels, 2) - 1)
+        hit_prob = hit_probability_func(
+            barrier_level,
+            brownian_levels[:, k],
+            brownian_levels[:, k+1],
+            brownian_variances[:, k],
+        )
+        no_hit = no_hit .* (1.0 .- hit_prob)
+    end
+    return no_hit
+end
diff --git a/test/componenttests/scenarios/asset_options.jl b/test/componenttests/scenarios/asset_options.jl
index 5e80d1ae..3a029b63 100644
--- a/test/componenttests/scenarios/asset_options.jl
+++ b/test/componenttests/scenarios/asset_options.jl
@@ -50,4 +50,30 @@ using UnicodePlots
         plot_scens(mv, "EUR-USD option mv")
     end
 
+    @testset "Test BarrierAssetOption simulation" begin
+        model = TestModels.hybrid_model_full
+        ch = TestModels.ch_full
+        times = 0.0:1.0:6.0
+        n_paths = 2^13
+        sim = DiffFusion.simple_simulation(model, ch, times, n_paths, with_progress_bar = false)
+        path = DiffFusion.path(sim, TestModels.ts_list, TestModels.context, DiffFusion.LinearPathInterpolation)
+        #
+        uoc_leg = DiffFusion.cashflow_leg(
+            "leg/UOC/EUR-USD/5y/1.25/1.85",
+            [ DiffFusion.BarrierAssetOptionFlow(5.0, 5.0, 1.25, 1.85, 0.0, "UOC", "EUR-USD", 0.5), ],
+            [ 1.0, ],
+            "USD"
+        )
+        uic_leg = DiffFusion.cashflow_leg(
+            "leg/UIC/EUR-USD/5y/1.25/1.85",
+            [ DiffFusion.BarrierAssetOptionFlow(5.0, 5.0, 1.25, 1.85, 0.0, "UIC", "EUR-USD", 0.5), ],
+            [ 1.0, ],
+            "USD"
+        )
+        scens = DiffFusion.scenarios([uoc_leg, uic_leg], times, path, "USD")
+        mv = DiffFusion.aggregate(scens, true, false)
+        #
+        plot_scens(mv, "EUR-USD barrier option mv")
+    end
+
 end
diff --git a/test/unittests/paths/path.jl b/test/unittests/paths/path.jl
index 65825b27..bf96d3fb 100644
--- a/test/unittests/paths/path.jl
+++ b/test/unittests/paths/path.jl
@@ -13,6 +13,7 @@ using Test
         @test_throws ErrorException DiffFusion.compounding_factor(NoPath(), 5.0, 8.0, 10.0, "Std")
         @test_throws ErrorException DiffFusion.asset(NoPath(), 5.0, "Std")
         @test_throws ErrorException DiffFusion.forward_asset(NoPath(), 5.0, 10.0, "Std")
+        @test_throws ErrorException DiffFusion.forward_asset_and_zero_bonds(NoPath(), 5.0, 10.0, "Std")
         @test_throws ErrorException DiffFusion.forward_index(NoPath(), 5.0, 10.0, "Std")
         @test_throws ErrorException DiffFusion.future_index(NoPath(), 5.0, 10.0, "Std")
         @test_throws ErrorException DiffFusion.fixing(NoPath(), 5.0, "Std")
@@ -302,11 +303,19 @@ end
         P_d = DiffFusion.zero_bond(p, t, T, "USD")
         P_f = DiffFusion.zero_bond(p, t, T, "EUR")
         @test isapprox(DiffFusion.forward_asset(p, t, T, "EUR-USD"), S_t .* P_f ./ P_d, atol=5.0e-15)
+        (S_t_, P_d_, P_f_) = DiffFusion.forward_asset_and_zero_bonds(p, t, T, "EUR-USD")
+        @test isapprox(S_t_, S_t, atol=1.0e-14)
+        @test isapprox(P_d_, P_d, atol=1.0e-14)
+        @test isapprox(P_f_, P_f, atol=1.0e-14)
         #
         S_t = DiffFusion.asset(p, t, "SXE50")
         P_d = DiffFusion.zero_bond(p, t, T, "EUR")
         P_f = DiffFusion.zero_bond(p, t, T, "SXE50")
         @test isapprox(DiffFusion.forward_asset(p, t, T, "SXE50"), S_t .* P_f ./ P_d, atol=5.0e-15)
+        (S_t_, P_d_, P_f_) = DiffFusion.forward_asset_and_zero_bonds(p, t, T, "SXE50")
+        @test isapprox(S_t_, S_t, atol=1.0e-11)
+        @test isapprox(P_d_, P_d, atol=1.0e-14)
+        @test isapprox(P_f_, P_f, atol=1.0e-14)
         #
         @test DiffFusion.fixing(p, -1.0, "SOFR") == 0.0123 * ones(5)
         @test DiffFusion.fixing(p,  0.0, "SOFR") == 0.0123 * ones(5)
@@ -399,6 +408,14 @@ end
         @test isapprox(DiffFusion.forward_asset(p, t, T, "EUR-USD"), 1.25 * exp(0.01*T) * ones(1), atol=5.0e-15)
         @test isapprox(DiffFusion.forward_asset(p, t, T, "SXE50"), 3750.00 * exp(0.01*T) * ones(1), atol=5.0e-15)
         #
+        (S_t_, P_d_, P_f_) = DiffFusion.forward_asset_and_zero_bonds(p, t, T, "EUR-USD")
+        S_T = DiffFusion.forward_asset(p, t, T, "EUR-USD")
+        @test isapprox(S_t_ .* P_f_ ./ P_d_, S_T, atol=1.0e-14)
+        #
+        (S_t_, P_d_, P_f_) = DiffFusion.forward_asset_and_zero_bonds(p, t, T, "SXE50")
+        S_T = DiffFusion.forward_asset(p, t, T, "SXE50")
+        @test isapprox(S_t_ .* P_f_ ./ P_d_, S_T, atol=1.0e-12)
+        #
         @test isapprox(DiffFusion.forward_index(p, 4.0, 5.0, "HICP"), 1.25 * ones(1), atol=5.0e-15)
         #
         @test isapprox(DiffFusion.future_index(p, 4.0, 5.0, "NIK"), 23776.5 * ones(1), atol=5.0e-15)
diff --git a/test/unittests/paths/simulated_path.jl b/test/unittests/paths/simulated_path.jl
index fc7bc430..12fd723b 100644
--- a/test/unittests/paths/simulated_path.jl
+++ b/test/unittests/paths/simulated_path.jl
@@ -134,6 +134,11 @@ using Test
                 zcb_f = DiffFusion.zero_bond(p, t, t+dT, "EUR")
                 one = fx_fwd ./ (fx .* zcb_f ./ zcb_d)
                 @test maximum(abs.(one .- 1.0)) < 1.5e-15
+                #
+                (S_t_, P_d_, P_f_) = DiffFusion.forward_asset_and_zero_bonds(p, t, t+dT, "EUR-USD")
+                @test maximum(abs.(S_t_ .- fx)) < 1.0e-15
+                @test maximum(abs.(P_d_ .- zcb_d)) < 1.0e-15
+                @test maximum(abs.(P_f_ .- zcb_f)) < 1.0e-15
             end
         end
     end
@@ -193,6 +198,11 @@ using Test
                 zcb_f = DiffFusion.zero_bond(p, t, t+dT, "EUR")
                 one = fx_fwd ./ (fx .* zcb_f ./ zcb_d)
                 @test maximum(abs.(one .- 1.0)) < 1.5e-15
+                #
+                (S_t_, P_d_, P_f_) = DiffFusion.forward_asset_and_zero_bonds(p, t, t+dT, "EUR-USD")
+                @test maximum(abs.(S_t_ .- fx)) < 1.0e-15
+                @test maximum(abs.(P_d_ .- zcb_d)) < 1.0e-15
+                @test maximum(abs.(P_f_ .- zcb_f)) < 1.0e-15
             end
         end
     end
diff --git a/test/unittests/payoffs/barrier_options.jl b/test/unittests/payoffs/barrier_options.jl
new file mode 100644
index 00000000..66675709
--- /dev/null
+++ b/test/unittests/payoffs/barrier_options.jl
@@ -0,0 +1,68 @@
+
+using DiffFusion
+
+using Test
+
+@testset "Test barrier option payoffs" begin
+
+    if !@isdefined(TestModels)
+        include("../../test_models.jl")
+    end
+
+    @testset "Barrier options setup" begin
+        F = DiffFusion.ForwardAsset(2.0, 5.0, "EUR-USD")
+        K = DiffFusion.Fixed(1.25)
+        H = DiffFusion.Fixed(1.75)
+        #
+        C = DiffFusion.BarrierAssetOption(F, K, H, "UOC", 0.0, 2)
+        P = DiffFusion.BarrierAssetOption(F, K, H, "UIP", 0.0, 2)
+        @test DiffFusion.obs_time(C) == 2.0
+        @test DiffFusion.obs_time(P) == 2.0
+        @test DiffFusion.obs_times(C) == union(Set(0.0), Set(2.0))
+        #
+        @test string(DiffFusion.BarrierAssetOption(F, K, H, "UOC", 0.0, 2)) == "UOCall(S(EUR-USD, 2.00, 5.00), X = 1.2500, H = 1.7500)"
+        @test string(DiffFusion.BarrierAssetOption(F, K, H, "UIC", 0.0, 2)) == "UICall(S(EUR-USD, 2.00, 5.00), X = 1.2500, H = 1.7500)"
+        @test string(DiffFusion.BarrierAssetOption(F, K, H, "DOC", 0.0, 2)) == "DOCall(S(EUR-USD, 2.00, 5.00), X = 1.2500, H = 1.7500)"
+        @test string(DiffFusion.BarrierAssetOption(F, K, H, "DIC", 0.0, 2)) == "DICall(S(EUR-USD, 2.00, 5.00), X = 1.2500, H = 1.7500)"
+        #
+        @test string(DiffFusion.BarrierAssetOption(F, K, H, "UOP", 0.0, 2)) == "UOPut(S(EUR-USD, 2.00, 5.00), X = 1.2500, H = 1.7500)"
+        @test string(DiffFusion.BarrierAssetOption(F, K, H, "UIP", 0.0, 2)) == "UIPut(S(EUR-USD, 2.00, 5.00), X = 1.2500, H = 1.7500)"
+        @test string(DiffFusion.BarrierAssetOption(F, K, H, "DOP", 0.0, 2)) == "DOPut(S(EUR-USD, 2.00, 5.00), X = 1.2500, H = 1.7500)"
+        @test string(DiffFusion.BarrierAssetOption(F, K, H, "DIP", 0.0, 2)) == "DIPut(S(EUR-USD, 2.00, 5.00), X = 1.2500, H = 1.7500)"
+    end
+
+
+    @testset "Barrier options valuation" begin
+        model = TestModels.hybrid_model_full
+        ch = TestModels.ch_full
+        times = [0.0, 2.0, 5.0]
+        n_paths = 8
+        sim = DiffFusion.simple_simulation(model, ch, times, n_paths, with_progress_bar = false)
+        path = DiffFusion.path(sim, TestModels.ts_list, TestModels.context, DiffFusion.LinearPathInterpolation)
+        #
+        F = DiffFusion.ForwardAsset(2.0, 5.0, "EUR-USD")
+        K = DiffFusion.Fixed(1.25)
+        C = DiffFusion.VanillaAssetOption(F, K, +1)(path)
+        P = DiffFusion.VanillaAssetOption(F, K, -1)(path)
+        for barrier_level in [ 1.25, 1.30, 1.35, 1.45 ]
+            H = DiffFusion.Fixed(barrier_level)
+            UOC = DiffFusion.BarrierAssetOption(F, K, H, "UOC", 0.0, 4)(path)
+            UIC = DiffFusion.BarrierAssetOption(F, K, H, "UIC", 0.0, 4)(path)
+            UOP = DiffFusion.BarrierAssetOption(F, K, H, "UOP", 0.0, 4)(path)
+            UIP = DiffFusion.BarrierAssetOption(F, K, H, "UIP", 0.0, 4)(path)
+            @test isapprox(UOC + UIC, C, atol=1.0e-14)
+            @test isapprox(UOP + UIP, P, atol=1.0e-14)
+        end
+        for barrier_level in [ 1.25, 1.20, 1.15, 1.05 ]
+            H = DiffFusion.Fixed(barrier_level)
+            DOC = DiffFusion.BarrierAssetOption(F, K, H, "DOC", 0.0, 4)(path)
+            DIC = DiffFusion.BarrierAssetOption(F, K, H, "DIC", 0.0, 4)(path)
+            DOP = DiffFusion.BarrierAssetOption(F, K, H, "DOP", 0.0, 4)(path)
+            DIP = DiffFusion.BarrierAssetOption(F, K, H, "DIP", 0.0, 4)(path)
+            @test isapprox(DOC + DIC, C, atol=1.0e-14)
+            @test isapprox(DOP + DIP, P, atol=1.0e-14)
+        end
+    end
+
+
+end
diff --git a/test/unittests/payoffs/payoffs.jl b/test/unittests/payoffs/payoffs.jl
index 1037bcc0..3eb4b6af 100644
--- a/test/unittests/payoffs/payoffs.jl
+++ b/test/unittests/payoffs/payoffs.jl
@@ -18,6 +18,7 @@ using Test
 
     include("amc_payoffs.jl")
     include("asset_options.jl")
+    include("barrier_options.jl")
     include("nodes_and_leafs.jl")
     include("rates_payoffs.jl")
     include("rates_options.jl")
diff --git a/test/unittests/products/asset_option_flows.jl b/test/unittests/products/asset_option_flows.jl
index cedaa509..30b772e8 100644
--- a/test/unittests/products/asset_option_flows.jl
+++ b/test/unittests/products/asset_option_flows.jl
@@ -25,4 +25,22 @@ using Test
         # println(string(DiffFusion.expected_amount(cf, 0.0)))
     end
 
+    @testset "Test BarrierAssetOptionFlow" begin
+        cf = DiffFusion.BarrierAssetOptionFlow(5.0, 6.0, 1.25, 1.35, 0.0, "UOC", "EUR-USD", 0.5)
+        @test DiffFusion.pay_time(cf) == 6.0
+        @test string(DiffFusion.expected_amount(cf, 0.0)) == "UOCall(S(EUR-USD, 0.00, 5.00), X = 1.2500, H = 1.3500)"
+        @test string(DiffFusion.expected_amount(cf, 2.0)) == "UOCall(S(EUR-USD, 2.00, 5.00), X = 1.2500, H = 1.3500)"
+        @test string(DiffFusion.expected_amount(cf, 5.0)) == "UOCall(S(EUR-USD, 5.00, 5.00), X = 1.2500, H = 1.3500)"
+        @test string(DiffFusion.expected_amount(cf, 6.0)) == "UOCall(S(EUR-USD, 5.00, 5.00), X = 1.2500, H = 1.3500)"
+        @test string(DiffFusion.expected_amount(cf, 7.0)) == "UOCall(S(EUR-USD, 5.00, 5.00), X = 1.2500, H = 1.3500)"
+        #
+        cf = DiffFusion.BarrierAssetOptionFlow(5.0, 6.0, 1.25, 1.35, 0.0, "DIP", "EUR-USD", 0.5)
+        @test DiffFusion.pay_time(cf) == 6.0
+        @test string(DiffFusion.expected_amount(cf, 0.0)) == "DIPut(S(EUR-USD, 0.00, 5.00), X = 1.2500, H = 1.3500)"
+        @test string(DiffFusion.expected_amount(cf, 2.0)) == "DIPut(S(EUR-USD, 2.00, 5.00), X = 1.2500, H = 1.3500)"
+        @test string(DiffFusion.expected_amount(cf, 5.0)) == "DIPut(S(EUR-USD, 5.00, 5.00), X = 1.2500, H = 1.3500)"
+        @test string(DiffFusion.expected_amount(cf, 6.0)) == "DIPut(S(EUR-USD, 5.00, 5.00), X = 1.2500, H = 1.3500)"
+        @test string(DiffFusion.expected_amount(cf, 7.0)) == "DIPut(S(EUR-USD, 5.00, 5.00), X = 1.2500, H = 1.3500)"
+    end
+
 end
\ No newline at end of file
diff --git a/test/unittests/utils/barriers.jl b/test/unittests/utils/barriers.jl
new file mode 100644
index 00000000..0d197af8
--- /dev/null
+++ b/test/unittests/utils/barriers.jl
@@ -0,0 +1,175 @@
+
+using DiffFusion
+
+using Test
+
+@testset "Black-Scholes single barrier pricing." begin
+# We test versus QuantLib values.
+
+    @testset "Vanilla Black-Scholes formula" begin
+        S = 100.0
+        X = [ 90.0, 100.0, 110.0 ]
+        T = 2.0
+        r = 0.08
+        b = 0.04
+        σ = 0.25
+        #
+        DF_r = exp(-r*T)
+        DF_b = exp(-b*T)
+        F = S ./ DF_b
+        #
+        C_BK = DF_r * DiffFusion.black_price(X, F, σ, T, +1)
+        P_BK = DF_r * DiffFusion.black_price(X, F, σ, T, -1)
+        C_BS = DiffFusion.black_scholes_vanilla_price(X, +1, S, DF_r, DF_b, σ * sqrt(T))
+        P_BS = DiffFusion.black_scholes_vanilla_price(X, -1, S, DF_r, DF_b, σ * sqrt(T))
+        @test C_BS == C_BK
+        @test P_BS == P_BK
+    end
+
+    @testset "Haug Barrier options" begin
+        # Haug, Table 4-13, p. 154
+        S = 100  # spot
+        K = 3    # rebate
+        T = 0.5  # time to maturity
+        r = 0.08
+        b = 0.04
+        σ1 = 0.25
+        σ2 = 0.30
+        #
+        DF_r = exp(-r*T)
+        DF_b = exp(-b*T)
+        #
+        data = [
+            # type,   X,   H,  ---  Price  ---
+            #                  σ=0.25,  σ=0.30
+            #
+            ["DOC",  90,  95,  9.0246,  8.8334, ],
+            ["DOC", 100,  95,  6.7924,  7.0285, ],
+            ["DOC", 110,  95,  4.8759,  5.4137, ],
+            ["DOC",  90, 100,  3.0000,  3.0000, ],
+            ["DOC", 100, 100,  3.0000,  3.0000, ],
+            ["DOC", 110, 100,  3.0000,  3.0000, ],
+            #
+            ["UOC",  90, 105,  2.6789,  2.6341, ],
+            ["UOC", 100, 105,  2.3580,  2.4389, ],
+            ["UOC", 110, 105,  2.3453,  2.4315, ],
+            #
+            #
+            ["DOP",  90,  95,  2.2798,  2.4170, ],
+            ["DOP", 100,  95,  2.2947,  2.4258, ],
+            ["DOP", 110,  95,  2.6252,  2.6246, ],
+            ["DOP",  90, 100,  3.0000,  3.0000, ],
+            ["DOP", 100, 100,  3.0000,  3.0000, ],
+            ["DOP", 110, 100,  3.0000,  3.0000, ],
+            #
+            ["UOP",  90, 105,  3.7760,  4.2293, ],
+            ["UOP", 100, 105,  5.4932,  5.8032, ],
+            ["UOP", 110, 105,  7.5187,  7.5649, ],
+            #
+            #
+            ["DIC",  90,  95,  7.7627,  9.0093, ],
+            ["DIC", 100,  95,  4.0109,  5.1370, ],
+            ["DIC", 110,  95,  2.0576,  2.8517, ],
+            ["DIC",  90, 100, 13.8333, 14.8816, ],
+            ["DIC", 100, 100,  7.8494,  9.2045, ],
+            ["DIC", 110, 100,  3.9795,  5.3043, ],
+            #
+            ["UIC",  90, 105, 14.1112, 15.2098, ],
+            ["UIC", 100, 105,  8.4482,  9.7278, ],
+            ["UIC", 110, 105,  4.5910,  5.8350, ],
+            #
+            #
+            ["DIP",  90,  95,  2.9586,  3.8769, ],
+            ["DIP", 100,  95,  6.5677,  7.7989, ],
+            ["DIP", 110,  95, 11.9752, 13.3078, ],
+            ["DIP",  90, 100,  2.2845,  3.3328, ],
+            ["DIP", 100, 100,  5.9085,  7.2636, ],
+            ["DIP", 110, 100, 11.6465, 12.9713, ],
+            #
+            ["UIP",  90, 105,  1.4653,  2.0658, ],
+            ["UIP", 100, 105,  3.3721,  4.4226, ],
+            ["UIP", 110, 105,  7.0846,  8.3686, ],
+        ]
+        for d in data
+            V1 = DiffFusion.black_scholes_barrier_price(d[2], d[3], K, d[1], S, DF_r, DF_b, σ1, T)
+            V2 = DiffFusion.black_scholes_barrier_price(d[2], d[3], K, d[1], S, DF_r, DF_b, σ2, T)
+            @test isapprox(V1, d[4], atol=1.0e-4)
+            @test isapprox(V2, d[5], atol=1.0e-4)
+        end
+    end
+
+    @testset "Broadcasting" begin
+        S = [ 90, 95, 100, 110, 120 ]
+        K = 0    # rebate
+        T = 0.5  # time to maturity
+        r = 0.08
+        b = 0.04
+        σ1 = 0.25
+        σ2 = 0.30
+        #
+        DF_r = exp(-r*T)
+        DF_b = exp(-b*T)
+        #
+        X = 100
+        H = 90
+        V = DiffFusion.black_scholes_barrier_price(X, H, K, "DOP", S, DF_r, DF_b, σ1, T)
+        V_ref = [
+            0.0,
+            0.1290974204294102,
+            0.2201391042455878,
+            0.2636825755986192,
+            0.19373974928092275
+        ]
+        @test isapprox(V, V_ref, atol=1.e-14)
+        #
+        S = 100
+        X = 90   # strike
+        H = 105  # barrier
+        K = 3    # rebate
+        σ = [ 0.25, 0.30 ]
+        V = DiffFusion.black_scholes_barrier_price(X, H, K, "UIP", S, DF_r, DF_b, σ, T)
+        V_ref = [
+            1.4653126853069676,
+            2.0658325935176096,
+        ]
+        @test isapprox(V, V_ref, atol=1.e-14)
+    end
+
+    @testset "Consistency with Vanilla option" begin
+        # Haug, Table 4-13, p. 154
+        S = 100  # spot
+        X = 100  # strike
+        K = 0    # rebate
+        T = 0.5  # time to maturity
+        r = 0.08
+        b = 0.04
+        σ1 = 0.25
+        σ2 = 0.30
+        #
+        DF_r = exp(-r*T)
+        DF_b = exp(-b*T)
+        #
+        C = DiffFusion.black_scholes_vanilla_price(X, +1, S, DF_r, DF_b, σ1 * sqrt(T))
+        P = DiffFusion.black_scholes_vanilla_price(X, -1, S, DF_r, DF_b, σ1 * sqrt(T))
+        for H in [ 95, 100, 110, 120, 150, ]
+            UOC = DiffFusion.black_scholes_barrier_price(X, H, K, "UOC", S, DF_r, DF_b, σ1, T)
+            UIC = DiffFusion.black_scholes_barrier_price(X, H, K, "UIC", S, DF_r, DF_b, σ1, T)
+            UOP = DiffFusion.black_scholes_barrier_price(X, H, K, "UOP", S, DF_r, DF_b, σ1, T)
+            UIP = DiffFusion.black_scholes_barrier_price(X, H, K, "UIP", S, DF_r, DF_b, σ1, T)
+            # println(string(UOC) * ", " * string(UIC))
+            # println(string(UOP) * ", " * string(UIP))
+            @test isapprox(UOC + UIC, C, atol=1.0e-14)
+            @test isapprox(UOP + UIP, P, atol=1.0e-14)
+        end
+        for H in [ 105, 100, 90, 80, 50, ]
+            DOC = DiffFusion.black_scholes_barrier_price(X, H, K, "DOC", S, DF_r, DF_b, σ1, T)
+            DIC = DiffFusion.black_scholes_barrier_price(X, H, K, "DIC", S, DF_r, DF_b, σ1, T)
+            DOP = DiffFusion.black_scholes_barrier_price(X, H, K, "DOP", S, DF_r, DF_b, σ1, T)
+            DIP = DiffFusion.black_scholes_barrier_price(X, H, K, "DIP", S, DF_r, DF_b, σ1, T)
+            # println(string(DOC) * ", " * string(DIC))
+            # println(string(DOP) * ", " * string(DIP))
+            @test isapprox(DOC + DIC, C, atol=1.0e-14)
+            @test isapprox(DOP + DIP, P, atol=1.0e-14)
+        end
+    end
+end
diff --git a/test/unittests/utils/brownian_bridge.jl b/test/unittests/utils/brownian_bridge.jl
new file mode 100644
index 00000000..08b6256d
--- /dev/null
+++ b/test/unittests/utils/brownian_bridge.jl
@@ -0,0 +1,62 @@
+
+using DiffFusion
+
+using Test
+
+@testset "Barrier hit probabilities via Brownian Bridge." begin
+
+    @testset "Hit funtions" begin
+        σ = 0.5
+        T = 10
+        ν² = σ^2 * T
+        #
+        @test DiffFusion.up_hit_probability(2.0, 1.0, 1.5, ν²) > 0 && DiffFusion.up_hit_probability(2.0, 1.0, 1.5, ν²) < 1
+        @test DiffFusion.up_hit_probability(1.8, 1.0, 1.5, ν²) > 0 && DiffFusion.up_hit_probability(1.8, 1.0, 1.5, ν²) < 1
+        @test DiffFusion.up_hit_probability(2.0, 1.0, 1.5, ν²) < DiffFusion.up_hit_probability(1.8, 1.0, 1.5, ν²)
+        @test DiffFusion.up_hit_probability(2.0, 1.0, 1.5, ν²) == DiffFusion.up_hit_probability(2.0, 1.5, 1.0, ν²)
+        @test DiffFusion.up_hit_probability(1.5, 1.0, 1.5, ν²) == 1.0
+        @test DiffFusion.up_hit_probability(1.2, 1.0, 1.5, ν²) == 1.0
+        @test DiffFusion.up_hit_probability(1.0, 1.0, 1.5, ν²) == 1.0
+        @test DiffFusion.up_hit_probability(0.5, 1.0, 1.5, ν²) == 1.0
+        @test DiffFusion.up_hit_probability(2.0, 1.0, 1.5, 1.0e-3) == 0.0
+        @test DiffFusion.up_hit_probability(1.5, 1.0, 1.5, 1.0e-3) == 1.0
+        #
+        @test DiffFusion.down_hit_probability(0.5, 1.0, 1.5, ν²) == DiffFusion.down_hit_probability(0.5, 1.5, 1.0, ν²)
+        @test DiffFusion.down_hit_probability(0.5, 1.0, 1.5, ν²) == DiffFusion.up_hit_probability(2.0, 1.0, 1.5, ν²)
+        @test DiffFusion.down_hit_probability(1.0, 1.0, 1.5, ν²) == 1.0
+        @test DiffFusion.down_hit_probability(1.2, 1.0, 1.5, ν²) == 1.0
+        @test DiffFusion.down_hit_probability(1.5, 1.0, 1.5, ν²) == 1.0
+        @test DiffFusion.down_hit_probability(2.0, 1.0, 1.5, ν²) == 1.0
+        @test DiffFusion.down_hit_probability(0.5, 1.0, 1.5, 1.0e-3) == 0.0
+        @test DiffFusion.down_hit_probability(1.0, 1.0, 1.5, 1.0e-3) == 1.0
+        # broadcasting
+        w1 = [ 0.0, 1.0, 2.0, 3.0 ]
+        w2 = [ 0.5, 1.5, 2.5, 3.5 ]
+        pu = [0.09071795328941251, 0.6703200460356393, 1.0, 1.0]
+        pd = [1.0, 1.0, 1.0, 0.30119421191220214]
+        @test isapprox(DiffFusion.up_hit_probability(2.0, w1, w2, ν²), pu, atol=1.0e-14)
+        @test isapprox(DiffFusion.down_hit_probability(2.0, w1, w2, ν²), pd, atol=1.0e-14)
+    end
+
+    @testset "No-hit probabilites" begin
+        σ = 0.5
+        T = 10
+        ν² = σ^2 * T
+        #
+        h = 2.0 # barrier level
+        w = [ 0.0, 1.0, 2.0, 3.0 ]
+        pu = DiffFusion.up_hit_probability(h, w, w, ν²)
+        pd = DiffFusion.down_hit_probability(h, w, w, ν²)
+        #
+        W = w * ones(3)'
+        V = ν² * ones(4, 2)
+        p_no_up = DiffFusion.barrier_no_hit_probability(h, -1, W, V)
+        p_no_do = DiffFusion.barrier_no_hit_probability(h, +1, W, V)
+        @test p_no_up == (1.0 .- pu) .* (1.0 .- pu)
+        @test p_no_do == (1.0 .- pd) .* (1.0 .- pd)
+        # broadcasting
+        @test DiffFusion.barrier_no_hit_probability(h, -1, W, ν² * ones(1, 2)) == p_no_up
+        @test DiffFusion.barrier_no_hit_probability(h, +1, W, ν² * ones(1, 2)) == p_no_do
+    end
+
+end
\ No newline at end of file
diff --git a/test/unittests/utils/utils.jl b/test/unittests/utils/utils.jl
index 980fe942..ba98b86b 100644
--- a/test/unittests/utils/utils.jl
+++ b/test/unittests/utils/utils.jl
@@ -4,7 +4,9 @@ using Test
 @testset "Commonly used utility methods." begin
 
     include("bachelier.jl")
+    include("barriers.jl")
     include("black.jl")
+    include("brownian_bridge.jl")
     include("polynomialregression.jl")
     include("piecewiseregression.jl")