I am trying to abstract some components of the Yields.jl package and having a hard time with a particular problem, which I’ve tried to distill down to a MWE. If it helps, the context is fitting a Nelson-Siegel-Svensson yield curve to observed prices.
My issue is that for a six-parameter function, I am not able to satisfactorily solve for the known parameters I used to generate the noise-less data. Ideally I would be able to deterministically achieve the target parameters to a very tight tolerance (e.g. parameters within something like 1e-8 or better of the known targets) but have found no routine that accomplishes this.
I have tried a variety of algorithms (
NewtonIP) with a variety of algorithm parameters. Bounded and unbounded as well, with autodiff and without.
Is this just a particularly difficult function to fit or do folks have recommendations?
using Optimization, OptimizationOptimJL, OptimizationMetaheuristics # the model, p are the parameters trying to be solved for f(t,p) = p + p * (1.0 - exp(-t / p)) / (t / p) + p * ((1.0 - exp(-t / p)) / (t / p) - exp(-t / p)) + p * ((1.0 - exp(-t / p)) / (t / p) - exp(-t / p)) # sample "real" parameters p_target = [0.777730,12.625607,1.193712/100,1.996922/100,2.953628/100,4.251043/100] # sample data ts = vcat([3/12, 6/12, 9/12], collect(1:30)) # goal is to mimimize difference between observed data which comes in the transformed form of `f` observed_data = map(t -> exp(-t*f(t, p_target)), ts) # untransformed data (note `_r` suffix): observed_data_r = map(t -> f(t, p_target), ts) x0 = [1.,1.,0.,0.,0.,0.] # upper and lower bounds if needed u = [100.,100.,10.,10.,10.,10.] l = [0.,0.,-10.,-10.,-10.,-10.] #minimize least squares on both transformed and untransformed data function objective(p,q) sum((observed_data[i] - exp(-ts[i]*f(ts[i], p)))^2 for i in 1:length(ts)) end function objective_r(p,q) sum((observed_data[i] - f(ts[i], p))^2 for i in 1:length(ts)) end fo = OptimizationFunction(objective,Optimization.AutoForwardDiff()) prob = OptimizationProblem(fo, x0) sol = solve(prob, BFGS()) sol.objective fo_r = OptimizationFunction(objective_r,Optimization.AutoForwardDiff()) prob_r = OptimizationProblem(fo_r, x0;lb=l, ub=u) sol_r = solve(prob_r, ECA(); maxiters=10000) sol_r.objective