diff --git a/Project.toml b/Project.toml index aa9d4dad..56b87f15 100644 --- a/Project.toml +++ b/Project.toml @@ -4,26 +4,27 @@ authors = ["Alec Loudenback and contributors"] version = "3.3.1" [deps] -BSplineKit = "093aae92-e908-43d7-9660-e50ee39d5a0a" +DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" FinanceCore = "b9b1ffdd-6612-4b69-8227-7663be06e089" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LsqFit = "2fda8390-95c7-5789-9bda-21331edee243" +NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" Optim = "429524aa-4258-5aef-a3af-852621145aeb" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228" [compat] -BSplineKit = "^0.8" FinanceCore = "^1" ForwardDiff = "^0.10" LsqFit = "^0.12" Optim = "^1" +Reexport = "^1.2" Roots = "^1" UnicodePlots = "^2" julia = "^1.6" -Reexport = "^1.2" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/src/Yields.jl b/src/Yields.jl index 9f540638..8419b76f 100644 --- a/src/Yields.jl +++ b/src/Yields.jl @@ -1,11 +1,13 @@ module Yields using Reexport +using StaticArrays using FinanceCore using FinanceCore: Rate, rate, discount, accumulation, Periodic, Continuous, forward @reexport using FinanceCore: Rate, rate, discount, accumulation, Periodic, Continuous, forward -import BSplineKit +import DataInterpolations import ForwardDiff +using NonlinearSolve using LinearAlgebra using UnicodePlots import LsqFit diff --git a/src/utils.jl b/src/utils.jl index cd502339..8005ee55 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -50,7 +50,7 @@ struct LinearSpline <: InterpolationKind end Bootstrap the rates with the given maturities, treating the rates according to the periodic frequencies in settlement_frequency. `interpolator` is any function that will take two vectors of inputs and output points and return a function that will estimate an output given a scalar input. That is -`interpolator` should be: `interpolator(xs, ys) -> f(x)` where `f(x)` is the interpolated value of `y` at `x`. +`interpolator` should be: `interpolator(ys, xs) -> f(x)` where `f(x)` is the interpolated value of `y` at `x`. Built in `interpolator`s in Yields are: - `QuadraticSpline()`: Quadratic spline interpolation. @@ -70,11 +70,11 @@ end # dispatch on the user-exposed InterpolationKind to the right # internally named interpolation function function _bootstrap_choose_interp(rates, maturities, settlement_frequency, i::QuadraticSpline) - return _bootstrap_inner(rates, maturities, settlement_frequency, cubic_interp) + return _bootstrap_inner(rates, maturities, settlement_frequency, DataInterpolations.QuadraticSpline) end function _bootstrap_choose_interp(rates, maturities, settlement_frequency, i::LinearSpline) - return _bootstrap_inner(rates, maturities, settlement_frequency, linear_interp) + return _bootstrap_inner(rates, maturities, settlement_frequency, DataInterpolations.LinearInterpolation) end @@ -85,6 +85,7 @@ function _bootstrap_inner(rates, maturities, settlement_frequency, interpolation discount_vec[1] = discount(Constant(rates[1]), maturities[1]) for t = 2:length(maturities) + @show t if isnothing(settlement_frequency[t]) # no settlement before maturity discount_vec[t] = discount(Constant(rates[t]), maturities[t]) @@ -95,13 +96,27 @@ function _bootstrap_inner(rates, maturities, settlement_frequency, interpolation cfs[end] += 1 function pv(v_guess) - v = interpolation_function([[0.0]; maturities[1:t]], vcat(1.0, discount_vec[1:t-1], v_guess...)) + @show "here1" + @show typeof(v_guess) + pv_inner(v_guess) + end + + function pv(v_guess::U{T}) where {T<:NonlinearSolve.ForwardDiff.Dual,U<:AbstractArray} + @show "here2" + pv_inner(v_guess.value) + end + + function pv_inner(v_guess) + v = interpolation_function(vcat(1.0, discount_vec[1:t-1], only(v_guess)),[[0.0]; maturities[1:t]]) return sum(v.(times) .* cfs) end target_pv = sum(map(t2 -> discount(Constant(rates[t]), t2), times) .* cfs) - root_func(v_guess) = pv(v_guess) - target_pv - root_func′(v_guess) = ForwardDiff.derivative(root_func, v_guess) - discount_vec[t] = solve(root_func, root_func′, rate(rates[t])) + f(u,p) = pv(u) .- p + # root_func(v_guess) = pv(v_guess) - target_pv + # root_func′(v_guess) = ForwardDiff.derivative(root_func, v_guess) + # discount_vec[t] = solve(root_func, root_func′, rate(rates[t])) + probN = NonlinearProblem{false}(f,@SVector[rate(rates[t])],target_pv) + discount_vec[t] = NonlinearSolve.solve(probN,NewtonRaphson(),tol=1e-9) end end @@ -109,56 +124,6 @@ function _bootstrap_inner(rates, maturities, settlement_frequency, interpolation return interpolation_function([0.0; maturities], [first(zero_vec); zero_vec]) end -# the ad-hoc approach to extrapoliatons is based on suggestion by author of -# BSplineKit at https://github.com/jipolanco/BSplineKit.jl/issues/19 -# this should not be exposed directly to user -struct _Extrap{I,L,R} - int::I # the BSplineKit interpolation - left::L # a tuple of (boundary, extrapolation function) - right::R # a tuple of (boundary, extrapolation function) -end - -function _wrap_spline(itp) - - S = BSplineKit.spline(itp) # spline passing through data points - B = BSplineKit.basis(S) # B-spline basis - - a, b = BSplineKit.boundaries(B) # left and right boundaries - - # For now, we construct the full spline S′(x). - # There are faster ways of doing this that should be implemented... - S′ = diff(S, BSplineKit.Derivative(1)) - - return _Extrap(itp, - (boundary = a, func = x->S(a) + S′(a)*(x-a)), - (boundary = b, func = x->S(b) + S′(b)*(x-b)), - - ) -end - -function _interp(e::_Extrap,x) - if x <= e.left.boundary - return e.left.func(x) - elseif x >= e.right.boundary - return e.right.func(x) - else - return e.int(x) - end -end - -function linear_interp(xs, ys) - int = BSplineKit.interpolate(xs, ys, BSplineKit.BSplineOrder(2)) - e = _wrap_spline(int) - return x -> _interp(e, x) -end - -function cubic_interp(xs, ys) - order = min(length(xs),3) # in case the length of xs is less than the spline order - int = BSplineKit.interpolate(xs, ys, BSplineKit.BSplineOrder(order)) - e = _wrap_spline(int) - return x -> _interp(e, x) -end - # used to display simple type name in show method # https://stackoverflow.com/questions/70043313/get-simple-name-of-type-in-julia?noredirect=1#comment123823820_70043313 name(::Type{T}) where {T} = (isempty(T.parameters) ? T : T.name.wrapper)