Optimization Challenges in matching simulated data

Specific Domains Optimization (Mathematical)
1

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 (Newton, ECA, BFGS, NelderMead, DE, 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[3] + p[4] * (1.0 - exp(-t / p[1])) / (t / p[1]) + p[5] * ((1.0 - exp(-t / p[1])) / (t / p[1]) - exp(-t / p[1])) + p[6] * ((1.0 - exp(-t / p[2])) / (t / p[2]) - exp(-t / p[2]))

# 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
2

I didn’t run the code, but does it even find solutions that are approximately close? e.g., to 1e-3 or so?

Your function is non-convex, which seems like it’d make things very tricky.