基于遗传算法的KMV模型最优违约点的确定

作者: hwang_zhic | 来源:发表于2019-04-11 10:34 被阅读0次

本文描述的是基于经典KMV模型企业信用风险评估，使用遗传算法确定适合不同市场的最有违约点。

1.KMV模型简介：主要说明KMV模型是什么，解决什么问题。
2.遗传算法：主要说明遗传算法的思想，内容和使用步骤。
3. 基于遗传算法的KMV模型，主要解释该模型解决什么问题遗传算法的描述和模型的流程
4. GA-KMV模型代码实现,主要介绍Geatpy库内容，数据特征的获取，具体代码说明和结果分析
5. 总结，总结一下该模型。

1. KMV模型

1.1 KMV模型简介

KMV模型是美国旧金山市KMV公司于1997年建立的用来估计借款企业违约概率的方法。该模型认为，贷款的信用风险是在给定负债的情况下由债务人的资产市场价值决定的。但资产并没有真实地在市场交易，资产的市场价值不能直接观测到。为此，模型将银行的贷款问题倒转一个角度，从借款企业所有者的角度考虑贷款归还的问题。在债务到期日，如果公司资产的市场价值高于公司债务值（违约点），则公司股权价值为公司资产市场价值与债务值之间的差额；如果此时公司资产价值低于公司债务值，则公司变卖所有资产用以偿还债务，股权价值变为零。

1.2 KMV模型的运用步骤

首先，它利用Black-Scholes期权定价公式，根据企业资产的市场价值、资产价值的波动性、到期时间、无风险借贷利率及负债的账面价值估计出企业股权的市场价值及其波动性。

其次根据公司的负债计算出公司的违约实施点 (default exercise point，为企业1年以下短期债务的价值加上未清偿长期债务账面价值的一半)，计算借款人的违约距离。

最后，根据企业的违约距离与预期违约率(EDF) 之间的对应关系，求出企业的预期违约率。

KMV模型计算和参数具体可以参考：https://www.jianshu.com/p/ee33c4375a8e

2. 遗传算法

2.1 主要思想

遗传算法是根据达尔文的“适者生存，优胜劣汰”的思想来找到最优解，其特点是所找到的解是全局最优解。

2.1 主要名词

个体（染色体）：一个染色体代表一个具体问题的一个解，一个染色体包含若干基因。

基因：一个基因代表具体问题解的一个决策变量。种群：多个个体（染色体）构成一个种群。即一个问题的多组解构成了解的种群。

目的就是让种群中”优胜劣汰“，最终只剩下一个最优解。接下来介绍最基本遗传算法，只用了选择，交叉，变异三种遗传算子。

2.2 主要步骤

1.种群初始化：我们需要首先通过随机生成的方式来创造一个种群，一般该种群的数量为100~500，这里我们采用二进制将一个染色体(解)编码为基因型。随后用进制转化，将二进制的基因型转化成十进制的表现型。
2.适应度计算(种群评估)：这里我们直接将目标函数值作为个体的适应度。
3.选择(复制)操作：根据种群中个体的适应度大小，通过轮盘赌等方式将适应度高的个体从当前种群中选择出来。其中轮盘赌即是与适应度成正比的概率来确定各个个体遗传到下一代群体中的数量。
具体步骤如下：
(1)首先计算出所有个体的适应度总和Σfi。
(2)其次计算出每个个体的相对适应度大小fi/Σfi，类似于softmax。
(3)再产生一个0到1之间的随机数，依据随机数出现在上述哪个概率区域内来确定各个个体被选中的次数。
4.交叉(交配)运算：该步骤是遗传算法中产生新的个体的主要操作过程，它用一定的交配概率阈值(pc，一般是0.4到0.99)来控制是否采取单点交叉，多点交叉等方式生成新的交叉个体。
具体步骤如下：
(1)先对群体随机配对。
(2)再随机设定交叉点的位置。
(3)再互换配对染色体间的部分基因。
5.变异运算：该步骤是产生新的个体的另一种操作。一般先随机产生变异点，再根据变异概率阈值(pm,一般是0.0001到0.1)将变异点的原有基因取反。
6.终止判断：如果满足条件(迭代次数，一般是200~500)则终止算法，否则返回step2。

3. 基于遗传算法的KMV模型

3.1 问题描述

