some progress

src/fit.jl

@@ -130,7 +130,8 @@ __default_optic(m::MyModel) = OptArgs([
 
 """
 __default_optic(m::Yield.Constant) = OptArgs(@optic(_.rate.value) => -1.0 .. 1.0)
-__default_optic(m::Yield.MonotoneConvex) = OptArgs(@optic(_.rates[*]) => -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
@@ -282,18 +283,19 @@ function fit(mod0, quotes, method::F=Fit.Loss(x -> x^2);
 
 end
 
-function fit(mod0::Yield.MonotoneConvexUnInit, quotes, method;
+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.contract) for q in quotes]
+    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)
+    rates = fit(mc, quotes, method; variables, optimizer)
 end
 
 

src/model/Yield/MonotoneConvex.jl

@@ -19,13 +19,19 @@ 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)
@@ -148,25 +154,35 @@ function __monotone_convex_init(t, rates, times)
 
     return t, i_time, rates, times
 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ᵈ = map(2:length(times)) do i
-        (times[i] * rates[i] - times[i-1] * rates[i-1]) / (times[i] - times[i-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
-    pushfirst!(fᵈ, 0)
     # step 2
-    f = map(2:length(times)-1) do i
-        (times[i] - times[i-1]) / (times[i+1] - times[i-1]) * fᵈ[i+1] +
-        (times[i+1] - times[i]) / (times[i+1] - times[i-1]) * fᵈ[i]
+    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)
-    pushfirst!(f, fᵈ[2] - 0.5 * (f[1] - fᵈ[2]))
-    fᵈ[end], f[end-1], fᵈ[end]
-    push!(f, fᵈ[end] - 0.5 * (f[end] - fᵈ[end]))
+    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])
 
@@ -176,7 +192,8 @@ function __monotone_convex_fs(rates, times)
 
     return f, fᵈ
 end
-function myzero(t, rates, times)
+
+function myzero(t, rates, times, f, fᵈ)
     lt = last(times)
     # if the time is greater than the last input time then extrapolate using the forwards
     if t > lt
@@ -185,7 +202,6 @@ function myzero(t, rates, times)
     end
 
     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])
     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)
@@ -196,22 +212,43 @@ function myzero(t, rates, times)
 end
 
 function Base.zero(mc::MonotoneConvex, t)
+    lt = last(mc.times)
+    f, fᵈ = mc.f, mc.fᵈ
+    t, i_time, rates, times = __monotone_convex_init(t, mc.rates, mc.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
+
+    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 Continuous(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::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)
+        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ᵈ = mc.f, mc.fᵈ
+    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::MonotoneConvex, t)
+function FinanceCore.discount(mc::MonotoneConvexUnoptimized, t)
     r = zero(mc, t)
     return exp(-r * t)
 end

test/MonotoneConvex.jl

@@ -0,0 +1,53 @@
+@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 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
+
+    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
+
+
+    # https://repository.up.ac.za/bitstream/handle/2263/25882/dissertation.pdf?sequence=1&isAllowed=y
+end