Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

switch to DataInterpolations instead of BSplineKit #141

Draft
wants to merge 2 commits into
base: master
Choose a base branch
from
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 4 additions & 3 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -4,26 +4,27 @@ authors = ["Alec Loudenback <[email protected]> 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"
Expand Down
4 changes: 3 additions & 1 deletion src/Yields.jl
Original file line number Diff line number Diff line change
@@ -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
Expand Down
79 changes: 22 additions & 57 deletions src/utils.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand All @@ -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


Expand All @@ -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])
Expand All @@ -95,70 +96,34 @@ 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
zero_vec = -log.(clamp.(discount_vec,0.00001,1)) ./ maturities
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)
Expand Down