diff --git a/Project.toml b/Project.toml
index 2d672346..315cc3ff 100644
--- a/Project.toml
+++ b/Project.toml
@@ -14,6 +14,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
 OptimizationMetaheuristics = "3aafef2f-86ae-4776-b337-85a36adf0b55"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
@@ -30,17 +31,17 @@ FinanceModelsMakieCoreExt = "MakieCore"
 AccessibleOptimization = "^0.1.1"
 Accessors = "^0.1"
 BSplineKit = "^0.16"
+Dates = "^1.6"
 FinanceCore = "^2.1"
 IntervalSets = "^0.7"
 LinearAlgebra = "^1.6"
-Dates = "^1.6"
 MakieCore = "0.6"
 Optimization = "^3.15"
 OptimizationMetaheuristics = "^0.1.2"
 PrecompileTools = "^1.1"
 Reexport = "^1.2"
-StaticArrays = "^1.6"
 SpecialFunctions = "2"
+StaticArrays = "^1.6"
 Transducers = "^0.4"
 UnicodePlots = "^3.6"
 julia = "1.9"
diff --git a/src/FinanceModels.jl b/src/FinanceModels.jl
index 812f5d96..0f81d88d 100644
--- a/src/FinanceModels.jl
+++ b/src/FinanceModels.jl
@@ -15,6 +15,7 @@ import BSplineKit
 import UnicodePlots
 using Transducers: @next, complete, __foldl__, asfoldable
 import SpecialFunctions
+import QuadGK
 
 
 