基于遗传算法的KMV模型（GA-KMV）模型的基本思路是以我国上市公司为样本、通过最小化样本公司类型判定的失误率来求解适用于我国上市公司的最优违约点。
样本公司可以分为 ST（包含 ST 和 ST*公司）和非 ST公司两类。当上市公司连续两年亏损时就会被证券交易所 ST 处理，连续三年亏损就会被 ST*处理，因此 ST 公司具有较高的违约风险，而非 ST 公司财务状况较好、违约风险相对较小。另一方面，可以通过违约距离来判断企业的违约风险。理论上当预期企业资产大于违约点即违约距离大于零时，企业不会违约，反之则会发生违约。将违约距离大于零的公司判定为 H（Healthy）类公司，违约距离小于零的公司判定为D（Distressed）类公司，则 H 类公司的违约风险应该较小，D 类公司的违约风险应该较大。根据以上两种分类结果，可以观察到以下四种情况：

①ST公司被判定为D类公司；
②ST公司被判定为H类公司；
③非ST公司被判定为H类公司；
④非ST公司被判定为D类公司。

在运用违约距离判定公司类型时，以上四种情况都有可能发生。当第一或第三种情况发生时，模型预测正确；而第二或第四种情况发生时，模型预测失误。GA-KMV模型就是要寻求这样一个违约点，使得通过其判定样本公司类型时，第二和第四种情况发生概率的平均值即平均误判率最小。设样本公司总数量为N，ST和非ST公司数均为N/2。将违约点设定为短、长期债务的线性组合即：

图3.1.1 对于每一组(α,β)，根据违约距离公式：

图3.1.2

可以求得各公司的违约距离，再根据以上判定规则就可以判定各公司类型，这样对于每个违约点样本公司都有一个分类。假定 N / 2 家ST公司中被判定为D类公司的数量为m，N / 2 家非ST公司中被判定为H类公司的数量n，则第二和第四种情况发生的概率分别为

求解最优违约点的问题就转化为求解(α,β)使得平均误判率最小：

3.2 算法描述

对于以上问题，本文运用遗传算法进行求解。以α、β为自变量，将平均误判率设定为适应度函数，通过选择、交叉、突变等操作，使解逐代优化，最终得到一个适用于我国上市公司的最优违约点。

（1）编码表示

编码是把问题的可行解从其解空间转换为算法所能处理的搜索空间的转换方法。对一些多极值函数和多维、高精度要求的连续函数优化问题，使用二进制编码难以得到最优解。浮点数编码以决策变量个数为编码长度，以一个浮点数作为基因值，它克服了二进制编码方法的缺点，能够在较大空间进行搜索，适合于精度要求较高的遗传算法。因此本文采用浮点数编码方式，对短、长期债务系数α、β进行编码。

（2）创建初始种群

编码表示之后，需要初始化一定数目的个体作为初始种群。本文将种群规模设定为50，以KMV模型中短、长期债务系数组(1,0.5)作为初始值，将初始种群范围设定为[0,10]，根据均匀分布在初始范围内创建初始种群。

（3）适应性评价

适应度值是评价个体优劣的标准。在本文中，平均误判率是评价一个违约点优劣的标准且平均误判率的取值范围为[0,1]，因此可以将适应度函数设定为

3.14

（4）遗传算子

选择算子是从父代种群中以一定的概率选择出个体到新种群中的操作。锦标赛选择策略每次从种群中选择一定数量的个体，再取出其中的最优个体进入下一代，重复进行直至达到种群规模。锦标赛选择策略可以使种群多样性较为丰富，同时保证了被选个体较优。

交叉算子是将种群中的个体两两分组，按一定概率和方式交换部分基因的操作。交叉操作能够组合出新的个体，在串空间进行有效搜索，同时降低对种群有效模式的破坏概率。本文采用离散交叉方式，在个体间交换变量的值，子个体的每个变量按等概率随机挑选父个体。

变异算子是对个体编码串中以变异概率和随机指定的某一位或某几位基因值做取反运算的操作。变异操作可以改善遗传算法的局部搜索能力，丰富种群多样性。本文选用均匀突变策略，每一个实数元素以相同的概率在域内变动。同时为了防止早熟收敛问题，

（5）终止条件本文的终止条件设定为：迭代代数为200

3.3 GA-KMV 模型流程

