Python 中的最大化对数似然估计

Maximizing Log Likelihood Estimation in Python

目前我正在写关于绿色债券流动性的论文。我正在尝试根据 Lesmond 等人计算 LOT 度量。 (1999) 和 Chen 等人。 (2007)。因此,我需要最大化以下对数似然函数

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

in python(附上屏幕截图以提高可读性):

Φ i,j

表示在

评估的每个债券年的累积分布函数
L(a i,j− β j,1D j,t∗ ∆R f,t− β j,2D j,t∗ ∆Index t)/σ j

到目前为止,我尝试应用以下代码,但不幸的是收到了错误的估计:

 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)

我不确定我的代码的实现是否正确。

与您提供的解析表达式相比,您的公式有几个错误:

in term_2 = np.sum(np.log(1/(2*np.pi*sigma)**0.5)) 西格玛应该平方 -> term_2 = np.sum(np.log(1/(2*np.pi*sigma**2)**0.5)).

term_0也是这样,把term_0改成term_0 = np.sum(np.log(1/(2*np.pi*sigma**2)**0.5))

您似乎也在最小化 对数似然。这与您想要做的相反。最小化负 LL 变化

return(LL1)

return(-LL1)

您也没有正确评估 term_4 中的 cdfs。在 Lesmond 1999 中,它指出 Φi,j 是在 (,−1Duration,*Δ−2Duration,*ΔS&P Index)/ 处评估的标准正态分布的 cdf。但是,您编写的程序是评估具有均值 (,−1Duration,∗Δ−2Duration,∗ΔS&P Index)/ 和方差 sigma^2 在某个值 y 的非标准正态分布。要评估正确的 cdf,您需要将 term_4 更改为

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)
            ))

,其中在第一个术语中使用 alpha_2,在第二个术语中使用 alpha_1,如 Chen 2007 中所述。

在 Chen 2007 中,它还指出 α1 j、α2 j、βj 和 σj 是通过最大化似然来求解的。您尚未在优化中包含 σj (我猜您已将其定义为某处的标量?)。您需要将 σj 添加到 params 参数,以便在最小化中也对其进行优化。

此外,总和没有使用正确的索引: n Chen 2007 说:∑1(区域 1)表示测量的负非零值 returns,∑2(区域 2)表示测量的正非零值 returns,∑0(区域 0)表示零测量 returns。但是,您对所有术语都使用 rjt。

您需要根据总和中每一项之前的索引对 rjt(以及总和中依赖于 t 的所有其他变量)进行子集化。在与 ∑1 的总和中,您需要检索负非零 return 区域的索引 j,对于 ∑2,正非零测量 return 区域的索引 j 等等。

您可以尝试按如下方式对条款进行子集化:

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)
          ))