Maximizing Log Likelihood Estimation in Python

Currently I am writing my thesis about the liquidity of green bonds. I am trying to compute the LOT measure based on Lesmond et al. (1999) and Chen et al. (2007). Therefore, I need to maximize the following log likelihood function

L(a 1,j,a 2,j, β j,1, β j,2,σ j|R j,t,∆Index)

in python (a screenshot is appended for better readability):

Φ i,j

denotes the cumulative distribution function for each bond-year evaluated at

L(a i,j− β j,1D j,t∗ ∆R f,t− β j,2D j,t∗ ∆Index t)/σ j

So far, I tried to apply the following code, but unfortunately received wrong estimates yet:

 import numpy as np 
 from scipy.optimize import minimize
 from scipy.stats import norm

  x = np.random.rand(200,7)
  y = np.random.rand(200,1)

  def negative_loglikehood(params,x):

     alpha_1 = params[0] 
     alpha_2 = params[1]
     beta_1 = params[2]
     beta_2 = params[4]

     rjt = x[:,0]
     djt = x[:,1]
     index = x[:,2]
     rft = x[:,3]

     term_0 = np.sum(np.log(1/(2*np.pi*sigma)**0.5))

     term_1 = np.sum( (1/ (2*sigma**2) ) * (np.sum(rjt+alpha_1-beta_1*djt*rft) - beta_2*djt*index )**2)

     term_2 = np.sum(np.log(1/(2*np.pi*sigma)**0.5))

     term_3 = np.sum( (1/ (2*sigma**2) ) * (np.sum(rjt+alpha_2-beta_1*djt*rft) - beta_2 *djt*index )**2)

     term_4 = np.sum(
            np.log(
              norm.cdf( y, loc=(alpha_1 - beta_1*djt*rft - beta_2*djt*index)/sigma, scale=sigma)
            - norm.cdf(y, loc=(alpha_1 - beta_1*djt*rft - beta_2*djt*index)/sigma, scale=sigma)
            ))

     LL1 = (term_0 - term_1 + term_2 - term_3 + term_4)
     return(LL1)
 
guess = [1,1,1,1,1]
results = minimize(negative_loglikehood, guess,args =(x), method = 'Nelder-Mead', options={'disp': True})
print(results.x)

I am not sure of the implementation of my code is correct.

There are several errors in your formula when compared to the analytic expression you have provided:

in term_2 = np.sum(np.log(1/(2*np.pi*sigma)**0.5)) sigma should be squared -> term_2 = np.sum(np.log(1/(2*np.pi*sigma**2)**0.5)).

This is also the case for term_0, change term_0 to term_0 = np.sum(np.log(1/(2*np.pi*sigma**2)**0.5))

You also seem to be minimizing the positive log-likelihood. This is the opposite of what you want to do. To minimize the negative LL change

return(LL1)

to

return(-LL1)

You also do not correctly evaluate the cdfs in term_4. In Lesmond 1999 it states that Φi,j are the cdfs of a standard normal distribution evaluated at (𝛼𝑖,𝑗−𝛽𝑗1Duration𝑗,𝑡∗Δ𝑅𝑓𝑡−𝛽𝑗2Duration𝑗,𝑡∗ΔS&P Index𝑡)/𝜎𝑗. However, what you have programmed is to evaluate a non-standard normal distribution with mean (𝛼𝑖,𝑗−𝛽𝑗1Duration𝑗,𝑡∗Δ𝑅𝑓𝑡−𝛽𝑗2Duration𝑗,𝑡∗ΔS&P Index𝑡)/𝜎𝑗 and variance sigma^2 at some value y. To evaluate the correct cdf, you need to change term_4 to

term_4 = np.sum(
            np.log(
              norm.cdf((alpha_2 - beta_1*djt*rft - beta_2*djt*index)/sigma, loc = 0, scale = 1)
            - norm.cdf((alpha_1 - beta_1*djt*rft - beta_2*djt*index)/sigma, loc = 0, scale = 1)
            ))

, where in the first term alpha_2 is used and in the second one alpha_1 as stated in Chen 2007.

In Chen 2007 it also states that α1 j, α2 j, βj, and σj are solved by maximizing the likelihood. You have yet not included σj in your optimization ( i guess you have defined it as a scalar somewhere?). You need to add σj to the params argument so that is also optimized for in the minimization.

Additionally, the sums do not use the correct indices: n Chen 2007 it says: ∑1 (region 1) represents the negative nonzero measured returns, ∑2 (region 2) represents the positive nonzero measured returns, and ∑0 (region 0) represents the zero measured returns. However, you use rjt for all terms.

You need to subset rjt (and all other variables in the sum that depend on t) according to the indices in the sum before each term. In the sums with ∑1, you need to retrieve the indices j for negative nonzero return regions, for ∑2 the indices j for the positive nonzero measured return regions and so on.

You can try to subset the terms as follows:

term_0 = np.sum(np.log(1/(2*np.pi*sigma**2)**0.5))

term_1 = np.sum( (1/ (2*sigma**2) ) * (np.sum(rjt[rjt < 0] +alpha_1-beta_1*djt[rjt < 0]*rft[rjt < 0]) - beta_2*djt[rjt < 0]*index[rjt < 0] )**2)

term_2 = np.sum(np.log(1/(2*np.pi*sigma**2)**0.5))

term_3 = np.sum( (1/ (2*sigma**2) ) * (np.sum(rjt[rjt > 0]+alpha_2-beta_1*djt[rjt > 0]*rft[rjt > 0]) - beta_2 *djt[rjt > 0]*index )**2)

term_4 = np.sum(
          np.log(
            norm.cdf( (alpha_2 - beta_1*djt[rjt == 0]*rft[rjt == 0] - beta_2*djt[rjt == 0]*index[rjt == 0])/sigma, loc=0, scale=1)
          - norm.cdf( (alpha_1 - beta_1*djt[rjt == 0]*rft[rjt == 0] - beta_2*djt[rjt == 0]*index[rjt == 0])/sigma, loc=0, scale=1)
          ))