在求解最优违约点之前，首先需计算样本公司资产价值AV和其波动率Aσ。然后将它们作为模型的输入参数。GA-KMV模型具体执行过程如下（如图3.3所示）：
（1）编码：对违约点的短、长期债务系数α、β进行编码。
（2）创建初始种群，群体中每个个体就是一组(),αβ值，代表了一个可能的违约点。
（3）计算违约距离DD：对于每个违约点，运用式（3.11）计算每个样本公司违约距离。
（4）判定公司类型：对于每个违约点，分别对ST和非ST样本公司类型进行判定，违约距离大于0则该公司被判定为H类公司，否则为D类公司；统计ST公司中被判定为D类公司的数量m和非ST公司中被判定为H类公司的数量n。
（5）计算每个个体适应度值，对每一组(),αβ值即每个违约点，计算第二和第四情况发生频率的平均值即误判率（1-（m+n）/N）将误判率作为该违约点的适应度值。
（6）进行选择、交叉、变异操作。
（7）检验是否达到终止条件，如果达到终止条件，停止算法，得到最优解即短、长期债务的最优系数，否则返回步骤（3）。

图3.3

4. GA-KMV模型代码实现

4.1 Python遗传和进化算法框架——Geatpy

4.1.1 Geatpy安装

这里为了简单代码的实现，使用了Python遗传和进化算法框架——Geatpy
这里安装很简单，只用打开pycham里的settings，步骤如下就安装完成。

4.1.2 Geatpy主要数据类型

1）种群染色体的数据结构：

Geatpy 中，种群染色体是一个二维矩阵，简称“种群矩阵”。一般所说的“种群”是特指种群染色体矩阵。

种群矩阵一般用Chrom 命名，是一个是numpy 的array 类型的，每一行对应一条染色体，同时也对应着一个个体。染色体的每个元素是染色体上的基因。

一般把种群的规模(即种群的个体数) 用Nind 命名；把种群个体的染色体长度用Lind 命名。其结构如下图所示：

2）种群表现型的数据结构：

种群表现型的数据结构跟种群染色体基本一致，也是numpy 的array 类型。我们一般用Phen 来命名(也可以命名为variable)。它是种群矩阵Chrom 经过解码操作后得到的基因表现型矩阵，每一行对应一个个体，每行中每个元素都代表着一个控制变量，并用Nvar 表示控制变量的个数。如下图：

3）个体适应度的数据结构：

每一行对应种群每个个体的适应度值。

4）个体可行性的数据结构：

Geatpy 采用列向量来存储种群个体可行性。一般命名为LegV ，它同样是numpy 的array 类型，每一行对应种群矩阵的每一个个体，表示对应的个体是否是可行解，0 表示非可行解，1 表示可行解。因此它拥有与Chrom 相同的行数。

5）区域描述器的数据结构：
区域描述器是用来描述种群染色体的特征，比如染色体中基因所表达的控制变量的范围、是否包含范围的边界、采用什么编码方式，是否使用对数刻度等等。
lens 包含染色体的每个子染色体的长度。sum(lens) 等于染色体长度。

lb 和ub 分别代表每个变量的上界和下界

codes 指明染色体子串用的是标准二进制编码还是格雷编码。codes[i] = 0 表示第i个变量使用的是标准二进制编码；codes[i] = 1 表示使用格雷编码。

scales 指明每个子串用的是算术刻度还是对数刻度。scales[i] = 0 为算术刻度，scales[i] = 1 为对数刻度。对数刻度可以用于变量的范围较大而且不确定的情况，对于大范围的参数边界，对数刻度让搜索可用较少的位数，从而减少了遗传算法的计算量。

lbin 和ubin 指明了变量是否包含其范围的边界。0 表示不包含边界；1 表示包含边界。

4.2 模型数据特征

这次的所有上市公司数据均来自锐思数据库并进行相关整理得到。
首先介绍模型所需和得出特征或数据，具体如下：
VE：当期的股权价值
VA：资产价值
σE ：股权收益率的波动率
σA：资产收益率的波动率
D：负债的账面价值
T：债务期限
r：无风险利率

下面介绍以上所需特征的获取方法

4.2.1 直接获取

预测周期T：通常设置为1年

年收益率和无风险收益率r可以直接得出，如下图，可以直接读出这两个特征的数值：

图4.2.1 （ST上市公司股价数据）

企业总债务D直接从图中读取（负债总计）：

图4.2.2 （ST上市公司负债数据）

4.2.2 计算获取

