Find parameters from the Vasicek model

Question

I am given the following bond: and need to fit the Vasicek model to this data.

My attempt is the following:

# ...  imports
years = np.array([1, 2, 3, 4, 7, 10])
pric = np.array([0, .93, .85, .78, .65, .55, .42])


X = sympy.symbols("a b sigma")
a, b, s = X
rt1_rt = np.diff(pric)
ab_rt = np.array([a*(b-r) for r in pric[1:] ])
term = rt1_rt - ab_rt

def normpdf(x, mean, sd):
    var = sd**2
    denom = (2*sym.pi*var)**.5
    num = sym.E**(-(x-mean)**2/(2*var))
    return num/denom

pdfs = np.array([sym.log(normpdf(x, 0, s)) for x in term])
func = 0
for el in pdfs:
    func += el
func = func.factor()
lmd = sym.lambdify(X, func)
def target_fun(params):
    return lmd(*params)
result = scipy.optimize.least_squares(target_fun, [10, 10, 10])

I don't think that it outputs correct solution.

Answer 1

Your code is almost correct. You want to maximize your function, therefore you need to place minus sign in front of lmd in your function.

def target_fun(params):
        return -lmd(*params)

Additionally, the initial values are usually set to less than 1. Picking 10 is not the best choice as the algorithm might converge to a saddle point. Consider [0.01, 0.01, 0.01].