diff --git a/src/fit.jl b/src/fit.jl
index b80e07f6..76489279 100644
--- a/src/fit.jl
+++ b/src/fit.jl
@@ -130,6 +130,8 @@ __default_optic(m::MyModel) = OptArgs([
 
 """
 __default_optic(m::Yield.Constant) = OptArgs(@optic(_.rate.value) => -1.0 .. 1.0)
+__default_optic(m::Yield.MonotoneConvexUnInit) = OptArgs(@optic(_.rates) .=> -1.0 .. 1.0)
+__default_optic(m::Yield.MonotoneConvex) = OptArgs(@optic(_.rates) .=> -1.0 .. 1.0)
 __default_optic(m::Yield.IntermediateYieldCurve) = OptArgs(@optic(_.ys[end]) => 0.0 .. 1.0)
 __default_optic(m::Yield.NelsonSiegel) = OptArgs([
     @optic(_.τ₁) => 0.0 .. 100.0
@@ -222,7 +224,7 @@ The default solver is `ECA()` from Metahueristics.jl. This is a stochastic globa
 ## Defining the variables
 
 An arbitrarily complex model may be the object we intend to fit - how does `fit` know what free variables are able to be solved for within the given model?
-`variables` is a singlular or vector optic argument. What does this mean?
+`variables` is a singular or vector optic argument. What does this mean?
 - An optic (or "lens") is a way to define an accessor to a given object. Example:
 
 ```julia-repl
@@ -237,7 +239,7 @@ julia> lens(obj)
 "AA"
 ```
 An optic argument is a singular or vector of lenses with an optional range of acceptable parameters. For example, we might have a model as follows where we want 
-`fit` to optize parameters `a` and `b`:
+`fit` to optimize parameters `a` and `b`:
 
 ```julia
 struct MyModel <:FinanceModels.AbstractModel
@@ -252,7 +254,7 @@ __default_optic(m::MyModel) = OptArgs([
 ```
 In this way, fit know which arbitrary parameters in a given object may be modified. Technically, we are not modifying the immutable `MyModel`, but instead efficiently creating a new instance. This is enabled by [AccessibleOptimization.jl](https://gitlab.com/aplavin/AccessibleOptimization.jl).
 
-Note that not all opitmization algorithms want a bounded interval. In that case, simply leave off the paired range. The prior example would then become:
+Note that not all optimization algorithms want a bounded interval. In that case, simply leave off the paired range. The prior example would then become:
 
 ```julia
 __default_optic(m::MyModel) = OptArgs([
@@ -265,7 +267,7 @@ __default_optic(m::MyModel) = OptArgs([
 
 ## Additional Examples
 
-See the tutorials in the package documentation for FinanceModels.jl or the docstrings of FinanceModels.jl's avaiable model types.
+See the tutorials in the package documentation for FinanceModels.jl or the docstrings of FinanceModels.jl's available model types.
 """
 function fit(mod0, quotes, method::F=Fit.Loss(x -> x^2);
     variables=__default_optic(mod0),
@@ -281,6 +283,22 @@ function fit(mod0, quotes, method::F=Fit.Loss(x -> x^2);
 
 end
 
+function fit(mod0::Yield.MonotoneConvexUnInit, quotes, method::F=Fit.Loss(x -> x^2);
+    variables=__default_optic(mod0),
+    optimizer=__default_optim(mod0)
+) where
+{F<:Fit.Loss}
+    # use maturities as the times
+    # fit a vector of rates to the times
+    times = [maturity(q.instrument) for q in quotes]
+    if !iszero(first(times))
+        pushfirst!(times, zero(eltype(times)))
+    end
+    mc = mod0(times)
+    rates = fit(mc, quotes, method; variables, optimizer)
+end
+
+
 function fit(mod0::Spline.BSpline, quotes, method::Fit.Bootstrap)
     discount_vector = [0.0]
     times = [maturity(quotes[1])]
diff --git a/src/model/Yield.jl b/src/model/Yield.jl
index 89477892..e907e227 100644
--- a/src/model/Yield.jl
+++ b/src/model/Yield.jl
@@ -68,6 +68,7 @@ end
 
 include("Yield/SmithWilson.jl")
 include("Yield/NelsonSiegelSvensson.jl")
+include("Yield/MonotoneConvex.jl")
 
 ## Generic and Fallbacks
 """
diff --git a/src/model/Yield/MonotoneConvex.jl b/src/model/Yield/MonotoneConvex.jl
new file mode 100644
index 00000000..be23d5b7
--- /dev/null
+++ b/src/model/Yield/MonotoneConvex.jl
@@ -0,0 +1,263 @@
+"""
+An unfit Monotone Convex Yield Curve Model will simply have 
+
+
+"""
+struct MonotoneConvex{T,U} <: AbstractYieldModel
+    f::Vector{T}
+    fᵈ::Vector{T}
+    rates::Vector{T}
+    times::Vector{U}
+    # inner constructor ensures f consistency with rates at construction
+    function MonotoneConvex(rates::Vector{T}, times::Vector{U}) where {T,U}
+        f, fᵈ = __monotone_convex_fs(rates, times)
+        new{T,U}(f, fᵈ, rates, times)
+    end
+end
+
+
+struct MonotoneConvexUnInit
+end
+
+struct MonotoneConvexUnoptimized{T,U}
+    rates::Vector{T}
+    times::Vector{U}
+end
+
+MonotoneConvex() = MonotoneConvexUnInit()
+function (m::MonotoneConvexUnInit)(times)
+    rates = zeros(length(times))
+    MonotoneConvex(rates, times)
+end
+
+
+
+function __issector1(g0, g1)
+    a = (g0 > 0) && (g1 >= -2 * g0) && (-0.5 * g0 >= g1)
+    b = (g0 < 0) && (g1 <= -2 * g0) && (-0.5 * g0 <= g1)
+    a || b
+end
+
+function __issector2(g0, g1)
+    a = (g0 > 0) && (g1 < -2 * g0)
+    b = (g0 < 0) && (g1 > -2 * g0)
+    a || b
+end
+
+
+# Hagan West - WILMOTT magazine pgs 75-81
+function g(f⁻, f, fᵈ)
+    @show f⁻, f, fᵈ
+    @show g0 = f⁻ - fᵈ
+    @show g1 = f - fᵈ
+    if sign(g0) == sign(g1)
+        # sector (iv)
+        η = g1 / (g1 + g0)
+        α = -g0 * g1 / (g1 + g0)
+
+        if x < η
+            return x -> α + (g0 - α) * ((η - x) / η)^2
+        else
+            return x -> α + (g1 - α) * ((x - η) / (1 - η))^2
+        end
+
+
+    elseif __issector1(g0, g1)
+        # sector (i)
+        x -> g0 * (1 - 4 * x + 3 * x^2) + g1 * (-2 * x + 3 * x^2)
+    elseif __issector2(g0, g1)
+        # sector (ii)
+        η = (g1 + 2 * g0) / (g1 - g0)
+        if x < η
+            return x -> g0
+        else
+            return x -> g0 + (g1 - g0) * ((x - η) / (1 - η))^2
+        end
+    else
+        # sector (iii)
+        η = 3 * g1 / (g1 - g0)
+        if x > η
+            return x -> g1
+        else
+            return x -> g1 + (g0 - g1) * ((η - x) / η)^2
+        end
+
+    end
+end
+
+function g_rate(x, f⁻, f, fᵈ)
+    @show x, f⁻, f, fᵈ
+    @show g0 = f⁻ - fᵈ
+    @show g1 = f - fᵈ
+    if sign(g0) == sign(g1)
+        # sector (iv)
+        η = g1 / (g1 + g0)
+        α = -g0 * g1 / (g1 + g0)
+
+        if x < η
+            return α * x - (g0 - α) * ((η - x)^3 / η^2 - η) / 3
+        else
+            return (2 * α + g0) / 3 * η + α * (x - η) + (g1 - α) / 3 * (x - η)^3 / (1 - η)^2
+        end
+
+
+    elseif __issector1(g0, g1)
+        # sector (i)
+        @show "(i)"
+        g0 * (x - 2 * x^2 + x^3) + g1 * (-x^2 + x^3)
+    elseif __issector2(g0, g1)
+        # sector (ii)
+        η = (g1 + 2 * g0) / (g1 - g0)
+        if x < η
+            return g0 * x
+        else
+            return g0 * x + ((g1 - g0) * (x - η)^3 / (1 - η)^2) / 3
+        end
+    else
+        # sector (iii)
+        η = 3 * g1 / (g1 - g0)
+        if x > η
+            return (2 * g1 + g0) / 3 * η + g1 * (x - η)
+        else
+            return g1 * x - (g0 - g1) * ((η - x)^3 / η^2 - η) / 3
+        end
+
+    end
+end
+
+function forward(t, rates, times)
+
+    t, i_time, rates, times = __monotone_convex_init(t, rates, times)
+    f, fᵈ = __monotone_convex_fs(rates, times)
+
+    x = (t - times[i_time]) / (times[i_time+1] - times[i_time])
+    return fᵈ[i_time+1] + g(x, f[i_time], f[i_time+1], fᵈ[i_time+1])
+end
+
+"""
+    returns the index associated with the time t, an initial rate vector, and a time vector
+"""
+function __i_time(t, times)
+    i_time = findfirst(x -> x > t, times)
+    if i_time == nothing
+        i_time = lastindex(times)
+    end
+    return i_time
+end
+
+"""
+    returns a pair of vectors (f and fᵈ) used in Monotone Convex Yield Curve fitting
+"""
+function __monotone_convex_fs(rates, times)
+    # step 1
+    fᵈ = deepcopy(rates)
+    for i in 2:length(times)
+        fᵈ[i] = (times[i] * rates[i] - times[i-1] * rates[i-1]) / (times[i] - times[i-1])
+    end
+    # step 2
+    f = similar(rates, length(rates) + 1)
+    # fill in middle elements first, then do 1st and last
+    for i in 1:length(rates)-1
+        t_prior = if i == 1
+            0
+        else
+            times[i-1]
+        end
+
+        weight1 = (times[i] - t_prior) / (times[i+1] - t_prior)
+        weight2 = (times[i+1] - times[i]) / (times[i+1] - t_prior)
+        f[i+1] = weight1 * fᵈ[i+1] + weight2 * fᵈ[i]
+    end
+    # step 3
+    # collar(a,b,c) = clamp(b, a, c)
+    f[1] = fᵈ[1] - 0.5 * (f[2] - fᵈ[1])
+
+    f[end] = fᵈ[end] - 0.5 * (f[end-1] - fᵈ[end])
+    f[1] = clamp(f[1], 0, 2 * fᵈ[2])
+    f[end] = clamp(f[end], 0, 2 * fᵈ[end])
+
+    for j in 2:(length(times)-1)
+        f[j] = clamp(f[j], 0, 2 * min(fᵈ[j], fᵈ[j+1]))
+    end
+
+    return f, fᵈ
+end
+
+
+function Base.zero(mc::MonotoneConvex, t)
+    lt = last(mc.times)
+    f, fᵈ = mc.f, mc.fᵈ
+    if t > lt
+        r = Base.zero(mc, lt)
+        i_time = __i_time(t, mc.times)
+        return r * lt / t + forward(lt, mc.rates, mc.times) * (1 - lt / t)
+    end
+    @show i_time = __i_time(t, mc.times)
+    # if the time is greater than the last input time then extrapolate using the forwards
+
+    x = if i_time == 1
+        x = t / times[i_time]
+    else
+        x = (t - times[i_time-1]) / (times[i_time] - times[i_time-1])
+    end
+    G = g(f[i_time], f[i_time+1], fᵈ[i_time])
+    return Continuous(1 / t * (times[i_time] * rates[i_time] + (t - times[i_time]) * fᵈ[i_time] + (times[i_time] - times[i_time-1]) * G))
+
+    # STATUS:
+    # intermediate G/other results OK
+    # need to get the right rate. Instead of last formula, trying the approach on pg 39 of 
+    # http://uu.diva-portal.org/smash/get/diva2:1477828/FULLTEXT01.pdf 
+    # and have added QuadGK and converted the G function to return a function of x instead of a calculated value 
+    # (also need to change the signature of associated test cases)
+    # QuadGK.quadgk(G,t)
+
+end
+
+function FinanceCore.discount(mc::MonotoneConvex, t)
+    r = zero(mc, t)
+    return discount(r, t)
+end
+
+function Base.zero(mc::MonotoneConvexUnoptimized, t)
+    lt = last(times)
+    # if the time is greater than the last input time then extrapolate using the forwards
+    if t > lt
+        r = myzero(lt, rates, times, f, fᵈ)
+        return r * lt / t + forward(lt, rates, times) * (1 - lt / t)
+    end
+
+    t, i_time, rates, times = __monotone_convex_init(t, rates, times)
+    f, fᵈ = __monotone_convex_fs(mc.rates, mc.times)
+    x = (t - times[i_time]) / (times[i_time+1] - times[i_time])
+    G = g_rate(x, f[i_time], f[i_time+1], fᵈ[i_time+1])
+    return 1 / t * (times[i_time] * rates[i_time] + (t - times[i_time]) * fᵈ[i_time+1] + (times[i_time+1] - times[i_time]) * G)
+
+end
+
+function FinanceCore.discount(mc::MonotoneConvexUnoptimized, t)
+    r = zero(mc, t)
+    return exp(-r * t)
+end
+
+times = 1:5
+rates = [0.03, 0.04, 0.047, 0.06, 0.06]
+# forward(5.19, rates, times)
+# myzero(1, rates, times)
+
+
+# using Test
+# @test forward(0.5, rates, times) ≈ 0.02875
+# @test forward(1, rates, times) ≈ 0.04
+# @test forward(2, rates, times) ≈ 0.0555
+# @test forward(2.5, rates, times) ≈ 0.0571254591368226
+# @test forward(5, rates, times) ≈ 0.05025
+# @test forward(5.2, rates, times) ≈ 0.05025
+
+# @test myzero(0.5, rates, times) ≈ 0.02625
+# @test myzero(1, rates, times) ≈ 0.03
+# @test myzero(2, rates, times) ≈ 0.04
+# @test myzero(2.5, rates, times) ≈ 0.0431375956535047
+# @test myzero(5, rates, times) ≈ 0.06
+# @test myzero(5.2, rates, times) ≈ 0.059625
+
+
diff --git a/test/MonotoneConvex.jl b/test/MonotoneConvex.jl
new file mode 100644
index 00000000..7ffbe918
--- /dev/null
+++ b/test/MonotoneConvex.jl
@@ -0,0 +1,65 @@
+@testset "Monotone Convex" begin
+
+    # Interpolation of the yield curve, Gustaf Dehlbom
+    # http://uu.diva-portal.org/smash/get/diva2:1477828/FULLTEXT01.pdf
+
+    prices = [0.98, 0.955, 0.92, 0.88, 0.830]
+    times = [1, 2, 3, 4, 5]
+    quotes = ZCBPrice.(prices, times)
+    rates = @. -log(prices) / times
+    c = Yield.MonotoneConvex(rates, times)
+    # fit(Yield.MonotoneConvex(), quotes)
+
+    f, fᵈ = Yield.__monotone_convex_fs(rates, times)
+
+    @test all(f .== c.f)
+    @test all(fᵈ .== c.fᵈ)
+
+    @test fᵈ[1] ≈ 0.0202 atol = 0.0001
+    @test fᵈ[2] ≈ 0.0258 atol = 0.0001
+    @test fᵈ[3] ≈ 0.0373 atol = 0.0001
+    @test fᵈ[4] ≈ 0.0445 atol = 0.0001
+    @test fᵈ[5] ≈ 0.0585 atol = 0.0001
+
+    @test f[1] ≈ 0.0188 atol = 0.0001
+    @test f[2] ≈ 0.023 atol = 0.0001
+    @test f[3] ≈ 0.0316 atol = 0.0001
+    @test f[4] ≈ 0.0409 atol = 0.0001
+    @test f[5] ≈ 0.0515 atol = 0.0001
+    @test f[6] ≈ 0.0620 atol = 0.0001
+
+    @test FinanceModels.Yield.g(0.5, 0.018793076350927487, 0.023021969250703423, 0.020202707317519466) ≈ 0.0042 * (0.5)^2 - 0.0014 atol = 0.0001
+    @test FinanceModels.Yield.g(0, 0.023021969250703423, 0.03158945081076577, 0.02584123118388738) ≈ -0.0028 atol = 0.0001
+    @test FinanceModels.Yield.g(0.5, 0.023021969250703423, 0.03158945081076577, 0.02584123118388738) ≈ 0.0087 * (0.5)^2 − 0.0002 * 0.5 - 0.0028 atol = 0.0001
+    @test FinanceModels.Yield.g(0.5, 0.03158945081076577, 0.04089471650423902, 0.03733767043764417) ≈ -0.0063 * (0.5)^2 + 0.0156 * 0.5 - 0.0057 atol = 0.0001
+    @test FinanceModels.Yield.g(0.5, 0.04089471650423902, 0.051473984626221235, 0.04445176257083387) ≈ 0.0102 * (0.5)^2 + 0.0004 * 0.5 - 0.0036 atol = 0.0001
+    @test FinanceModels.Yield.g(0.5, 0.04089471650423902, 0.051473984626221235, 0.044451762570833877) ≈ -0.0105 * (0.5)^2 + 0.021 * 0.5 - 0.007 atol = 0.0001
+
+
+    function r(t)
+        if 0 <= t <= 1
+            return 0.0014t^2 + 0.0188
+        elseif 1 <= t <= 1.0233
+            return -0.0028 / t + 0.0230
+        elseif 1.0233 <= t <= 2
+            return 0.0029t^2 - 0.0088t - 0.0058 / t + 0.0319
+        elseif 2 <= t <= 3
+            return -0.0022t^2 + 0.0212t + 0.0324 / t - 0.0268
+        elseif 3 <= t <= 4
+            return 0.0031t^2 - 0.0274t - 0.1188 / t + 0.1217
+        elseif 4 <= t <= 5
+            return -0.0035t^2 + 0.0525t + 0.314 / t - 0.2005
+        else
+            error("t is out of the defined range")
+        end
+    end
+
+
+    @testset for t in range(0, 5, 30)
+        @test zero(c, t) ≈ r(t)
+    end
+
+
+    # TODO: More tests from
+    # https://repository.up.ac.za/bitstream/handle/2263/25882/dissertation.pdf?sequence=1&isAllowed=y
+end
\ No newline at end of file