违约触发点DPT：DPT = aSDT + bLDT；其中a，b为遗传算法中的种群中的每个不同个体，SDT为短期借款、LDT为长期借款，都可从表中直接读取，如图4.2.2，分别对应流动负债和非流动负债

企业股权价值VE：
中国证券市场发展历史较为特殊，上市公司股票被人为分割为上市流通股票和暂不上市流通股票两种。在计算上市公司股权市场价值时需要考虑以什么样的价格来计算非流通股市场价值/由于非流通股没有市场交易价格，因此如何给非流通股定价是一件困难的事情。参考上市公司股票全流通研究中非流通股定价，以每股净资产计算非流通股的价格。公式如下：
流通股市场价值=12月份周平均收盘价格×流通股股数
非流通股市场价值=每股净资产×非流通股股数
上市公司股权市场价值=流通股市场价值+非流通股市场价值

以上数据”收盘价格“、“流通股股数”、“每股净资产”和“非流通股股数”都可从表中读取，如图4.2.1所示

4.3 GA-KMV模型代码

该部分直接贴出所以代码，并逐个分析

1.首先是import库，问题不大

2.接下来是数据的读取，首先设置总公司数和计算st和非st的公司数，然后st和st*的公司数，由于st公司数据跟st*公司数据分开，所以要这样设置好变量方便使用。
然后是读取excel表，表格总共有6个因此要逐一读取和提取其中的需要的列的数据。
提取完后还要继续细分给公司并进行相关计算，最终得到了ST、ST*和非ST公司的“无风险列表”、“股权收益率波动率列表”、“股权市场价值列表”、“短期负债列表”、“长期负债列表”、“负债合计列表”

3.KMVOptSearch()计算资产价值Va和资产收益率波动率σ。返回这两个值的列表

4.calculate_DD计算DD距离并返归m，n值。i为a，b的当前个体下标；counter为当前迭代公司的次数；IS_ST_bool为判断是否为ST公司，当判断为1时，计算m值，否则计算n值。其他为相关数据。

5.aim()定义目标适应函数，输出个体适应度矩阵和可行矩阵；这是模板中必须调用目标函数。里面分别计算了每个个体的a，b的适应值(或者说误判率)并添加到diff的适应值矩阵中，最终返回。

6.main()函数，因为是两个变量a和b，因此设置了两个列表x1和x2描述这两个变量的变化范围并设置b1,b2描述是否包含边界(0为不包含，1为包含)，跟着生成相关的矩阵。然后生成区域描述器FieldDR，参数为上述的范围和边界矩阵。
最后调用sga_new_real_templet()（改进的单目标进化算法模板（实值编码））参数分别为：
（1）AIM_M - 目标函数的地址；
（2）AIM_F : str - 目标函数名；
（3）AFieldDR : array - 实际值种群区域描述器；
（4）problem : str - 表明是整数问题还是实数问题，'I'表示是整数问题，'R'表示是实数问题；
（5）maxormin int - 最小最大化标记，1表示目标函数最小化；-1表示目标函数最大化；
（6）MAXGEN : int - 最大遗传代数；
（7）NIND : int - 种群规模，即种群中包含多少个个体；
（8） SUBPOP : int - 子种群数量，即对一个种群划分多少个子种群；
（9）selectStyle : str - 指代所采用的低级选择算子的名称；
（10）recombinStyle: str - 指代所采用的低级重组算子的名称，如'xovsp'(单点交叉)；
（11）recopt : float - 交叉概率；
（12）pm : float - 重组概率，即为变异概率；
（13）distribute : bool - 是否增强种群的分布性（可能会造成收敛慢）；
（14）drawing : int - (可选参数)，0表示不绘图，1表示绘制最终结果图，2表示绘制动画。默认drawing为1；
其他未说明的为用不到的。
该程序里的sga_new_real_templet()直接实现了遗传算法中的编码表示、种群生成、适应性评价、遗传算子和迭代终止，只要相关函数和参数设置好了就能直接调用了。

完整代码：

# -*- coding: utf-8 -*-
"""
    功能：基于遗传算法确定最优违约触发点
    版本：v1.0
"""

from scipy import stats
from scipy.optimize import fsolve
import numpy as np
import pandas as pd
import geatpy as ga

# 公司总样本数 / 2
N = 64
# ST 和非 ST 公司数均为 N / 2
num_com = int((N / 2))
# ST 和 ST*公司数均为
count = int((num_com / 2))

# 读取ST上市公司数据
ST_data_xls_1 = pd.read_excel("ST公司数据.xls")
ST_data_xls_2 = pd.read_excel("ST公司数据负债.xls")
# 读取ST*上市公司数据
ST_Plus_data_xls_1 = pd.read_excel("ST+公司数据.xls")
ST_Plus_data_xls_2 = pd.read_excel("ST+公司数据负债.xls")
# 读取非ST上市公司数据
NST_data_xls_1 = pd.read_excel("非ST公司数据.xls")
NST_data_xls_2 = pd.read_excel("非ST公司数据负债.xls")

# 读取ST上市公司数据
ST_data_1 = ST_data_xls_1.ix[: count,
            ['收盘价_Clpr', '流通股_Trdshr', '已上市流通股_Lsttrdshr', '年收益率_Yrret', '年无风险收益率_Yrrfret', '每股净资产(元/股)_NAPS']].values
ST_data_2 = ST_data_xls_2.ix[: count,
            ['流动负债合计(元)_TotCurLia', '非流动负债合计(元)_TotNCurLia', '负债合计(元)_TotLiab']].values

# 读取ST*上市公司数据
ST_Plus_data_1 = ST_Plus_data_xls_1.ix[: count,
                 ['收盘价_Clpr', '流通股_Trdshr', '已上市流通股_Lsttrdshr', '年收益率_Yrret', '年无风险收益率_Yrrfret',
                  '每股净资产(元/股)_NAPS']].values
ST_Plus_data_2 = ST_Plus_data_xls_2.ix[: count,
                 ['流动负债合计(元)_TotCurLia', '非流动负债合计(元)_TotNCurLia', '负债合计(元)_TotLiab']].values

# 读取非ST上市公司数据
NST_data_1 = NST_data_xls_1.ix[: num_com,
             ['收盘价_Clpr', '流通股_Trdshr', '已上市流通股_Lsttrdshr', '年收益率_Yrret', '年无风险收益率_Yrrfret', '每股净资产(元/股)_NAPS']].values
NST_data_2 = NST_data_xls_2.ix[: num_com,
             ['流动负债合计(元)_TotCurLia', '非流动负债合计(元)_TotNCurLia', '负债合计(元)_TotLiab']].values

# ST公司数据
# 收盘价_Clpr
ST_Clpr_list = ST_data_1[:, 0]
# 流通股_Trdshr
ST_Trdshr_list = ST_data_1[:, 1]
# 已上市流通股_Lsttrdshr
ST_Lsttrdshr_list = ST_data_1[:, 2]
# 每股净资产
ST_NAPS_list = ST_data_1[:, 5]

# 无风险列表
ST_r_list = ST_data_1[:, 4]
# 股权收益率波动率列表
ST_PriceTheta_list = ST_data_1[:, 3]
# 股权市场价值列表
ST_E_list = ST_Clpr_list * ST_Lsttrdshr_list + (ST_Trdshr_list - ST_Lsttrdshr_list) * ST_NAPS_list
# Short-term Debt (in 10,000)
ST_SD_list = ST_data_2[:, 0]
# print(ST_SD_list)
# Long-term Debt (in 10,000)
ST_LD_list = ST_data_2[:, 1]
# print(ST_LD_list)
# 负债合计
ST_D_list = ST_data_2[:, 2]

# ST*公司数据
# 收盘价_Clpr
ST_Plus_Clpr_list = ST_Plus_data_1[:, 0]
# 流通股_Trdshr
ST_Plus_Trdshr_list = ST_Plus_data_1[:, 1]
# 已上市流通股_Lsttrdshr
ST_Plus_Lsttrdshr_list = ST_Plus_data_1[:, 2]
# 每股净资产
ST_Plus_NAPS_list = ST_Plus_data_1[:, 5]

# 无风险列表
ST_Plus_r_list = ST_Plus_data_1[:, 4]
# 股权收益率波动率列表
ST_Plus_PriceTheta_list = ST_data_1[:, 3]
# 股权市场价值列表
ST_Plus_E_list = ST_Plus_Clpr_list * ST_Plus_Lsttrdshr_list + (
        ST_Plus_Trdshr_list - ST_Plus_Lsttrdshr_list) * ST_Plus_NAPS_list
# Short-term Debt (in 10,000)
ST_Plus_SD_list = ST_Plus_data_2[:, 0]
# print(ST_SD_list)
# Long-term Debt (in 10,000)
ST_Plus_LD_list = ST_Plus_data_2[:, 1]
# print(ST_LD_list)
# 负债合计
ST_Plus_D_list = ST_Plus_data_2[:, 2]

# 非ST公司数据
# 收盘价_Clpr
NST_Clpr_list = NST_data_1[:, 0]
# 流通股_Trdshr
NST_Trdshr_list = NST_data_1[:, 1]
# 已上市流通股_Lsttrdshr
NST_Lsttrdshr_list = NST_data_1[:, 2]
# 每股净资产
NST_NAPS_list = NST_data_1[:, 5]

# 无风险列表
NST_r_list = NST_data_1[:, 4]
# 股权收益率波动率列表
NST_PriceTheta_list = NST_data_1[:, 3]
# 股权市场价值列表
NST_E_list = NST_Clpr_list * NST_Lsttrdshr_list + (NST_Trdshr_list - NST_Lsttrdshr_list) * NST_NAPS_list
# Short-term Debt (in 10,000)
NST_SD_list = NST_data_2[:, 0]
# print(ST_SD_list)
# Long-term Debt (in 10,000)
NST_LD_list = NST_data_2[:, 1]
# 负债合计
NST_D_list = NST_data_2[:, 2]


def KMVOptSearch(E, D, r, T, EquityTheta):
    """
       KMV模型计算
    """
    EtoD = float(E) / float(D)

    def KMVfun(x, EtoD, r, T, EquityTheta):
        """
            解d1,d2方程组
        """
        EtoD, r, T, EquityTheta = float(EtoD), float(r), float(T), float(EquityTheta)
        d1 = (np.log(x[0] * EtoD) + (r + 0.5 * x[1] ** 2) * T) / (x[1] * np.sqrt(T))
        d2 = d1 - x[1] * np.sqrt(T)

        return [
            x[0] * stats.norm.cdf(d1, 0.0, 1.0) - np.exp(-r * T) * stats.norm.cdf(d2, 0.0, 1.0) / EtoD - 1,
            stats.norm.cdf(d1, 0.0, 1.0) * x[0] * x[1] - EquityTheta
        ]

    VaThetaX = fsolve(KMVfun, [1, 0.1], args=(EtoD, r, T, EquityTheta))
    Va = VaThetaX[0] * E
    AssetTheta = VaThetaX[1]
    return Va, AssetTheta


def calculate_DD(i, counter, r_list, SD_list, LD_list, a, b, D_list, PriceTheta_list, E_list, IS_ST_bool):
    """
        计算DD距离并返归m，n值
    """
    m = 0
    n = 0
    for time in range(counter):
        # 利用KMV模型计算
        r = r_list[time]  # 无风险利率
        T = 1  # 期权到期时间，以年为单位
        SD = SD_list[time]  # SD为长期负债 （单位：元）
        LD = LD_list[time]  # LD短期负债 （单位：元）
        DP = a[i] * SD + b[i] * LD  # 违约触发点DPT # DP=（TD-LD） +0.5LD （7） 其中：TD为总负债；LD为长期负债。 （单位：元）
        D = D_list[time]  # D=TD-LD+0.5∗LD  债券的面值是D （单位：元）
        PriceTheta = PriceTheta_list[time]  # 交易波动率，可以是日月年的数据 （单位：%）
        EquityTheta = PriceTheta  # 股权收益率的波动率σE
        E = E_list[time]  # 公司股票的市场价值 （单位：元）
        Va, AssetTheta = KMVOptSearch(E, D, r, T, EquityTheta)  # 资产价值VA（单位：元）资产收益率波动率σA
        AssetTheta = AssetTheta  # 资产收益率波动率σA （单位：%）
        DD = (Va - DP) / (Va * AssetTheta)  # 违约距离DD
        EDF = stats.norm.cdf(-DD)  # 预期违约概率EDF（单位：%）
        # print("ST_DD = {}".format(DD))
        if IS_ST_bool == 1:
            if DD >= 0:
                m += 1
        elif IS_ST_bool == 0:
            if DD < 0:
                n += 1

    return [m, n]


# Attention: you had better put the aim-function in another file.
def aim(Phen, LegV):  # define the aim function
    """
        定义目标函数，输出个体适应度矩阵和可行矩阵
    """
    global counter
    counter += 1
    print("第{}次迭代开始：".format(counter))

    for i in range(Phen.shape[0]):
        n = 0
        m = 0
        a = Phen[:, 0]
        b = Phen[:, 1]

        # 计算ST公司的m，n值
        mn_list = calculate_DD(i, count, ST_r_list, ST_SD_list, ST_LD_list, a, b, ST_D_list, ST_PriceTheta_list,
                               ST_E_list,
                               1)

        m += mn_list[0]
        n += mn_list[1]

        # 计算*ST公司的m，n值
        mn_list = calculate_DD(i, count, ST_Plus_r_list, ST_Plus_SD_list, ST_Plus_LD_list, a, b, ST_Plus_D_list,
                               ST_Plus_PriceTheta_list, ST_Plus_E_list, 1)
        m += mn_list[0]
        n += mn_list[1]

        # 计算非ST公司的m，n值
        mn_list = calculate_DD(i, N, NST_r_list, NST_SD_list, NST_LD_list, a, b, NST_D_list, NST_PriceTheta_list,
                               NST_E_list, 0)
        m += mn_list[0]
        n += mn_list[1]

        print("i:", i + 1)
        print("a = {}, b = {}".format(a[i], b[i]))
        print("m = {}, n = {}".format(m, n))
        Mis_rate = 1 - ((m + n) / N)
        print("Mis_rate：", Mis_rate)
        if i == 0:
            diff = np.array([[Mis_rate]])
        else:
            diff = np.append(diff, [[Mis_rate]], axis=0)
        print("diff.shape:", diff.shape)
        print()

    print("[diff, LegV][diff, LegV][diff, LegV][diff, LegV]:", [diff, LegV])
    return [diff, LegV]


global counter
counter = 0


def main():
    """
        主函数：主要设定遗传算法变量和调用
    """

    AIM_M = __import__('aimfuc')
    AIM_M = __import__('main2')  # get the handle of aim-function
    # variables setting

    # 变量设置
    x1 = [0, 10]  # 自变量1的范围
    x2 = [0, 10]  # 自变量2的范围
    b1 = [0, 1]  # 自变量1是否包含下界
    b2 = [0, 1]  # 自变量2是否包含上界

    ranges = np.vstack([x1, x2]).T  # 生成自变量的范围矩阵
    borders = np.vstack([b1, b2]).T  # 生成自变量的边界矩阵
    # 生成区域描述器
    FieldDR = ga.crtfld(ranges, borders)  # create FieldDR

    # call the GEA algorithm template
    [pop_trace, var_trace, times] = ga.sga_new_real_templet(AIM_M, 'aim', None, None, FieldDR, problem='R', maxormin=1,
                                                            MAXGEN=200, NIND=40, SUBPOP=1, GGAP=0.9,
                                                            selectStyle='etour',
                                                            recombinStyle='xovdprs', recopt=0.9, pm=None,
                                                            distribute=True, drawing=1)

    # output results
    for num in var_trace[np.argmin(pop_trace[:, 1]), :]:
        print(chr(int(num)), end='')


# 若调用该文件则调用main()主函数
if __name__ == "__main__":
    main()

运行结果：

图4.3.1

图4.3.2

4.4 结果分析

图4.3.2展示了一次运行的迭代收敛过程。由图3.3可知，适应值函数在前4代内迅速下降，在第5代开始收敛于0.375，收敛速度很快，最终算法在200代终止。图4.3.2还显示前10代内种群平均适应度值也快速下降，之后一直在0.375上下浮动，算法终止时种群的平均适应度值为0.387，说明在进化初期种群整体适应度迅速得到优化，之后种群平均适应度值比较平稳。
从图4.3.1里得到最优违约点的短、长期负债系数分别为5.289437953188504、10.0。
但是整体的运算速度很慢，一共用了7400多秒。

5.0 总结

本次使用遗传算法优化KMV模型中的违约触发点，得到最优违约点还是有用的。因为传统KMV模型的系数是KMV公司根据他们那边市场的历史数据得出来的，而使用遗传算法可以适当优化并得到适用于中国市场的最优违约点参数设置。

Geatpy库使用参考：https://blog.csdn.net/weixin_37790882/article/details/84034956
https://blog.csdn.net/qq_33353186/article/details/82020507
参考代码、数据文件和文章或论文的码云地址：https://gitee.com/hwang_zc/GA_KMV