liaozhaorun/NewStock
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
4090b3c5befde9a1e8dcf8e02ccefc7a6b587d14
NewStock/main/train/Classify2.ipynb
liaozhaorun 4090b3c5be refactor(qmt): 优化配置模型和百分比交易逻辑 
- 统一使用配置模型属性访问替代字典下标
- 完善百分比模式买卖的日志记录和错误处理
- 代码格式化和死代码清理
- 更新 notebook 数据及测试脚本
2026-03-03 15:16:37 +08:00

2517 lines
180 KiB
Plaintext
Raw Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"id": "79a7758178bafdd3",
"metadata": {
"ExecuteTime": {
"end_time": "2025-04-03T12:46:06.987506Z",
"start_time": "2025-04-03T12:46:06.259551Z"
},
"jupyter": {
"source_hidden": true
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"/mnt/d/PyProject/NewStock\n"
]
}
],
"source": [
"%load_ext autoreload\n",
"%autoreload 2\n",
"# %load_ext cudf.pandas\n",
"\n",
"import gc\n",
"import os\n",
"import sys\n",
"sys.path.append('/mnt/d/PyProject/NewStock/')\n",
"print(os.getcwd())\n",
"import pandas as pd\n",
"from main.factor.factor import get_rolling_factor, get_simple_factor\n",
"from main.utils.factor import read_industry_data\n",
"from main.utils.factor_processor import calculate_score\n",
"from main.utils.utils import read_and_merge_h5_data, merge_with_industry_data\n",
"\n",
"import warnings\n",
"\n",
"warnings.filterwarnings(\"ignore\")\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "4a481c60",
"metadata": {},
"outputs": [],
"source": [
"# 设置使用核心\n",
"import os\n",
"os.environ[\"MODIN_CPUS\"] = \"4\"\n"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "a79cafb06a7e0e43",
"metadata": {
"ExecuteTime": {
"end_time": "2025-04-03T12:47:00.212859Z",
"start_time": "2025-04-03T12:46:06.998047Z"
},
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"daily data\n",
"daily basic\n",
"inner merge on ['ts_code', 'trade_date']\n",
"stk limit\n",
"left merge on ['ts_code', 'trade_date']\n",
"money flow\n",
"left merge on ['ts_code', 'trade_date']\n",
"cyq perf\n",
"left merge on ['ts_code', 'trade_date']\n",
"<class 'pandas.core.frame.DataFrame'>\n",
"RangeIndex: 9267817 entries, 0 to 9267816\n",
"Data columns (total 33 columns):\n",
" # Column Dtype \n",
"--- ------ ----- \n",
" 0 ts_code object \n",
" 1 trade_date datetime64[ns]\n",
" 2 open float64 \n",
" 3 close float64 \n",
" 4 high float64 \n",
" 5 low float64 \n",
" 6 vol float64 \n",
" 7 amount float64 \n",
" 8 pct_chg float64 \n",
" 9 turnover_rate float64 \n",
" 10 pe_ttm float64 \n",
" 11 circ_mv float64 \n",
" 12 total_mv float64 \n",
" 13 volume_ratio float64 \n",
" 14 is_st bool \n",
" 15 up_limit float64 \n",
" 16 down_limit float64 \n",
" 17 buy_sm_vol float64 \n",
" 18 sell_sm_vol float64 \n",
" 19 buy_lg_vol float64 \n",
" 20 sell_lg_vol float64 \n",
" 21 buy_elg_vol float64 \n",
" 22 sell_elg_vol float64 \n",
" 23 net_mf_vol float64 \n",
" 24 his_low float64 \n",
" 25 his_high float64 \n",
" 26 cost_5pct float64 \n",
" 27 cost_15pct float64 \n",
" 28 cost_50pct float64 \n",
" 29 cost_85pct float64 \n",
" 30 cost_95pct float64 \n",
" 31 weight_avg float64 \n",
" 32 winner_rate float64 \n",
"dtypes: bool(1), datetime64[ns](1), float64(30), object(1)\n",
"memory usage: 2.2+ GB\n",
"None\n"
]
}
],
"source": [
"from main.utils.utils import read_and_merge_h5_data\n",
"\n",
"print('daily data')\n",
"df = read_and_merge_h5_data('/mnt/d/PyProject/NewStock/data/daily_data.h5', key='daily_data',\n",
" columns=['ts_code', 'trade_date', 'open', 'close', 'high', 'low', 'vol', 'amount', 'pct_chg'],\n",
" df=None)\n",
"\n",
"print('daily basic')\n",
"df = read_and_merge_h5_data('/mnt/d/PyProject/NewStock/data/daily_basic.h5', key='daily_basic',\n",
" columns=['ts_code', 'trade_date', 'turnover_rate', 'pe_ttm', 'circ_mv', 'total_mv', 'volume_ratio',\n",
" 'is_st'], df=df, join='inner')\n",
"\n",
"print('stk limit')\n",
"df = read_and_merge_h5_data('/mnt/d/PyProject/NewStock/data/stk_limit.h5', key='stk_limit',\n",
" columns=['ts_code', 'trade_date', 'pre_close', 'up_limit', 'down_limit'],\n",
" df=df)\n",
"print('money flow')\n",
"df = read_and_merge_h5_data('/mnt/d/PyProject/NewStock/data/money_flow.h5', key='money_flow',\n",
" columns=['ts_code', 'trade_date', 'buy_sm_vol', 'sell_sm_vol', 'buy_lg_vol', 'sell_lg_vol',\n",
" 'buy_elg_vol', 'sell_elg_vol', 'net_mf_vol'],\n",
" df=df)\n",
"print('cyq perf')\n",
"df = read_and_merge_h5_data('/mnt/d/PyProject/NewStock/data/cyq_perf.h5', key='cyq_perf',\n",
" columns=['ts_code', 'trade_date', 'his_low', 'his_high', 'cost_5pct', 'cost_15pct',\n",
" 'cost_50pct',\n",
" 'cost_85pct', 'cost_95pct', 'weight_avg', 'winner_rate'],\n",
" df=df)\n",
"print(df.info())"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "cac01788dac10678",
"metadata": {
"ExecuteTime": {
"end_time": "2025-04-03T12:47:10.527104Z",
"start_time": "2025-04-03T12:47:00.488715Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"industry\n"
]
}
],
"source": [
"print('industry')\n",
"industry_df = read_and_merge_h5_data('/mnt/d/PyProject/NewStock/data/industry_data.h5', key='industry_data',\n",
" columns=['ts_code', 'l2_code', 'in_date'],\n",
" df=None, on=['ts_code'], join='left')\n",
"\n",
"\n",
"def merge_with_industry_data(df, industry_df):\n",
" # 确保日期字段是 datetime 类型\n",
" df['trade_date'] = pd.to_datetime(df['trade_date'])\n",
" industry_df['in_date'] = pd.to_datetime(industry_df['in_date'])\n",
"\n",
" # 对 industry_df 按 ts_code 和 in_date 排序\n",
" industry_df_sorted = industry_df.sort_values(['in_date', 'ts_code'])\n",
"\n",
" # 对原始 df 按 ts_code 和 trade_date 排序\n",
" df_sorted = df.sort_values(['trade_date', 'ts_code'])\n",
"\n",
" # 使用 merge_asof 进行向后合并\n",
" merged = pd.merge_asof(\n",
" df_sorted,\n",
" industry_df_sorted,\n",
" by='ts_code', # 按 ts_code 分组\n",
" left_on='trade_date',\n",
" right_on='in_date',\n",
" direction='backward'\n",
" )\n",
"\n",
" # 获取每个 ts_code 的最早 in_date 记录\n",
" min_in_date_per_ts = (industry_df_sorted\n",
" .groupby('ts_code')\n",
" .first()\n",
" .reset_index()[['ts_code', 'l2_code']])\n",
"\n",
" # 填充未匹配到的记录（trade_date 早于所有 in_date 的情况）\n",
" merged['l2_code'] = merged['l2_code'].fillna(\n",
" merged['ts_code'].map(min_in_date_per_ts.set_index('ts_code')['l2_code'])\n",
" )\n",
"\n",
" # 保留需要的列并重置索引\n",
" result = merged.reset_index(drop=True)\n",
" return result\n",
"\n",
"\n",
"# 使用示例\n",
"df = merge_with_industry_data(df, industry_df)\n",
"# print(mdf[mdf['ts_code'] == '600751.SH'][['ts_code', 'trade_date', 'l2_code']])"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "c4e9e1d31da6dba6",
"metadata": {
"ExecuteTime": {
"end_time": "2025-04-03T12:47:10.719252Z",
"start_time": "2025-04-03T12:47:10.541247Z"
},
"jupyter": {
"source_hidden": true
}
},
"outputs": [],
"source": [
"from main.factor.factor import *\n",
"\n",
"def calculate_indicators(df):\n",
" \"\"\"\n",
" 计算四个指标：当日涨跌幅、5日移动平均、RSI、MACD。\n",
" \"\"\"\n",
" df = df.sort_values('trade_date')\n",
" df['daily_return'] = (df['close'] - df['pre_close']) / df['pre_close'] * 100\n",
" # df['5_day_ma'] = df['close'].rolling(window=5).mean()\n",
" delta = df['close'].diff()\n",
" gain = delta.where(delta > 0, 0)\n",
" loss = -delta.where(delta < 0, 0)\n",
" avg_gain = gain.rolling(window=14).mean()\n",
" avg_loss = loss.rolling(window=14).mean()\n",
" rs = avg_gain / avg_loss\n",
" df['RSI'] = 100 - (100 / (1 + rs))\n",
"\n",
" # 计算MACD\n",
" ema12 = df['close'].ewm(span=12, adjust=False).mean()\n",
" ema26 = df['close'].ewm(span=26, adjust=False).mean()\n",
" df['MACD'] = ema12 - ema26\n",
" df['Signal_line'] = df['MACD'].ewm(span=9, adjust=False).mean()\n",
" df['MACD_hist'] = df['MACD'] - df['Signal_line']\n",
"\n",
" # 4. 情绪因子1：市场上涨比例（Up Ratio）\n",
" df['up_ratio'] = df['daily_return'].apply(lambda x: 1 if x > 0 else 0)\n",
" df['up_ratio_20d'] = df['up_ratio'].rolling(window=20).mean() # 过去20天上涨比例\n",
"\n",
" # 5. 情绪因子2：成交量变化率（Volume Change Rate）\n",
" df['volume_mean'] = df['vol'].rolling(window=20).mean() # 过去20天的平均成交量\n",
" df['volume_change_rate'] = (df['vol'] - df['volume_mean']) / df['volume_mean'] * 100 # 成交量变化率\n",
"\n",
" # 6. 情绪因子3：波动率（Volatility）\n",
" df['volatility'] = df['daily_return'].rolling(window=20).std() # 过去20天的日收益率标准差\n",
"\n",
" # 7. 情绪因子4：成交额变化率（Amount Change Rate）\n",
" df['amount_mean'] = df['amount'].rolling(window=20).mean() # 过去20天的平均成交额\n",
" df['amount_change_rate'] = (df['amount'] - df['amount_mean']) / df['amount_mean'] * 100 # 成交额变化率\n",
"\n",
" # df = sentiment_panic_greed_index(df)\n",
" # df = sentiment_market_breadth_proxy(df)\n",
" # df = sentiment_reversal_indicator(df)\n",
"\n",
" return df\n",
"\n",
"\n",
"def generate_index_indicators(h5_filename):\n",
" df = pd.read_hdf(h5_filename, key='index_data')\n",
" df['trade_date'] = pd.to_datetime(df['trade_date'], format='%Y%m%d')\n",
" df = df.sort_values('trade_date')\n",
"\n",
" # 计算每个ts_code的相关指标\n",
" df_indicators = []\n",
" for ts_code in df['ts_code'].unique():\n",
" df_index = df[df['ts_code'] == ts_code].copy()\n",
" df_index = calculate_indicators(df_index)\n",
" df_indicators.append(df_index)\n",
"\n",
" # 合并所有指数的结果\n",
" df_all_indicators = pd.concat(df_indicators, ignore_index=True)\n",
"\n",
" # 保留trade_date列，并将同一天的数据按ts_code合并成一行\n",
" df_final = df_all_indicators.pivot_table(\n",
" index='trade_date',\n",
" columns='ts_code',\n",
" values=['daily_return', \n",
" 'RSI', 'MACD', 'Signal_line', 'MACD_hist', \n",
" # 'sentiment_panic_greed_index',\n",
" 'up_ratio_20d', 'volume_change_rate', 'volatility',\n",
" 'amount_change_rate', 'amount_mean'],\n",
" aggfunc='last'\n",
" )\n",
"\n",
" df_final.columns = [f\"{col[1]}_{col[0]}\" for col in df_final.columns]\n",
" df_final = df_final.reset_index()\n",
"\n",
" return df_final\n",
"\n",
"\n",
"# 使用函数\n",
"h5_filename = '/mnt/d/PyProject/NewStock/data/index_data.h5'\n",
"index_data = generate_index_indicators(h5_filename)\n",
"index_data = index_data.dropna()\n"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "a735bc02ceb4d872",
"metadata": {
"ExecuteTime": {
"end_time": "2025-04-03T12:47:10.821169Z",
"start_time": "2025-04-03T12:47:10.751831Z"
}
},
"outputs": [],
"source": [
"import talib\n",
"import numpy as np"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "53f86ddc0677a6d7",
"metadata": {
"ExecuteTime": {
"end_time": "2025-04-03T12:47:15.944254Z",
"start_time": "2025-04-03T12:47:10.826179Z"
},
"jupyter": {
"source_hidden": true
},
"scrolled": true
},
"outputs": [],
"source": [
"from main.utils.factor import get_act_factor\n",
"\n",
"\n",
"def read_industry_data(h5_filename):\n",
" # 读取 H5 文件中所有的行业数据\n",
" industry_data = pd.read_hdf(h5_filename, key='sw_daily', columns=[\n",
" 'ts_code', 'trade_date', 'open', 'close', 'high', 'low', 'pe', 'pb', 'vol'\n",
" ]) # 假设 H5 文件的键是 'industry_data'\n",
" industry_data = industry_data.sort_values(by=['ts_code', 'trade_date'])\n",
" industry_data = industry_data.reindex()\n",
" industry_data['trade_date'] = pd.to_datetime(industry_data['trade_date'], format='%Y%m%d')\n",
"\n",
" grouped = industry_data.groupby('ts_code', group_keys=False)\n",
" industry_data['obv'] = grouped.apply(\n",
" lambda x: pd.Series(talib.OBV(x['close'].values, x['vol'].values), index=x.index)\n",
" )\n",
" industry_data['return_5'] = grouped['close'].apply(lambda x: x / x.shift(5) - 1)\n",
" industry_data['return_20'] = grouped['close'].apply(lambda x: x / x.shift(20) - 1)\n",
"\n",
" industry_data = get_act_factor(industry_data, cat=False)\n",
" industry_data = industry_data.sort_values(by=['trade_date', 'ts_code'])\n",
"\n",
" # # 计算每天每个 ts_code 的因子和当天所有 ts_code 的中位数的偏差\n",
" # factor_columns = ['obv', 'return_5', 'return_20', 'act_factor1', 'act_factor2', 'act_factor3', 'act_factor4'] # 因子列\n",
" # \n",
" # for factor in factor_columns:\n",
" # if factor in industry_data.columns:\n",
" # # 计算每天每个 ts_code 的因子值与当天所有 ts_code 的中位数的偏差\n",
" # industry_data[f'{factor}_deviation'] = industry_data.groupby('trade_date')[factor].transform(\n",
" # lambda x: x - x.mean())\n",
"\n",
" industry_data['return_5_percentile'] = industry_data.groupby('trade_date')['return_5'].transform(\n",
" lambda x: x.rank(pct=True))\n",
" industry_data['return_20_percentile'] = industry_data.groupby('trade_date')['return_20'].transform(\n",
" lambda x: x.rank(pct=True))\n",
"\n",
" # cs_rank_intraday_range(industry_data)\n",
" # cs_rank_close_pos_in_range(industry_data)\n",
"\n",
" industry_data = industry_data.drop(columns=['open', 'close', 'high', 'low', 'pe', 'pb', 'vol'])\n",
"\n",
" industry_data = industry_data.rename(\n",
" columns={col: f'industry_{col}' for col in industry_data.columns if col not in ['ts_code', 'trade_date']})\n",
"\n",
" industry_data = industry_data.rename(columns={'ts_code': 'cat_l2_code'})\n",
" return industry_data\n",
"\n",
"\n",
"industry_df = read_industry_data('/mnt/d/PyProject/NewStock/data/sw_daily.h5')\n"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "dbe2fd8021b9417f",
"metadata": {
"ExecuteTime": {
"end_time": "2025-04-03T12:47:15.969344Z",
"start_time": "2025-04-03T12:47:15.963327Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"['ts_code', 'open', 'close', 'high', 'low', 'amount', 'circ_mv', 'total_mv', 'is_st', 'up_limit', 'down_limit', 'buy_sm_vol', 'sell_sm_vol', 'buy_lg_vol', 'sell_lg_vol', 'buy_elg_vol', 'sell_elg_vol', 'net_mf_vol', 'his_low', 'his_high', 'cost_5pct', 'cost_15pct', 'cost_50pct', 'cost_85pct', 'cost_95pct', 'weight_avg', 'in_date']\n"
]
}
],
"source": [
"origin_columns = df.columns.tolist()\n",
"origin_columns = [col for col in origin_columns if\n",
" col not in ['turnover_rate', 'pe_ttm', 'volume_ratio', 'vol', 'pct_chg', 'l2_code', 'winner_rate']]\n",
"origin_columns = [col for col in origin_columns if col not in index_data.columns]\n",
"origin_columns = [col for col in origin_columns if 'cyq' not in col]\n",
"print(origin_columns)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "85c3e3d0235ffffa",
"metadata": {
"ExecuteTime": {
"end_time": "2025-04-03T12:47:16.089879Z",
"start_time": "2025-04-03T12:47:15.990101Z"
}
},
"outputs": [],
"source": [
"fina_indicator_df = read_and_merge_h5_data('/mnt/d/PyProject/NewStock/data/fina_indicator.h5', key='fina_indicator',\n",
" columns=['ts_code', 'ann_date', 'undist_profit_ps', 'ocfps', 'bps', 'roa', 'roe'],\n",
" df=None)\n",
"cashflow_df = read_and_merge_h5_data('/mnt/d/PyProject/NewStock/data/cashflow.h5', key='cashflow',\n",
" columns=['ts_code', 'ann_date', 'n_cashflow_act'],\n",
" df=None)\n",
"balancesheet_df = read_and_merge_h5_data('/mnt/d/PyProject/NewStock/data/balancesheet.h5', key='balancesheet',\n",
" columns=['ts_code', 'ann_date', 'money_cap', 'total_liab'],\n",
" df=None)\n",
"top_list_df = read_and_merge_h5_data('/mnt/d/PyProject/NewStock/data/top_list.h5', key='top_list',\n",
" columns=['ts_code', 'trade_date', 'reason'],\n",
" df=None)\n",
"\n",
"top_list_df = top_list_df.sort_values(by='trade_date', ascending=False).drop_duplicates(subset=['ts_code', 'trade_date'], keep='first').sort_values(by='trade_date')\n",
"\n",
"stk_holdertrade_df = read_and_merge_h5_data('/mnt/d/PyProject/NewStock/data/stk_holdertrade.h5', key='stk_holdertrade',\n",
" columns=['ts_code', 'ann_date', 'in_de', 'change_ratio', 'after_ratio'],\n",
" df=None)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "92d84ce15a562ec6",
"metadata": {
"ExecuteTime": {
"end_time": "2025-04-03T13:08:01.612695Z",
"start_time": "2025-04-03T12:47:16.121802Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"使用 'ann_date' 作为财务数据生效日期。\n",
"警告: 从 financial_data_subset 中移除了 366 行，因为其 'ts_code' 或 'ann_date' 列存在空值。\n",
"使用 'ann_date' 作为财务数据生效日期。\n",
"警告: 从 financial_data_subset 中移除了 366 行，因为其 'ts_code' 或 'ann_date' 列存在空值。\n",
"使用 'ann_date' 作为财务数据生效日期。\n",
"警告: 从 financial_data_subset 中移除了 366 行，因为其 'ts_code' 或 'ann_date' 列存在空值。\n",
"使用 'ann_date' 作为财务数据生效日期。\n",
"警告: 从 financial_data_subset 中移除了 366 行，因为其 'ts_code' 或 'ann_date' 列存在空值。\n",
"开始计算因子: AR, BR (原地修改)...\n",
"因子 AR, BR 计算成功。\n",
"因子 AR, BR 计算流程结束。\n",
"使用 'ann_date' 作为财务数据生效日期。\n",
"使用 'ann_date' 作为财务数据生效日期。\n",
"使用 'ann_date' 作为财务数据生效日期。\n",
"使用 'ann_date' 作为财务数据生效日期。\n",
"警告: 从 financial_data_subset 中移除了 366 行，因为其 'ts_code' 或 'ann_date' 列存在空值。\n",
"计算 BBI...\n",
"--- 计算日级别偏离度 (使用 pct_chg) ---\n",
"--- 计算日级别动量基准 (使用 pct_chg) ---\n",
"日级别动量基准计算完成 (使用 pct_chg)。\n",
"日级别偏离度计算完成 (使用 pct_chg)。\n",
"--- 计算日级别行业偏离度 (使用 pct_chg 和行业基准) ---\n",
"--- 计算日级别行业动量基准 (使用 pct_chg 和 cat_l2_code) ---\n",
"错误: 计算日级别行业动量基准需要以下列: ['pct_chg', 'cat_l2_code', 'trade_date', 'ts_code']。\n",
"错误: 计算日级别行业偏离度需要以下列: ['pct_chg', 'daily_industry_positive_benchmark', 'daily_industry_negative_benchmark']。请先运行 daily_industry_momentum_benchmark(df)。\n",
"Index(['ts_code', 'trade_date', 'open', 'close', 'high', 'low', 'vol',\n",
" 'amount', 'pct_chg', 'turnover_rate', 'pe_ttm', 'circ_mv', 'total_mv',\n",
" 'volume_ratio', 'is_st', 'up_limit', 'down_limit', 'buy_sm_vol',\n",
" 'sell_sm_vol', 'buy_lg_vol', 'sell_lg_vol', 'buy_elg_vol',\n",
" 'sell_elg_vol', 'net_mf_vol', 'his_low', 'his_high', 'cost_5pct',\n",
" 'cost_15pct', 'cost_50pct', 'cost_85pct', 'cost_95pct', 'weight_avg',\n",
" 'winner_rate', 'l2_code', 'undist_profit_ps', 'ocfps', 'roa', 'roe',\n",
" 'AR', 'BR', 'AR_BR', 'log_circ_mv', 'cashflow_to_ev_factor',\n",
" 'book_to_price_ratio', 'turnover_rate_mean_5', 'variance_20',\n",
" 'bbi_ratio_factor', 'daily_deviation', 'lg_elg_net_buy_vol',\n",
" 'flow_lg_elg_intensity', 'sm_net_buy_vol', 'flow_divergence_diff',\n",
" 'flow_divergence_ratio', 'total_buy_vol', 'lg_elg_buy_prop',\n",
" 'flow_struct_buy_change', 'lg_elg_net_buy_vol_change',\n",
" 'flow_lg_elg_accel', 'chip_concentration_range', 'chip_skewness',\n",
" 'floating_chip_proxy', 'cost_support_15pct_change',\n",
" 'cat_winner_price_zone', 'flow_chip_consistency',\n",
" 'profit_taking_vs_absorb', '_is_positive', '_is_negative',\n",
" 'cat_is_positive', '_pos_returns', '_neg_returns', '_pos_returns_sq',\n",
" '_neg_returns_sq', 'upside_vol', 'downside_vol', 'vol_ratio',\n",
" 'return_skew', 'return_kurtosis', 'volume_change_rate',\n",
" 'cat_volume_breakout', 'turnover_deviation', 'cat_turnover_spike',\n",
" 'avg_volume_ratio', 'cat_volume_ratio_breakout', 'vol_spike',\n",
" 'vol_std_5', 'atr_14', 'atr_6', 'obv'],\n",
" dtype='object')\n",
"Calculating lg_flow_mom_corr_20_60...\n",
"Finished lg_flow_mom_corr_20_60.\n",
"Calculating lg_flow_accel...\n",
"Finished lg_flow_accel.\n",
"Calculating profit_pressure...\n",
"Finished profit_pressure.\n",
"Calculating underwater_resistance...\n",
"Finished underwater_resistance.\n",
"Calculating cost_conc_std_20...\n",
"Finished cost_conc_std_20.\n",
"Calculating profit_decay_20...\n",
"Finished profit_decay_20.\n",
"Calculating vol_amp_loss_20...\n",
"Finished vol_amp_loss_20.\n",
"Calculating vol_drop_profit_cnt_5...\n",
"Finished vol_drop_profit_cnt_5.\n",
"Calculating lg_flow_vol_interact_20...\n",
"Finished lg_flow_vol_interact_20.\n",
"Calculating cost_break_confirm_cnt_5...\n",
"Finished cost_break_confirm_cnt_5.\n",
"Calculating atr_norm_channel_pos_14...\n",
"Finished atr_norm_channel_pos_14.\n",
"Calculating turnover_diff_skew_20...\n",
"Finished turnover_diff_skew_20.\n",
"Calculating lg_sm_flow_diverge_20...\n",
"Finished lg_sm_flow_diverge_20.\n",
"Calculating pullback_strong_20_20...\n",
"Finished pullback_strong_20_20.\n",
"Calculating vol_wgt_hist_pos_20...\n",
"Finished vol_wgt_hist_pos_20.\n",
"Calculating vol_adj_roc_20...\n",
"Finished vol_adj_roc_20.\n",
"Calculating cs_rank_net_lg_flow_val...\n",
"Finished cs_rank_net_lg_flow_val.\n",
"Calculating cs_rank_flow_divergence...\n",
"Finished cs_rank_flow_divergence.\n",
"Calculating cs_rank_ind_adj_lg_flow...\n",
"Finished cs_rank_ind_adj_lg_flow.\n",
"Calculating cs_rank_elg_buy_ratio...\n",
"Finished cs_rank_elg_buy_ratio.\n",
"Calculating cs_rank_rel_profit_margin...\n",
"Finished cs_rank_rel_profit_margin.\n",
"Calculating cs_rank_cost_breadth...\n",
"Finished cs_rank_cost_breadth.\n",
"Calculating cs_rank_dist_to_upper_cost...\n",
"Finished cs_rank_dist_to_upper_cost.\n",
"Calculating cs_rank_winner_rate...\n",
"Finished cs_rank_winner_rate.\n",
"Calculating cs_rank_intraday_range...\n",
"Finished cs_rank_intraday_range.\n",
"Calculating cs_rank_close_pos_in_range...\n",
"Finished cs_rank_close_pos_in_range.\n",
"Calculating cs_rank_opening_gap...\n",
"Error calculating cs_rank_opening_gap: Missing 'pre_close' column. Assigning NaN.\n",
"Calculating cs_rank_pos_in_hist_range...\n",
"Finished cs_rank_pos_in_hist_range.\n",
"Calculating cs_rank_vol_x_profit_margin...\n",
"Finished cs_rank_vol_x_profit_margin.\n",
"Calculating cs_rank_lg_flow_price_concordance...\n",
"Finished cs_rank_lg_flow_price_concordance.\n",
"Calculating cs_rank_turnover_per_winner...\n",
"Finished cs_rank_turnover_per_winner.\n",
"Calculating cs_rank_ind_cap_neutral_pe (Placeholder - requires statsmodels)...\n",
"Finished cs_rank_ind_cap_neutral_pe (Placeholder).\n",
"Calculating cs_rank_volume_ratio...\n",
"Finished cs_rank_volume_ratio.\n",
"Calculating cs_rank_elg_buy_sell_sm_ratio...\n",
"Finished cs_rank_elg_buy_sell_sm_ratio.\n",
"Calculating cs_rank_cost_dist_vol_ratio...\n",
"Finished cs_rank_cost_dist_vol_ratio.\n",
"Calculating cs_rank_size...\n",
"Finished cs_rank_size.\n",
"<class 'pandas.core.frame.DataFrame'>\n",
"RangeIndex: 4987750 entries, 0 to 4987749\n",
"Columns: 181 entries, ts_code to cs_rank_size\n",
"dtypes: bool(10), datetime64[ns](1), float64(165), int64(3), object(2)\n",
"memory usage: 6.4+ GB\n",
"None\n",
"['ts_code', 'trade_date', 'open', 'close', 'high', 'low', 'vol', 'amount', 'pct_chg', 'turnover_rate', 'pe_ttm', 'circ_mv', 'total_mv', 'volume_ratio', 'is_st', 'up_limit', 'down_limit', 'buy_sm_vol', 'sell_sm_vol', 'buy_lg_vol', 'sell_lg_vol', 'buy_elg_vol', 'sell_elg_vol', 'net_mf_vol', 'his_low', 'his_high', 'cost_5pct', 'cost_15pct', 'cost_50pct', 'cost_85pct', 'cost_95pct', 'weight_avg', 'winner_rate', 'cat_l2_code', 'undist_profit_ps', 'ocfps', 'roa', 'roe', 'AR', 'BR', 'AR_BR', 'log_circ_mv', 'cashflow_to_ev_factor', 'book_to_price_ratio', 'turnover_rate_mean_5', 'variance_20', 'bbi_ratio_factor', 'daily_deviation', 'lg_elg_net_buy_vol', 'flow_lg_elg_intensity', 'sm_net_buy_vol', 'flow_divergence_diff', 'flow_divergence_ratio', 'total_buy_vol', 'lg_elg_buy_prop', 'flow_struct_buy_change', 'lg_elg_net_buy_vol_change', 'flow_lg_elg_accel', 'chip_concentration_range', 'chip_skewness', 'floating_chip_proxy', 'cost_support_15pct_change', 'cat_winner_price_zone', 'flow_chip_consistency', 'profit_taking_vs_absorb', 'cat_is_positive', 'upside_vol', 'downside_vol', 'vol_ratio', 'return_skew', 'return_kurtosis', 'volume_change_rate', 'cat_volume_breakout', 'turnover_deviation', 'cat_turnover_spike', 'avg_volume_ratio', 'cat_volume_ratio_breakout', 'vol_spike', 'vol_std_5', 'atr_14', 'atr_6', 'obv', 'maobv_6', 'rsi_3', 'return_5', 'return_20', 'std_return_5', 'std_return_90', 'std_return_90_2', 'act_factor1', 'act_factor2', 'act_factor3', 'act_factor4', 'rank_act_factor1', 'rank_act_factor2', 'rank_act_factor3', 'cov', 'delta_cov', 'alpha_22_improved', 'alpha_003', 'alpha_007', 'alpha_013', 'vol_break', 'weight_roc5', 'price_cost_divergence', 'smallcap_concentration', 'cost_stability', 'high_cost_break_days', 'liquidity_risk', 'turnover_std', 'mv_volatility', 'volume_growth', 'mv_growth', 'momentum_factor', 'resonance_factor', 'log_close', 'cat_vol_spike', 'up', 'down', 'obv_maobv_6', 'std_return_5_over_std_return_90', 'std_return_90_minus_std_return_90_2', 'cat_af2', 'cat_af3', 'cat_af4', 'act_factor5', 'act_factor6', 'active_buy_volume_large', 'active_buy_volume_big', 'active_buy_volume_small', 'buy_lg_vol_minus_sell_lg_vol', 'buy_elg_vol_minus_sell_elg_vol', 'ctrl_strength', 'low_cost_dev', 'asymmetry', 'lock_factor', 'cat_vol_break', 'cost_atr_adj', 'cat_golden_resonance', 'mv_turnover_ratio', 'mv_adjusted_volume', 'mv_weighted_turnover', 'nonlinear_mv_volume', 'mv_volume_ratio', 'mv_momentum', 'lg_flow_mom_corr_20_60', 'lg_flow_accel', 'profit_pressure', 'underwater_resistance', 'cost_conc_std_20', 'profit_decay_20', 'vol_amp_loss_20', 'vol_drop_profit_cnt_5', 'lg_flow_vol_interact_20', 'cost_break_confirm_cnt_5', 'atr_norm_channel_pos_14', 'turnover_diff_skew_20', 'lg_sm_flow_diverge_20', 'pullback_strong_20_20', 'vol_wgt_hist_pos_20', 'vol_adj_roc_20', 'cs_rank_net_lg_flow_val', 'cs_rank_flow_divergence', 'cs_rank_ind_adj_lg_flow', 'cs_rank_elg_buy_ratio', 'cs_rank_rel_profit_margin', 'cs_rank_cost_breadth', 'cs_rank_dist_to_upper_cost', 'cs_rank_winner_rate', 'cs_rank_intraday_range', 'cs_rank_close_pos_in_range', 'cs_rank_opening_gap', 'cs_rank_pos_in_hist_range', 'cs_rank_vol_x_profit_margin', 'cs_rank_lg_flow_price_concordance', 'cs_rank_turnover_per_winner', 'cs_rank_ind_cap_neutral_pe', 'cs_rank_volume_ratio', 'cs_rank_elg_buy_sell_sm_ratio', 'cs_rank_cost_dist_vol_ratio', 'cs_rank_size']\n"
]
}
],
"source": [
"\n",
"import numpy as np\n",
"from main.factor.factor import *\n",
"from main.factor.money_factor import *\n",
"\n",
"\n",
"def filter_data(df):\n",
" # df = df.groupby('trade_date').apply(lambda x: x.nlargest(1000, 'act_factor1'))\n",
" # df = df[df['trade_date'] <= '2025-06-01']\n",
" df = df[~df['is_st']]\n",
" df = df[~df['ts_code'].str.endswith('BJ')]\n",
" df = df[~df['ts_code'].str.startswith('30')]\n",
" df = df[~df['ts_code'].str.startswith('68')]\n",
" df = df[~df['ts_code'].str.startswith('8')]\n",
" df = df[df['trade_date'] >= '2019-01-01']\n",
" if 'in_date' in df.columns:\n",
" df = df.drop(columns=['in_date'])\n",
" df = df.reset_index(drop=True)\n",
" return df\n",
"\n",
"gc.collect()\n",
"\n",
"df = filter_data(df)\n",
"df = df.sort_values(by=['ts_code', 'trade_date'])\n",
"\n",
"# df = price_minus_deduction_price(df, n=120)\n",
"# df = price_deduction_price_diff_ratio_to_sma(df, n=120)\n",
"# df = cat_price_vs_sma_vs_deduction_price(df, n=120)\n",
"# df = cat_reason(df, top_list_df)\n",
"# df = cat_is_on_top_list(df, top_list_df)\n",
"\n",
"# df = ts_turnover_rate_acceleration_5_20(df)\n",
"# df = ts_vol_sustain_10_30(df)\n",
"# df = cs_turnover_rate_relative_strength_20(df)\n",
"# df = cs_amount_outlier_10(df)\n",
"# df = holder_trade_factors(stk_holdertrade_df, df)\n",
"\n",
"df = add_financial_factor(df, fina_indicator_df, factor_value_col='undist_profit_ps')\n",
"df = add_financial_factor(df, fina_indicator_df, factor_value_col='ocfps')\n",
"df = add_financial_factor(df, fina_indicator_df, factor_value_col='roa')\n",
"df = add_financial_factor(df, fina_indicator_df, factor_value_col='roe')\n",
"\n",
"calculate_arbr(df, N=26)\n",
"df['log_circ_mv'] = np.log(df['circ_mv'])\n",
"df = calculate_cashflow_to_ev_factor(df, cashflow_df, balancesheet_df)\n",
"df = caculate_book_to_price_ratio(df, fina_indicator_df)\n",
"\n",
"df = turnover_rate_n(df, n=5)\n",
"df = variance_n(df, n=20)\n",
"df = bbi_ratio_factor(df)\n",
"df = daily_deviation(df)\n",
"df = daily_industry_deviation(df)\n",
"df, _ = get_rolling_factor(df)\n",
"df, _ = get_simple_factor(df)\n",
"\n",
"df = df.rename(columns={'l1_code': 'cat_l1_code'})\n",
"df = df.rename(columns={'l2_code': 'cat_l2_code'})\n",
"\n",
"lg_flow_mom_corr(df, N=20, M=60)\n",
"lg_flow_accel(df)\n",
"profit_pressure(df)\n",
"underwater_resistance(df)\n",
"cost_conc_std(df, N=20)\n",
"profit_decay(df, N=20)\n",
"vol_amp_loss(df, N=20)\n",
"vol_drop_profit_cnt(df, N=20, M=5)\n",
"lg_flow_vol_interact(df, N=20)\n",
"cost_break_confirm_cnt(df, M=5)\n",
"atr_norm_channel_pos(df, N=14)\n",
"turnover_diff_skew(df, N=20)\n",
"lg_sm_flow_diverge(df, N=20)\n",
"pullback_strong(df, N=20, M=20)\n",
"vol_wgt_hist_pos(df, N=20)\n",
"vol_adj_roc(df, N=20)\n",
"\n",
"cs_rank_net_lg_flow_val(df)\n",
"cs_rank_flow_divergence(df)\n",
"cs_rank_industry_adj_lg_flow(df) # Needs cat_l2_code\n",
"cs_rank_elg_buy_ratio(df)\n",
"cs_rank_rel_profit_margin(df)\n",
"cs_rank_cost_breadth(df)\n",
"cs_rank_dist_to_upper_cost(df)\n",
"cs_rank_winner_rate(df)\n",
"cs_rank_intraday_range(df)\n",
"cs_rank_close_pos_in_range(df)\n",
"cs_rank_opening_gap(df) # Needs pre_close\n",
"cs_rank_pos_in_hist_range(df) # Needs his_low, his_high\n",
"cs_rank_vol_x_profit_margin(df)\n",
"cs_rank_lg_flow_price_concordance(df)\n",
"cs_rank_turnover_per_winner(df)\n",
"cs_rank_ind_cap_neutral_pe(df) # Placeholder - needs external libraries\n",
"cs_rank_volume_ratio(df) # Needs volume_ratio\n",
"cs_rank_elg_buy_sell_sm_ratio(df)\n",
"cs_rank_cost_dist_vol_ratio(df) # Needs volume_ratio\n",
"cs_rank_size(df) # Needs circ_mv\n",
"\n",
"\n",
"# df = df.merge(index_data, on='trade_date', how='left')\n",
"\n",
"print(df.info())\n",
"print(df.columns.tolist())"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "b87b938028afa206",
"metadata": {
"ExecuteTime": {
"end_time": "2025-04-03T13:08:03.658725Z",
"start_time": "2025-04-03T13:08:02.469611Z"
}
},
"outputs": [],
"source": [
"from scipy.stats import ks_2samp, wasserstein_distance\n",
"\n",
"\n",
"def remove_shifted_features(train_data, test_data, feature_columns, ks_threshold=0.05, wasserstein_threshold=0.1,\n",
" importance_threshold=0.05):\n",
" dropped_features = []\n",
"\n",
" # **统计数据漂移**\n",
" numeric_columns = train_data.select_dtypes(include=['float64', 'int64']).columns\n",
" numeric_columns = [col for col in numeric_columns if col in feature_columns]\n",
" for feature in numeric_columns:\n",
" ks_stat, p_value = ks_2samp(train_data[feature], test_data[feature])\n",
" wasserstein_dist = wasserstein_distance(train_data[feature], test_data[feature])\n",
"\n",
" if p_value < ks_threshold or wasserstein_dist > wasserstein_threshold:\n",
" dropped_features.append(feature)\n",
"\n",
" print(f\"检测到 {len(dropped_features)} 个可能漂移的特征: {dropped_features}\")\n",
"\n",
" # **应用阈值进行最终筛选**\n",
" filtered_features = [f for f in feature_columns if f not in dropped_features]\n",
"\n",
" return filtered_features, dropped_features\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "f4f16d63ad18d1bc",
"metadata": {
"ExecuteTime": {
"end_time": "2025-04-03T13:08:03.670700Z",
"start_time": "2025-04-03T13:08:03.665739Z"
}
},
"outputs": [],
"source": [
"import numpy as np\n",
"import statsmodels.api as sm # 用于中性化回归\n",
"from tqdm import tqdm # 可选，用于显示进度条\n",
"\n",
"# --- 常量 ---\n",
"epsilon = 1e-10 # 防止除零\n",
"\n",
"# --- 1. 中位数去极值 (MAD) ---\n",
"\n",
"def cs_mad_filter(df: pd.DataFrame,\n",
" features: list,\n",
" k: float = 3.0,\n",
" scale_factor: float = 1.4826):\n",
" \"\"\"\n",
" 对指定特征列进行截面 MAD 去极值处理 (原地修改)。\n",
"\n",
" 方法: 对每日截面数据，计算 median 和 MAD，\n",
" 将超出 [median - k * scale * MAD, median + k * scale * MAD] 范围的值\n",
" 替换为边界值 (Winsorization)。\n",
" scale_factor=1.4826 使得 MAD 约等于正态分布的标准差。\n",
"\n",
" Args:\n",
" df (pd.DataFrame): 输入 DataFrame，需包含 'trade_date' 和 features 列。\n",
" features (list): 需要处理的特征列名列表。\n",
" k (float): MAD 的倍数，用于确定边界。默认为 3.0。\n",
" scale_factor (float): MAD 的缩放因子。默认为 1.4826。\n",
"\n",
" WARNING: 此函数会原地修改输入的 DataFrame 'df'。\n",
" \"\"\"\n",
" print(f\"开始截面 MAD 去极值处理 (k={k})...\")\n",
" if not all(col in df.columns for col in features):\n",
" missing = [col for col in features if col not in df.columns]\n",
" print(f\"错误: DataFrame 中缺少以下特征列: {missing}。跳过去极值处理。\")\n",
" return\n",
"\n",
" grouped = df.groupby('trade_date')\n",
"\n",
" for col in tqdm(features, desc=\"MAD Filtering\"):\n",
" try:\n",
" # 计算截面中位数\n",
" median = grouped[col].transform('median')\n",
" # 计算截面 MAD (Median Absolute Deviation from Median)\n",
" mad = (df[col] - median).abs().groupby(df['trade_date']).transform('median')\n",
"\n",
" # 计算上下边界\n",
" lower_bound = median - k * scale_factor * mad\n",
" upper_bound = median + k * scale_factor * mad\n",
"\n",
" # 原地应用 clip\n",
" df[col] = np.clip(df[col], lower_bound, upper_bound)\n",
"\n",
" except KeyError:\n",
" print(f\"警告: 列 '{col}' 可能不存在或在分组中出错，跳过此列的 MAD 处理。\")\n",
" except Exception as e:\n",
" print(f\"警告: 处理列 '{col}' 时发生错误: {e}，跳过此列的 MAD 处理。\")\n",
"\n",
" print(\"截面 MAD 去极值处理完成。\")\n",
"\n",
"\n",
"# --- 2. 行业市值中性化 ---\n",
"\n",
"from tqdm import tqdm\n",
"\n",
"def cs_neutralize_market_cap_numpy(df: pd.DataFrame,\n",
" features: list,\n",
" market_cap_col: str = 'circ_mv'):\n",
" \"\"\"\n",
" 对 DataFrame 中的指定特征进行截面市值中性化 (NumPy 优化)。\n",
"\n",
" Args:\n",
" df (pd.DataFrame): 包含数据的 DataFrame，需要有 'trade_date' 和 market_cap_col 列。\n",
" features (list): 需要进行市值中性化的特征列名列表。\n",
" market_cap_col (str): 包含市值数据的列名，默认为 'circ_mv'。\n",
" \"\"\"\n",
" print(\"开始截面市值中性化 (NumPy 优化)...\")\n",
" required_cols = features + ['trade_date', market_cap_col]\n",
" if not all(col in df.columns for col in required_cols):\n",
" missing = [col for col in required_cols if col not in df.columns]\n",
" print(f\"错误: DataFrame 中缺少必需列: {missing}。无法进行中性化。\")\n",
" return\n",
"\n",
" df_copy = df\n",
" log_cap_col = '_log_market_cap'\n",
" df_copy[log_cap_col] = np.log1p(df_copy[market_cap_col])\n",
"\n",
" # 创建一个 DataFrame 来存储所有日期的残差结果\n",
" residuals_container = pd.DataFrame(index=df_copy.index, columns=features, dtype=float)\n",
"\n",
" for date, group_df in tqdm(df_copy.groupby('trade_date'), desc=\"Neutralizing by Date (NumPy)\"):\n",
" # 准备 X 矩阵 (自变量)：常数项和对数市值\n",
" X_daily = np.concatenate([np.ones((len(group_df), 1)), group_df[[log_cap_col]].values], axis=1)\n",
"\n",
" for feature_col in features:\n",
" Y_daily = group_df[feature_col].values\n",
"\n",
" # 处理 NaN：只对有效数据对进行回归\n",
" valid_mask_y = ~np.isnan(Y_daily)\n",
" valid_mask_x = ~np.isnan(X_daily).any(axis=1)\n",
" valid_mask = valid_mask_y & valid_mask_x\n",
"\n",
" current_feature_indices = group_df.index[valid_mask]\n",
"\n",
" if np.sum(valid_mask) < X_daily.shape[1] + 1:\n",
" # 有效数据不足，此特征在此日期保持 NaN\n",
" continue\n",
"\n",
" Y_valid = Y_daily[valid_mask]\n",
" X_valid = X_daily[valid_mask, :]\n",
"\n",
" try:\n",
" # 使用 np.linalg.lstsq 进行 OLS 计算\n",
" beta, sum_sq_resid, rank, s = np.linalg.lstsq(X_valid, Y_valid, rcond=None)\n",
"\n",
" # 计算预测值 Y_hat = X_valid @ beta\n",
" Y_hat_valid = X_valid @ beta\n",
"\n",
" # 计算残差 residuals = Y_valid - Y_hat_valid\n",
" residuals_valid = Y_valid - Y_hat_valid\n",
"\n",
" # 将计算得到的残差放回 residuals_container\n",
" residuals_container.loc[current_feature_indices, feature_col] = residuals_valid\n",
"\n",
" except np.linalg.LinAlgError:\n",
" pass\n",
" except Exception as e:\n",
" pass\n",
"\n",
" # 将所有计算得到的残差更新回原始的 df (原地修改)\n",
" for feature_col in features:\n",
" df[feature_col] = residuals_container[feature_col]\n",
"\n",
" # 清理临时列\n",
" df.drop(columns=[log_cap_col], inplace=True, errors='ignore')\n",
" print(\"截面市值中性化完成 (NumPy 优化)。\")\n",
"\n",
"# --- 3. Z-Score 标准化 ---\n",
"\n",
"def cs_zscore_standardize(df: pd.DataFrame, features: list, epsilon: float = 1e-10):\n",
" \"\"\"\n",
" 对指定特征列进行截面 Z-Score 标准化 (原地修改)。\n",
" 方法: Z = (value - cross_sectional_mean) / (cross_sectional_std + epsilon)\n",
"\n",
" Args:\n",
" df (pd.DataFrame): 输入 DataFrame，需包含 'trade_date' 和 features 列。\n",
" features (list): 需要处理的特征列名列表。\n",
" epsilon (float): 防止除以零的小常数。\n",
"\n",
" WARNING: 此函数会原地修改输入的 DataFrame 'df'。\n",
" \"\"\"\n",
" print(\"开始截面 Z-Score 标准化...\")\n",
" if not all(col in df.columns for col in features):\n",
" missing = [col for col in features if col not in df.columns]\n",
" print(f\"错误: DataFrame 中缺少以下特征列: {missing}。跳过标准化处理。\")\n",
" return\n",
"\n",
" grouped = df.groupby('trade_date')\n",
"\n",
" for col in tqdm(features, desc=\"Standardizing\"):\n",
" try:\n",
" # 使用 transform 计算截面均值和标准差\n",
" mean = grouped[col].transform('mean')\n",
" std = grouped[col].transform('std')\n",
"\n",
" # 计算 Z-Score 并原地赋值\n",
" df[col] = (df[col] - mean) / (std + epsilon)\n",
"\n",
" except KeyError:\n",
" print(f\"警告: 列 '{col}' 可能不存在或在分组中出错，跳过此列的标准化处理。\")\n",
" except Exception as e:\n",
" print(f\"警告: 处理列 '{col}' 时发生错误: {e}，跳过此列的标准化处理。\")\n",
"\n",
" print(\"截面 Z-Score 标准化完成。\")\n",
"\n",
"def fill_nan_with_daily_median(df: pd.DataFrame, feature_columns: list[str]) -> pd.DataFrame:\n",
" \"\"\"\n",
" 对指定特征列进行每日截面中位数填充缺失值 (NaN)。\n",
"\n",
" 参数:\n",
" df (pd.DataFrame): 包含多日数据的DataFrame，需要包含 'trade_date' 和 feature_columns 中的列。\n",
" feature_columns (list[str]): 需要进行缺失值填充的特征列名称列表。\n",
"\n",
" 返回:\n",
" pd.DataFrame: 包含缺失值填充后特征列的DataFrame。在输入DataFrame的副本上操作。\n",
" \"\"\"\n",
" processed_df = df.copy() # 在副本上操作，保留原始数据\n",
"\n",
" # 确保 trade_date 是 datetime 类型以便正确分组\n",
" processed_df['trade_date'] = pd.to_datetime(processed_df['trade_date'])\n",
"\n",
" def _fill_daily_nan(group):\n",
" # group 是某一个交易日的 DataFrame\n",
"\n",
" # 遍历指定的特征列\n",
" for feature_col in feature_columns:\n",
" # 检查列是否存在于当前分组中\n",
" if feature_col in group.columns:\n",
" # 计算当日该特征的中位数\n",
" median_val = group[feature_col].median()\n",
"\n",
" # 使用当日中位数填充该特征列的 NaN 值\n",
" # inplace=True 会直接修改 group DataFrame\n",
" group[feature_col].fillna(median_val, inplace=True)\n",
" # else:\n",
" # print(f\"Warning: Feature column '{feature_col}' not found in daily group for {group['trade_date'].iloc[0]}. Skipping.\")\n",
"\n",
" return group\n",
"\n",
" # 按交易日期分组，并应用每日填充函数\n",
" # group_keys=False 避免将分组键添加到结果索引中\n",
" filled_df = processed_df.groupby('trade_date', group_keys=False).apply(_fill_daily_nan)\n",
"\n",
" return filled_df"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "40e6b68a91b30c79",
"metadata": {
"ExecuteTime": {
"end_time": "2025-04-03T13:08:04.694262Z",
"start_time": "2025-04-03T13:08:03.694904Z"
}
},
"outputs": [],
"source": [
"def remove_outliers_label_percentile(label: pd.Series, lower_percentile: float = 0.01, upper_percentile: float = 0.99,\n",
" log=True):\n",
" if not (0 <= lower_percentile < upper_percentile <= 1):\n",
" raise ValueError(\"Percentile values must satisfy 0 <= lower_percentile < upper_percentile <= 1.\")\n",
"\n",
" # Calculate lower and upper bounds based on percentiles\n",
" lower_bound = label.quantile(lower_percentile)\n",
" upper_bound = label.quantile(upper_percentile)\n",
"\n",
" # Filter out values outside the bounds\n",
" filtered_label = label[(label >= lower_bound) & (label <= upper_bound)]\n",
"\n",
" # Print the number of removed outliers\n",
" if log:\n",
" print(f\"Removed {len(label) - len(filtered_label)} outliers.\")\n",
" return filtered_label\n",
"\n",
"\n",
"def calculate_risk_adjusted_target(df, days=5):\n",
" df = df.sort_values(by=['ts_code', 'trade_date'])\n",
"\n",
" df['future_close'] = df.groupby('ts_code')['close'].shift(-days)\n",
" df['future_open'] = df.groupby('ts_code')['open'].shift(-1)\n",
" df['future_return'] = (df['future_close'] - df['future_open']) / df['future_open']\n",
"\n",
" df['future_volatility'] = df.groupby('ts_code')['future_return'].rolling(days, min_periods=1).std().reset_index(\n",
" level=0, drop=True)\n",
" sharpe_ratio = df['future_return'] * df['future_volatility']\n",
" sharpe_ratio.replace([np.inf, -np.inf], np.nan, inplace=True)\n",
"\n",
" return sharpe_ratio\n",
"\n",
"\n",
"def calculate_score(df, days=5, lambda_param=1.0):\n",
" def calculate_max_drawdown(prices):\n",
" peak = prices.iloc[0] # 初始化峰值\n",
" max_drawdown = 0 # 初始化最大回撤\n",
"\n",
" for price in prices:\n",
" if price > peak:\n",
" peak = price # 更新峰值\n",
" else:\n",
" drawdown = (peak - price) / peak # 计算当前回撤\n",
" max_drawdown = max(max_drawdown, drawdown) # 更新最大回撤\n",
"\n",
" return max_drawdown\n",
"\n",
" def compute_stock_score(stock_df):\n",
" stock_df = stock_df.sort_values(by=['trade_date'])\n",
" future_return = stock_df['future_return']\n",
" # 使用已有的 pct_chg 字段计算波动率\n",
" volatility = stock_df['pct_chg'].rolling(days).std().shift(-days)\n",
" max_drawdown = stock_df['close'].rolling(days).apply(calculate_max_drawdown, raw=False).shift(-days)\n",
" score = future_return - lambda_param * max_drawdown\n",
" return score\n",
"\n",
" # # 确保 DataFrame 按照股票代码和交易日期排序\n",
" # df = df.sort_values(by=['ts_code', 'trade_date'])\n",
"\n",
" # 对每个股票分别计算 score\n",
" df['score'] = df.groupby('ts_code').apply(compute_stock_score).reset_index(level=0, drop=True)\n",
"\n",
" return df['score']\n",
"\n",
"\n",
"def remove_highly_correlated_features(df, feature_columns, threshold=0.9):\n",
" numeric_features = df[feature_columns].select_dtypes(include=[np.number]).columns.tolist()\n",
" if not numeric_features:\n",
" raise ValueError(\"No numeric features found in the provided data.\")\n",
"\n",
" corr_matrix = df[numeric_features].corr().abs()\n",
" upper = corr_matrix.where(np.triu(np.ones(corr_matrix.shape), k=1).astype(bool))\n",
" to_drop = [column for column in upper.columns if any(upper[column] > threshold)]\n",
" remaining_features = [col for col in feature_columns if col not in to_drop\n",
" or 'act' in col or 'af' in col]\n",
" return remaining_features\n",
"\n",
"\n",
"def cross_sectional_standardization(df, features):\n",
" df_sorted = df.sort_values(by='trade_date') # 按时间排序\n",
" df_standardized = df_sorted.copy()\n",
"\n",
" for date in df_sorted['trade_date'].unique():\n",
" # 获取当前时间点的数据\n",
" current_data = df_standardized[df_standardized['trade_date'] == date]\n",
"\n",
" # 只对指定特征进行标准化\n",
" scaler = StandardScaler()\n",
" standardized_values = scaler.fit_transform(current_data[features])\n",
"\n",
" # 将标准化结果重新赋值回去\n",
" df_standardized.loc[df_standardized['trade_date'] == date, features] = standardized_values\n",
"\n",
" return df_standardized\n",
"\n",
"\n",
"def neutralize_manual_revised(df: pd.DataFrame, features: list, industry_col: str, mkt_cap_col: str) -> pd.DataFrame:\n",
" \"\"\"\n",
" 手动实现简单回归以提升速度，通过构建 Series 确保索引对齐。\n",
" 对特征在行业内部进行市值中性化。\n",
"\n",
" Args:\n",
" df: 输入的 DataFrame，包含特征、行业分类和市值列。\n",
" features: 需要进行中性化的特征列名列表。\n",
" industry_col: 行业分类列的列名。\n",
" mkt_cap_col: 市值列的列名。\n",
"\n",
" Returns:\n",
" 中性化后的 DataFrame。\n",
" \"\"\"\n",
"\n",
" df[mkt_cap_col] = pd.to_numeric(df[mkt_cap_col], errors='coerce')\n",
" df_cleaned = df.dropna(subset=[mkt_cap_col]).copy()\n",
" df_cleaned = df_cleaned[df_cleaned[mkt_cap_col] > 0].copy()\n",
"\n",
" if df_cleaned.empty:\n",
" print(\"警告: 清理市值异常值后 DataFrame 为空。\")\n",
" return df # 返回原始或空df，取决于清理前的状态\n",
"\n",
" processed_df = df\n",
"\n",
" for col in features:\n",
" if col not in df_cleaned.columns:\n",
" print(f\"警告: 特征列 '{col}' 不存在于清理后的 DataFrame 中，已跳过。\")\n",
" # 对于原始 df 中该列不存在的，在结果 df 中也保持原样（可能全是NaN）\n",
" processed_df[col] = df[col] if col in df.columns else np.nan\n",
" continue\n",
"\n",
" # 跳过对控制变量本身进行中性化\n",
" if col == mkt_cap_col or col == industry_col:\n",
" print(f\"警告: 特征列 '{col}' 是控制变量或内部使用的列，跳过中性化。\")\n",
" # 在结果 df 中也保持原样\n",
" processed_df[col] = df[col] if col in df.columns else np.nan\n",
" continue\n",
"\n",
" residual_series = pd.Series(index=df_cleaned.index, dtype=float)\n",
"\n",
" # 在分组前处理特征列的 NaN，只对有因子值的行进行回归计算\n",
" df_subset_factor = df_cleaned.dropna(subset=[col]).copy()\n",
"\n",
" if not df_subset_factor.empty:\n",
" for industry, group in df_subset_factor.groupby(industry_col):\n",
" x = group[mkt_cap_col] # 市值对数\n",
" y = group[col] # 因子值\n",
"\n",
" # 确保有足够的数据点 (>1) 且市值对数有方差 (>0) 进行回归计算\n",
" # 检查 np.var > 一个很小的正数，避免浮点数误差导致的零方差判断问题\n",
" if len(group) > 1 and np.var(x) > 1e-9:\n",
" try:\n",
" beta = np.cov(y, x)[0, 1] / np.var(x)\n",
" alpha = np.mean(y) - beta * np.mean(x)\n",
"\n",
" # 计算残差\n",
" resid = y - (alpha + beta * x)\n",
"\n",
" # 将计算出的残差存储到 residual_series 中，通过索引自动对齐\n",
" residual_series.loc[resid.index] = resid\n",
"\n",
" except Exception as e:\n",
" # 捕获可能的计算异常，例如np.cov或np.var因为极端数据报错\n",
" print(f\"警告: 在行业 {industry} 计算回归时发生错误: {e}。该组残差将设为原始值或 NaN。\")\n",
" # 此时该组的残差会保持 residual_series 初始化时的 NaN 或后续处理\n",
" # 也可以选择保留原始值：residual_series.loc[group.index] = group[col]\n",
"\n",
" else:\n",
" residual_series.loc[group.index] = group[col] # 保留原始因子值\n",
" processed_df.loc[residual_series.index, col] = residual_series\n",
"\n",
"\n",
" else:\n",
" processed_df[col] = np.nan # 或 df[col] if col in df.columns else np.nan\n",
"\n",
" return processed_df\n",
"\n",
"\n",
"import gc\n",
"\n",
"gc.collect()\n",
"\n",
"\n",
"def mad_filter(df, features, n=3):\n",
" for col in features:\n",
" median = df[col].median()\n",
" mad = np.median(np.abs(df[col] - median))\n",
" upper = median + n * mad\n",
" lower = median - n * mad\n",
" df[col] = np.clip(df[col], lower, upper) # 截断极值\n",
" return df\n",
"\n",
"\n",
"def percentile_filter(df, features, lower_percentile=0.01, upper_percentile=0.99):\n",
" for col in features:\n",
" # 按日期分组计算上下百分位数\n",
" lower_bound = df.groupby('trade_date')[col].transform(\n",
" lambda x: x.quantile(lower_percentile)\n",
" )\n",
" upper_bound = df.groupby('trade_date')[col].transform(\n",
" lambda x: x.quantile(upper_percentile)\n",
" )\n",
" # 截断超出范围的值\n",
" df[col] = np.clip(df[col], lower_bound, upper_bound)\n",
" return df\n",
"\n",
"\n",
"from scipy.stats import iqr\n",
"\n",
"\n",
"def iqr_filter(df, features):\n",
" for col in features:\n",
" df[col] = df.groupby('trade_date')[col].transform(\n",
" lambda x: (x - x.median()) / iqr(x) if iqr(x) != 0 else x\n",
" )\n",
" return df\n",
"\n",
"\n",
"def quantile_filter(df, features, lower_quantile=0.01, upper_quantile=0.99, window=60):\n",
" df = df.copy()\n",
" for col in features:\n",
" # 计算 rolling 统计量，需要按日期进行 groupby\n",
" rolling_lower = df.groupby('trade_date')[col].transform(lambda x: x.rolling(window=min(len(x), window)).quantile(lower_quantile))\n",
" rolling_upper = df.groupby('trade_date')[col].transform(lambda x: x.rolling(window=min(len(x), window)).quantile(upper_quantile))\n",
"\n",
" # 对数据进行裁剪\n",
" df[col] = np.clip(df[col], rolling_lower, rolling_upper)\n",
" \n",
" return df\n",
"\n",
"def select_top_features_by_rankic(df: pd.DataFrame, feature_columns: list, n: int, target_column: str = 'future_return') -> list:\n",
" \"\"\"\n",
" 计算给定特征与目标列的 RankIC，并返回 RankIC 绝对值最高的 n 个特征。\n",
"\n",
" Args:\n",
" df: 包含特征列和目标列的 Pandas DataFrame。\n",
" feature_columns: 包含所有待评估特征列名的列表。\n",
" n: 希望选取的 RankIC 绝对值最高的特征数量。\n",
" target_column: 目标列的名称，用于计算 RankIC。默认为 'future_return'。\n",
"\n",
" Returns:\n",
" 包含 RankIC 绝对值最高的 n 个特征列名的列表。\n",
" \"\"\"\n",
" numeric_columns = df.select_dtypes(include=['float64', 'int64']).columns\n",
" numeric_columns = [col for col in numeric_columns if col in feature_columns]\n",
" if target_column not in df.columns:\n",
" raise ValueError(f\"目标列 '{target_column}' 不存在于 DataFrame 中。\")\n",
"\n",
" rankic_scores = {}\n",
" for feature in numeric_columns:\n",
" if feature not in df.columns:\n",
" print(f\"警告: 特征列 '{feature}' 不存在于 DataFrame 中，已跳过。\")\n",
" continue\n",
"\n",
" # 计算特征与目标列的 RankIC (斯皮尔曼相关系数)\n",
" # dropna() 是为了处理缺失值，确保相关性计算不失败\n",
" valid_data = df[[feature, target_column]].dropna()\n",
" if len(valid_data) > 1: # 确保有足够的数据点进行相关性计算\n",
" # 计算斯皮尔曼相关性\n",
" correlation = valid_data[feature].corr(valid_data[target_column], method='spearman')\n",
" rankic_scores[feature] = abs(correlation) # 使用绝对值来衡量相关性强度\n",
" else:\n",
" rankic_scores[feature] = 0 # 数据不足，RankIC设为0或跳过\n",
"\n",
" # 将 RankIC 分数转换为 Series 便于排序\n",
" rankic_series = pd.Series(rankic_scores)\n",
"\n",
" # 按 RankIC 绝对值降序排序，选取前 n 个特征\n",
" # handle case where n might be larger than available features\n",
" n_actual = min(n, len(rankic_series))\n",
" top_features = rankic_series.sort_values(ascending=False).head(n_actual).index.tolist()\n",
" top_features = [col for col in feature_columns if col in top_features or col not in numeric_columns]\n",
" return top_features\n",
"\n",
"def create_deviation_within_dates(df, feature_columns):\n",
" groupby_col = 'cat_l2_code' # 使用 trade_date 进行分组\n",
" new_columns = {}\n",
" ret_feature_columns = feature_columns[:]\n",
"\n",
" # 自动选择所有数值型特征\n",
" num_features = [col for col in feature_columns if 'cat' not in col and 'index' not in col]\n",
"\n",
" # num_features = ['vol', 'pct_chg', 'turnover_rate', 'volume_ratio', 'cat_vol_spike', 'obv', 'maobv_6', 'return_5', 'return_10', 'return_20', 'std_return_5', 'std_return_15', 'std_return_90', 'std_return_90_2', 'act_factor1', 'act_factor2', 'act_factor3', 'act_factor4', 'act_factor5', 'act_factor6', 'rank_act_factor1', 'rank_act_factor2', 'rank_act_factor3', 'active_buy_volume_large', 'active_buy_volume_big', 'active_buy_volume_small', 'alpha_022', 'alpha_003', 'alpha_007', 'alpha_013']\n",
" num_features = [col for col in num_features if 'cat' not in col and 'industry' not in col]\n",
" num_features = [col for col in num_features if 'limit' not in col]\n",
" num_features = [col for col in num_features if 'cyq' not in col]\n",
"\n",
" # 遍历所有数值型特征\n",
" for feature in num_features:\n",
" if feature == 'trade_date': # 不需要对 'trade_date' 计算偏差\n",
" continue\n",
"\n",
" # grouped_mean = df.groupby(['trade_date'])[feature].transform('mean')\n",
" # deviation_col_name = f'deviation_mean_{feature}'\n",
" # new_columns[deviation_col_name] = df[feature] - grouped_mean\n",
" # ret_feature_columns.append(deviation_col_name)\n",
"\n",
" grouped_mean = df.groupby(['trade_date', groupby_col])[feature].transform('mean')\n",
" deviation_col_name = f'deviation_mean_{feature}'\n",
" new_columns[deviation_col_name] = df[feature] - grouped_mean\n",
" ret_feature_columns.append(deviation_col_name)\n",
"\n",
" # 将新计算的偏差特征与原始 DataFrame 合并\n",
" df = pd.concat([df, pd.DataFrame(new_columns)], axis=1)\n",
"\n",
" # for feature in ['obv', 'return_20', 'act_factor1', 'act_factor2', 'act_factor3', 'act_factor4']:\n",
" # df[f'deviation_industry_{feature}'] = df[feature] - df[f'industry_{feature}']\n",
"\n",
" return df, ret_feature_columns\n"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "47c12bb34062ae7a",
"metadata": {
"ExecuteTime": {
"end_time": "2025-04-03T14:57:50.841165Z",
"start_time": "2025-04-03T14:49:25.889057Z"
}
},
"outputs": [],
"source": [
"days = 5\n",
"validation_days = 120\n",
"\n",
"import gc\n",
"\n",
"gc.collect()\n",
"\n",
"df = df.sort_values(by=['ts_code', 'trade_date'])\n",
"df['future_return'] = df.groupby('ts_code', group_keys=False)['close'].apply(lambda x: x.shift(-days) / x - 1)\n",
"# df['future_return'] = (df.groupby('ts_code')['close'].shift(-days) - df.groupby('ts_code')['open'].shift(-1)) / \\\n",
"# df.groupby('ts_code')['open'].shift(-1)\n",
"\n",
"df['cat_up_limit'] = df['pct_chg'] > 5\n",
"df['label'] = df.groupby('ts_code')['cat_up_limit'].rolling(window=5, min_periods=1).max().groupby('ts_code').shift(-5).fillna(0).astype(int).reset_index(level=0, drop=True)\n",
"\n",
"filter_index = df['future_return'].between(df['future_return'].quantile(0.01), df['future_return'].quantile(0.99))\n",
"\n",
"# for col in [col for col in df.columns]:\n",
"# train_data[col] = train_data[col].astype('str')\n",
"# test_data[col] = test_data[col].astype('str')"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "29221dde",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"191\n"
]
}
],
"source": [
"feature_columns = [col for col in df.head(10).merge(industry_df, on=['cat_l2_code', 'trade_date'], how='left').merge(index_data, on='trade_date', how='left').columns]\n",
"feature_columns = [col for col in feature_columns if col not in ['trade_date',\n",
" 'ts_code',\n",
" 'label']]\n",
"feature_columns = [col for col in feature_columns if 'future' not in col]\n",
"feature_columns = [col for col in feature_columns if 'label' not in col]\n",
"feature_columns = [col for col in feature_columns if 'score' not in col]\n",
"feature_columns = [col for col in feature_columns if 'gen' not in col]\n",
"feature_columns = [col for col in feature_columns if 'is_st' not in col]\n",
"feature_columns = [col for col in feature_columns if 'pe_ttm' not in col]\n",
"# feature_columns = [col for col in feature_columns if 'volatility' not in col]\n",
"feature_columns = [col for col in feature_columns if 'circ_mv' not in col]\n",
"feature_columns = [col for col in feature_columns if 'code' not in col]\n",
"feature_columns = [col for col in feature_columns if col not in origin_columns]\n",
"feature_columns = [col for col in feature_columns if not col.startswith('_')]\n",
"# feature_columns = [col for col in feature_columns if col not in ['ts_code', 'trade_date', 'vol_std_5', 'cov', 'delta_cov', 'alpha_22_improved', 'alpha_007', 'consecutive_up_limit', 'mv_volatility', 'volume_growth', 'mv_growth', 'arbr']]\n",
"feature_columns = [col for col in feature_columns if col not in ['intraday_lg_flow_corr_20', \n",
" 'cap_neutral_cost_metric', \n",
" 'hurst_net_mf_vol_60', \n",
" 'complex_factor_deap_1', \n",
" 'lg_buy_consolidation_20',\n",
" 'cs_rank_ind_cap_neutral_pe',\n",
" 'cs_rank_opening_gap',\n",
" 'cs_rank_ind_adj_lg_flow']]\n",
"feature_columns = [col for col in feature_columns if col not in ['roa', 'roe']]\n",
"print(len(feature_columns))"
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "03ee5daf",
"metadata": {},
"outputs": [],
"source": [
"# df = fill_nan_with_daily_median(df, feature_columns)\n",
"for feature_col in [col for col in feature_columns if col in df.columns]:\n",
" # median_val = df[feature_col].median()\n",
" df[feature_col].fillna(0, inplace=True)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b76ea08a",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" ts_code trade_date log_circ_mv\n",
"0 000001.SZ 2019-01-02 16.574219\n",
"1 000001.SZ 2019-01-03 16.583965\n",
"2 000001.SZ 2019-01-04 16.633371\n",
"['vol', 'pct_chg', 'turnover_rate', 'volume_ratio', 'winner_rate', 'undist_profit_ps', 'ocfps', 'AR', 'BR', 'AR_BR', 'cashflow_to_ev_factor', 'book_to_price_ratio', 'turnover_rate_mean_5', 'variance_20', 'bbi_ratio_factor', 'daily_deviation', 'lg_elg_net_buy_vol', 'flow_lg_elg_intensity', 'sm_net_buy_vol', 'total_buy_vol', 'lg_elg_buy_prop', 'flow_struct_buy_change', 'lg_elg_net_buy_vol_change', 'flow_lg_elg_accel', 'chip_concentration_range', 'chip_skewness', 'floating_chip_proxy', 'cost_support_15pct_change', 'cat_winner_price_zone', 'flow_chip_consistency', 'profit_taking_vs_absorb', 'cat_is_positive', 'upside_vol', 'downside_vol', 'vol_ratio', 'return_skew', 'return_kurtosis', 'volume_change_rate', 'cat_volume_breakout', 'turnover_deviation', 'cat_turnover_spike', 'avg_volume_ratio', 'cat_volume_ratio_breakout', 'vol_spike', 'vol_std_5', 'atr_14', 'atr_6', 'obv', 'maobv_6', 'rsi_3', 'return_5', 'return_20', 'std_return_5', 'std_return_90', 'std_return_90_2', 'act_factor1', 'act_factor2', 'act_factor3', 'act_factor4', 'rank_act_factor1', 'rank_act_factor2', 'rank_act_factor3', 'cov', 'delta_cov', 'alpha_22_improved', 'alpha_003', 'alpha_007', 'alpha_013', 'vol_break', 'weight_roc5', 'smallcap_concentration', 'cost_stability', 'high_cost_break_days', 'liquidity_risk', 'turnover_std', 'mv_volatility', 'volume_growth', 'mv_growth', 'momentum_factor', 'resonance_factor', 'log_close', 'cat_vol_spike', 'up', 'down', 'obv_maobv_6', 'std_return_5_over_std_return_90', 'std_return_90_minus_std_return_90_2', 'cat_af2', 'cat_af3', 'cat_af4', 'act_factor5', 'act_factor6', 'active_buy_volume_large', 'active_buy_volume_big', 'active_buy_volume_small', 'buy_lg_vol_minus_sell_lg_vol', 'buy_elg_vol_minus_sell_elg_vol', 'ctrl_strength', 'low_cost_dev', 'asymmetry', 'lock_factor', 'cat_vol_break', 'cost_atr_adj', 'cat_golden_resonance', 'mv_turnover_ratio', 'mv_adjusted_volume', 'mv_weighted_turnover', 'nonlinear_mv_volume', 'mv_volume_ratio', 'mv_momentum', 'lg_flow_mom_corr_20_60', 'lg_flow_accel', 'profit_pressure', 'underwater_resistance', 'cost_conc_std_20', 'profit_decay_20', 'vol_amp_loss_20', 'vol_drop_profit_cnt_5', 'lg_flow_vol_interact_20', 'cost_break_confirm_cnt_5', 'atr_norm_channel_pos_14', 'turnover_diff_skew_20', 'lg_sm_flow_diverge_20', 'pullback_strong_20_20', 'vol_wgt_hist_pos_20', 'vol_adj_roc_20', 'cs_rank_net_lg_flow_val', 'cs_rank_elg_buy_ratio', 'cs_rank_rel_profit_margin', 'cs_rank_cost_breadth', 'cs_rank_dist_to_upper_cost', 'cs_rank_winner_rate', 'cs_rank_intraday_range', 'cs_rank_close_pos_in_range', 'cs_rank_pos_in_hist_range', 'cs_rank_vol_x_profit_margin', 'cs_rank_lg_flow_price_concordance', 'cs_rank_turnover_per_winner', 'cs_rank_volume_ratio', 'cs_rank_elg_buy_sell_sm_ratio', 'cs_rank_cost_dist_vol_ratio', 'cs_rank_size', 'cat_up_limit', 'industry_obv', 'industry_return_5', 'industry_return_20', 'industry__ema_5', 'industry__ema_13', 'industry__ema_20', 'industry__ema_60', 'industry_act_factor1', 'industry_act_factor2', 'industry_act_factor3', 'industry_act_factor4', 'industry_act_factor5', 'industry_act_factor6', 'industry_rank_act_factor1', 'industry_rank_act_factor2', 'industry_rank_act_factor3', 'industry_return_5_percentile', 'industry_return_20_percentile', '000852.SH_MACD', '000905.SH_MACD', '399006.SZ_MACD', '000852.SH_MACD_hist', '000905.SH_MACD_hist', '399006.SZ_MACD_hist', '000852.SH_RSI', '000905.SH_RSI', '399006.SZ_RSI', '000852.SH_Signal_line', '000905.SH_Signal_line', '399006.SZ_Signal_line', '000852.SH_amount_change_rate', '000905.SH_amount_change_rate', '399006.SZ_amount_change_rate', '000852.SH_amount_mean', '000905.SH_amount_mean', '399006.SZ_amount_mean', '000852.SH_daily_return', '000905.SH_daily_return', '399006.SZ_daily_return', '000852.SH_up_ratio_20d', '000905.SH_up_ratio_20d', '399006.SZ_up_ratio_20d', '000852.SH_volatility', '000905.SH_volatility', '399006.SZ_volatility', '000852.SH_volume_change_rate', '000905.SH_volume_change_rate', '399006.SZ_volume_change_rate']\n",
"去除极值\n",
"开始截面 MAD 去极值处理 (k=3.0)...\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"MAD Filtering: 100%|██████████| 131/131 [00:09<00:00, 13.87it/s]\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"截面 MAD 去极值处理完成。\n",
"开始截面 MAD 去极值处理 (k=3.0)...\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"MAD Filtering: 100%|██████████| 131/131 [00:20<00:00, 6.50it/s]\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"截面 MAD 去极值处理完成。\n",
"开始截面 MAD 去极值处理 (k=3.0)...\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"MAD Filtering: 0it [00:00, ?it/s]\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"截面 MAD 去极值处理完成。\n",
"开始截面 MAD 去极值处理 (k=3.0)...\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"MAD Filtering: 0it [00:00, ?it/s]\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"截面 MAD 去极值处理完成。\n",
"feature_columns: ['vol', 'pct_chg', 'turnover_rate', 'volume_ratio', 'winner_rate', 'undist_profit_ps', 'ocfps', 'AR', 'BR', 'AR_BR', 'cashflow_to_ev_factor', 'book_to_price_ratio', 'turnover_rate_mean_5', 'variance_20', 'bbi_ratio_factor', 'daily_deviation', 'lg_elg_net_buy_vol', 'flow_lg_elg_intensity', 'sm_net_buy_vol', 'total_buy_vol', 'lg_elg_buy_prop', 'flow_struct_buy_change', 'lg_elg_net_buy_vol_change', 'flow_lg_elg_accel', 'chip_concentration_range', 'chip_skewness', 'floating_chip_proxy', 'cost_support_15pct_change', 'cat_winner_price_zone', 'flow_chip_consistency', 'profit_taking_vs_absorb', 'cat_is_positive', 'upside_vol', 'downside_vol', 'vol_ratio', 'return_skew', 'return_kurtosis', 'volume_change_rate', 'cat_volume_breakout', 'turnover_deviation', 'cat_turnover_spike', 'avg_volume_ratio', 'cat_volume_ratio_breakout', 'vol_spike', 'vol_std_5', 'atr_14', 'atr_6', 'obv', 'maobv_6', 'rsi_3', 'return_5', 'return_20', 'std_return_5', 'std_return_90', 'std_return_90_2', 'act_factor1', 'act_factor2', 'act_factor3', 'act_factor4', 'rank_act_factor1', 'rank_act_factor2', 'rank_act_factor3', 'cov', 'delta_cov', 'alpha_22_improved', 'alpha_003', 'alpha_007', 'alpha_013', 'vol_break', 'weight_roc5', 'smallcap_concentration', 'cost_stability', 'high_cost_break_days', 'liquidity_risk', 'turnover_std', 'mv_volatility', 'volume_growth', 'mv_growth', 'momentum_factor', 'resonance_factor', 'log_close', 'cat_vol_spike', 'up', 'down', 'obv_maobv_6', 'std_return_5_over_std_return_90', 'std_return_90_minus_std_return_90_2', 'cat_af2', 'cat_af3', 'cat_af4', 'act_factor5', 'act_factor6', 'active_buy_volume_large', 'active_buy_volume_big', 'active_buy_volume_small', 'buy_lg_vol_minus_sell_lg_vol', 'buy_elg_vol_minus_sell_elg_vol', 'ctrl_strength', 'low_cost_dev', 'asymmetry', 'lock_factor', 'cat_vol_break', 'cost_atr_adj', 'cat_golden_resonance', 'mv_turnover_ratio', 'mv_adjusted_volume', 'mv_weighted_turnover', 'nonlinear_mv_volume', 'mv_volume_ratio', 'mv_momentum', 'lg_flow_mom_corr_20_60', 'lg_flow_accel', 'profit_pressure', 'underwater_resistance', 'cost_conc_std_20', 'profit_decay_20', 'vol_amp_loss_20', 'vol_drop_profit_cnt_5', 'lg_flow_vol_interact_20', 'cost_break_confirm_cnt_5', 'atr_norm_channel_pos_14', 'turnover_diff_skew_20', 'lg_sm_flow_diverge_20', 'pullback_strong_20_20', 'vol_wgt_hist_pos_20', 'vol_adj_roc_20', 'cs_rank_net_lg_flow_val', 'cs_rank_elg_buy_ratio', 'cs_rank_rel_profit_margin', 'cs_rank_cost_breadth', 'cs_rank_dist_to_upper_cost', 'cs_rank_winner_rate', 'cs_rank_intraday_range', 'cs_rank_close_pos_in_range', 'cs_rank_pos_in_hist_range', 'cs_rank_vol_x_profit_margin', 'cs_rank_lg_flow_price_concordance', 'cs_rank_turnover_per_winner', 'cs_rank_volume_ratio', 'cs_rank_elg_buy_sell_sm_ratio', 'cs_rank_cost_dist_vol_ratio', 'cs_rank_size', 'cat_up_limit', 'industry_obv', 'industry_return_5', 'industry_return_20', 'industry__ema_5', 'industry__ema_13', 'industry__ema_20', 'industry__ema_60', 'industry_act_factor1', 'industry_act_factor2', 'industry_act_factor3', 'industry_act_factor4', 'industry_act_factor5', 'industry_act_factor6', 'industry_rank_act_factor1', 'industry_rank_act_factor2', 'industry_rank_act_factor3', 'industry_return_5_percentile', 'industry_return_20_percentile', '000852.SH_MACD', '000905.SH_MACD', '399006.SZ_MACD', '000852.SH_MACD_hist', '000905.SH_MACD_hist', '399006.SZ_MACD_hist', '000852.SH_RSI', '000905.SH_RSI', '399006.SZ_RSI', '000852.SH_Signal_line', '000905.SH_Signal_line', '399006.SZ_Signal_line', '000852.SH_amount_change_rate', '000905.SH_amount_change_rate', '399006.SZ_amount_change_rate', '000852.SH_amount_mean', '000905.SH_amount_mean', '399006.SZ_amount_mean', '000852.SH_daily_return', '000905.SH_daily_return', '399006.SZ_daily_return', '000852.SH_up_ratio_20d', '000905.SH_up_ratio_20d', '399006.SZ_up_ratio_20d', '000852.SH_volatility', '000905.SH_volatility', '399006.SZ_volatility', '000852.SH_volume_change_rate', '000905.SH_volume_change_rate', '399006.SZ_volume_change_rate']\n",
"df最小日期: 2019-01-02\n",
"df最大日期: 2025-12-26\n",
"1345068\n",
"train_data最小日期: 2020-01-02\n",
"train_data最大日期: 2021-12-31\n",
"2940685\n",
"test_data最小日期: 2022-01-04\n",
"test_data最大日期: 2025-12-26\n",
" ts_code trade_date log_circ_mv\n",
"0 000001.SZ 2019-01-02 16.574219\n",
"1 000001.SZ 2019-01-03 16.583965\n",
"2 000001.SZ 2019-01-04 16.633371\n"
]
}
],
"source": [
"split_date = '2023-01-01'\n",
"train_data = df[filter_index & (df['trade_date'] <= split_date) & (df['trade_date'] >= '2020-01-01')]\n",
"test_data = df[(df['trade_date'] >= split_date)]\n",
"\n",
"print(df[['ts_code', 'trade_date', 'log_circ_mv']].head(3))\n",
"\n",
"industry_df = industry_df.sort_values(by=['trade_date'])\n",
"index_data = index_data.sort_values(by=['trade_date'])\n",
"\n",
"# train_data = train_data.merge(industry_df, on=['cat_l2_code', 'trade_date'], how='left')\n",
"# train_data = train_data.merge(index_data, on='trade_date', how='left')\n",
"# test_data = test_data.merge(industry_df, on=['cat_l2_code', 'trade_date'], how='left')\n",
"# test_data = test_data.merge(index_data, on='trade_date', how='left')\n",
"\n",
"train_data, test_data = train_data.replace([np.inf, -np.inf], np.nan), test_data.replace([np.inf, -np.inf], np.nan)\n",
"\n",
"# feature_columns_new = feature_columns[:]\n",
"# train_data, _ = create_deviation_within_dates(train_data, [col for col in feature_columns if col in train_data.columns])\n",
"# test_data, _ = create_deviation_within_dates(test_data, [col for col in feature_columns if col in train_data.columns])\n",
"\n",
"# feature_columns = [\n",
"# 'undist_profit_ps', \n",
"# 'AR_BR',\n",
"# 'pe_ttm',\n",
"# 'alpha_22_improved', \n",
"# 'alpha_003', \n",
"# 'alpha_007', \n",
"# 'alpha_013', \n",
"# 'cat_up_limit', \n",
"# 'cat_down_limit', \n",
"# 'up_limit_count_10d', \n",
"# 'down_limit_count_10d', \n",
"# 'consecutive_up_limit', \n",
"# 'vol_break', \n",
"# 'weight_roc5', \n",
"# 'price_cost_divergence', \n",
"# 'smallcap_concentration', \n",
"# 'cost_stability', \n",
"# 'high_cost_break_days', \n",
"# 'liquidity_risk', \n",
"# 'turnover_std', \n",
"# 'mv_volatility', \n",
"# 'volume_growth', \n",
"# 'mv_growth', \n",
"# 'lg_flow_mom_corr_20_60', \n",
"# 'lg_flow_accel', \n",
"# 'profit_pressure', \n",
"# 'underwater_resistance', \n",
"# 'cost_conc_std_20', \n",
"# 'profit_decay_20', \n",
"# 'vol_amp_loss_20', \n",
"# 'vol_drop_profit_cnt_5', \n",
"# 'lg_flow_vol_interact_20', \n",
"# 'cost_break_confirm_cnt_5', \n",
"# 'atr_norm_channel_pos_14', \n",
"# 'turnover_diff_skew_20', \n",
"# 'lg_sm_flow_diverge_20', \n",
"# 'pullback_strong_20_20', \n",
"# 'vol_wgt_hist_pos_20', \n",
"# 'vol_adj_roc_20',\n",
"# 'cashflow_to_ev_factor',\n",
"# 'ocfps',\n",
"# 'book_to_price_ratio',\n",
"# 'turnover_rate_mean_5',\n",
"# 'variance_20',\n",
"# 'bbi_ratio_factor'\n",
"# ]\n",
"# feature_columns = [col for col in feature_columns if col in train_data.columns]\n",
"# feature_columns = [col for col in feature_columns if not col.startswith('_')]\n",
"\n",
"numeric_columns = df.select_dtypes(include=['float64', 'int64']).columns\n",
"numeric_columns = [col for col in numeric_columns if col in feature_columns]\n",
"# feature_columns = select_top_features_by_rankic(df, numeric_columns, n=10)\n",
"print(feature_columns)\n",
"\n",
"# train_data = fill_nan_with_daily_median(train_data, feature_columns)\n",
"# test_data = fill_nan_with_daily_median(test_data, feature_columns)\n",
"\n",
"train_data = train_data.dropna(subset=[col for col in feature_columns if col in train_data.columns])\n",
"train_data = train_data.dropna(subset=['label'])\n",
"train_data = train_data.reset_index(drop=True)\n",
"# print(test_data.tail())\n",
"test_data = test_data.dropna(subset=[col for col in feature_columns if col in train_data.columns])\n",
"# test_data = test_data.dropna(subset=['label'])\n",
"test_data = test_data.reset_index(drop=True)\n",
"\n",
"transform_feature_columns = feature_columns\n",
"transform_feature_columns = [col for col in transform_feature_columns if col in feature_columns and not col.startswith('cat') and col in train_data.columns]\n",
"# transform_feature_columns.remove('undist_profit_ps')\n",
"print('去除极值')\n",
"cs_mad_filter(train_data, transform_feature_columns)\n",
"# print('中性化')\n",
"# cs_neutralize_market_cap_numpy(train_data, transform_feature_columns)\n",
"# print('标准化')\n",
"# cs_zscore_standardize(train_data, transform_feature_columns)\n",
"\n",
"cs_mad_filter(test_data, transform_feature_columns)\n",
"# cs_neutralize_market_cap_numpy(test_data, transform_feature_columns)\n",
"# cs_zscore_standardize(test_data, transform_feature_columns)\n",
"\n",
"mad_filter_feature_columns = [col for col in feature_columns if col not in transform_feature_columns and not col.startswith('cat') and col in train_data.columns]\n",
"cs_mad_filter(train_data, mad_filter_feature_columns)\n",
"cs_mad_filter(test_data, mad_filter_feature_columns)\n",
"\n",
"\n",
"print(f'feature_columns: {feature_columns}')\n",
"\n",
"\n",
"print(f\"df最小日期: {df['trade_date'].min().strftime('%Y-%m-%d')}\")\n",
"print(f\"df最大日期: {df['trade_date'].max().strftime('%Y-%m-%d')}\")\n",
"print(len(train_data))\n",
"print(f\"train_data最小日期: {train_data['trade_date'].min().strftime('%Y-%m-%d')}\")\n",
"print(f\"train_data最大日期: {train_data['trade_date'].max().strftime('%Y-%m-%d')}\")\n",
"print(len(test_data))\n",
"print(f\"test_data最小日期: {test_data['trade_date'].min().strftime('%Y-%m-%d')}\")\n",
"print(f\"test_data最大日期: {test_data['trade_date'].max().strftime('%Y-%m-%d')}\")\n",
"\n",
"cat_columns = [col for col in feature_columns if col.startswith('cat')]\n",
"for col in cat_columns:\n",
" train_data[col] = train_data[col].astype('category')\n",
" test_data[col] = test_data[col].astype('category')\n",
"\n",
"print(df[['ts_code', 'trade_date', 'log_circ_mv']].head(3))\n"
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "3ff2d1c5",
"metadata": {},
"outputs": [],
"source": [
"from sklearn.preprocessing import StandardScaler\n",
"from sklearn.linear_model import LogisticRegression\n",
"import matplotlib.pyplot as plt # 保持 matplotlib 导入，尽管LightGBM的绘图功能已移除\n",
"from sklearn.decomposition import PCA\n",
"import datetime # 用于日期计算\n",
"from catboost import CatBoostClassifier\n",
"from catboost import Pool\n",
"import lightgbm as lgb\n",
"\n",
"def train_model(train_data_df, feature_columns,\n",
" print_info=True, # 调整参数名，更通用\n",
" validation_days=180, use_pca=False, split_date=None,\n",
" target_column='label', type='light'): # 增加目标列参数\n",
"\n",
" print('train data size: ', len(train_data_df))\n",
" print(train_data_df[['ts_code', 'trade_date', 'log_circ_mv']])\n",
" # 确保数据按时间排序\n",
" train_data_df = train_data_df.sort_values(by='trade_date')\n",
"\n",
" # 识别数值型特征列\n",
" numeric_feature_columns = train_data_df[feature_columns].select_dtypes(include=['float64', 'int64']).columns.tolist()\n",
"\n",
" # 去除标签为空的样本\n",
" initial_len = len(train_data_df)\n",
" train_data_df = train_data_df.dropna(subset=[target_column])\n",
"\n",
" if print_info:\n",
" print(f'原始样本数: {initial_len}, 去除标签为空后样本数: {len(train_data_df)}')\n",
"\n",
" # 提取特征和标签，只取数值型特征用于线性回归\n",
" \n",
" if split_date is None:\n",
" all_dates = train_data_df['trade_date'].unique() # 获取所有唯一的 trade_date\n",
" split_date = all_dates[-validation_days] # 划分点为倒数第 validation_days 天\n",
" train_data_split = train_data_df[train_data_df['trade_date'] < split_date] # 训练集\n",
" val_data_split = train_data_df[train_data_df['trade_date'] >= split_date] # 验证集\n",
" \n",
" X_train = train_data_split[feature_columns]\n",
" y_train = train_data_split[target_column]\n",
" \n",
" X_val = val_data_split[feature_columns]\n",
" y_val = val_data_split['label']\n",
"\n",
"\n",
" # # 标准化数值特征 (使用 StandardScaler 对训练集fit并transform, 对验证集只transform)\n",
" scaler = StandardScaler()\n",
" # X_train = scaler.fit_transform(X_train)\n",
"\n",
" # 训练线性回归模型\n",
" # model = LogisticRegression(random_state=42)\n",
" \n",
" # # 使用处理后的特征和样本权重进行训练\n",
" # model.fit(X_train, y_train)\n",
"\n",
"\n",
" if type == 'cat':\n",
" params = {\n",
" 'loss_function': 'Logloss', # 适用于二分类\n",
" 'eval_metric': 'Logloss', # 评估指标\n",
" 'iterations': 1500,\n",
" 'learning_rate': 0.01,\n",
" 'depth': 10, # 控制模型复杂度\n",
" 'l2_leaf_reg': 50, # L2 正则化\n",
" 'verbose': 5000,\n",
" 'early_stopping_rounds': 300,\n",
" # 'od_type': 'Iter', # Overfitting detector type\n",
" # 'od_wait': 300, # Number of iterations to wait after the bes\n",
" 'one_hot_max_size': 50,\n",
" 'class_weights': [0.6, 1.2],\n",
" 'task_type': 'GPU',\n",
" 'has_time': True,\n",
" 'random_seed': 7\n",
" }\n",
" cat_features = [i for i, col in enumerate(feature_columns) if col.startswith('cat')]\n",
" train_pool = Pool(data=X_train, label=y_train, cat_features=cat_features)\n",
" val_pool = Pool(data=X_val, label=y_val, cat_features=cat_features)\n",
"\n",
"\n",
" model = CatBoostClassifier(**params)\n",
" model.fit(train_pool,\n",
" eval_set=val_pool, \n",
" plot=True, \n",
" use_best_model=True\n",
" )\n",
" elif type == 'light':\n",
" params = {\n",
" 'objective': 'binary',\n",
" 'metric': 'average_precision',\n",
" 'learning_rate': 0.01,\n",
" 'is_unbalance': True,\n",
" 'num_leaves': 2048,\n",
" 'min_data_in_leaf': 1024,\n",
" 'max_depth': 32,\n",
" 'max_bin': 1024,\n",
" 'feature_fraction': 0.5,\n",
" 'bagging_fraction': 0.5,\n",
" 'bagging_freq': 1,\n",
" 'lambda_l1': 50,\n",
" 'lambda_l2': 50,\n",
" 'verbosity': -1,\n",
" 'num_threads' : 8\n",
" }\n",
" categorical_feature = [col for col in feature_columns if 'cat' in col]\n",
" train_dataset = lgb.Dataset(\n",
" X_train, label=y_train,\n",
" categorical_feature=categorical_feature\n",
" )\n",
" val_dataset = lgb.Dataset(\n",
" X_val, label=y_val,\n",
" categorical_feature=categorical_feature\n",
" )\n",
"\n",
" evals = {}\n",
" callbacks = [lgb.log_evaluation(period=1000),\n",
" lgb.callback.record_evaluation(evals),\n",
" lgb.early_stopping(100, first_metric_only=True)\n",
" ]\n",
" # 训练模型\n",
" model = lgb.train(\n",
" params, train_dataset, num_boost_round=1000,\n",
" valid_sets=[train_dataset, val_dataset], valid_names=['train', 'valid'],\n",
" callbacks=callbacks\n",
" )\n",
"\n",
" # 打印特征重要性（如果需要）\n",
" if True:\n",
" lgb.plot_metric(evals)\n",
" lgb.plot_importance(model, importance_type='split', max_num_features=20)\n",
" plt.show()\n",
"\n",
"\n",
" return model, scaler, None # 返回训练好的模型、scaler 和 pca 对象"
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "c6eb5cd4-e714-420a-ac48-39af3e11ee81",
"metadata": {
"ExecuteTime": {
"end_time": "2025-04-03T15:03:18.426481Z",
"start_time": "2025-04-03T15:02:19.926352Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"train data size: 24300\n",
" ts_code trade_date log_circ_mv\n",
"0 600306.SH 2020-01-02 11.552040\n",
"1 603269.SH 2020-01-02 11.324801\n",
"2 002633.SZ 2020-01-02 11.759023\n",
"3 603991.SH 2020-01-02 11.181150\n",
"4 000691.SZ 2020-01-02 11.677910\n",
"... ... ... ...\n",
"24295 002862.SZ 2021-12-31 11.678858\n",
"24296 002494.SZ 2021-12-31 11.776090\n",
"24297 002849.SZ 2021-12-31 11.778088\n",
"24298 000632.SZ 2021-12-31 12.045706\n",
"24299 603725.SH 2021-12-31 12.027914\n",
"\n",
"[24300 rows x 3 columns]\n",
"原始样本数: 24300, 去除标签为空后样本数: 24300\n"
]
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "9fe459ca63c24ce980e8163713effbab",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"MetricVisualizer(layout=Layout(align_self='stretch', height='500px'))"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"0:\tlearn: 0.6886545\ttest: 0.6879590\tbest: 0.6879590 (0)\ttotal: 515ms\tremaining: 12m 51s\n",
"bestTest = 0.4521785262\n",
"bestIteration = 511\n",
"Shrink model to first 512 iterations.\n"
]
}
],
"source": [
"\n",
"gc.collect()\n",
"\n",
"use_pca = False\n",
"type = 'cat'\n",
"# feature_contri = [2 if feat.startswith('act_factor') or 'buy' in feat or 'sell' in feat else 1 for feat in feature_columns]\n",
"# light_params['feature_contri'] = feature_contri\n",
"# print(f'feature_contri: {feature_contri}')\n",
"model, scaler, pca = train_model(train_data\n",
" .dropna(subset=['label']).groupby('trade_date', group_keys=False)\n",
" .apply(lambda x: x.nsmallest(50, 'total_mv'))\n",
" .merge(industry_df, on=['cat_l2_code', 'trade_date'], how='left')\n",
" .merge(index_data, on='trade_date', how='left'), feature_columns, type=type)\n"
]
},
{
"cell_type": "code",
"execution_count": 20,
"id": "5d1522a7538db91b",
"metadata": {
"ExecuteTime": {
"end_time": "2025-04-03T15:04:39.656944Z",
"start_time": "2025-04-03T15:04:39.298483Z"
}
},
"outputs": [],
"source": [
"score_df = test_data.groupby('trade_date', group_keys=False).apply(lambda x: x.nsmallest(300, 'total_mv'))\n",
"# score_df = fill_nan_with_daily_median(score_df, ['pe_ttm'])\n",
"# score_df = score_df[score_df['pe_ttm'] > 0]\n",
"score_df = score_df.merge(industry_df, on=['cat_l2_code', 'trade_date'], how='left')\n",
"score_df = score_df.merge(index_data, on='trade_date', how='left')\n",
"# score_df = score_df.groupby('trade_date', group_keys=False).apply(lambda x: x.nsmallest(50, 'total_mv')).reset_index()\n",
"numeric_columns = score_df.select_dtypes(include=['float64', 'int64']).columns\n",
"numeric_columns = [col for col in feature_columns if col in numeric_columns]\n",
"# score_df.loc[:, numeric_columns] = scaler.transform(score_df[numeric_columns])\n",
"# score_df = cross_sectional_standardization(score_df, numeric_columns)\n",
"\n",
"if type == 'cat':\n",
" score_df['score'] = model.predict_proba(score_df[feature_columns])[:, 1]\n",
"elif type == 'light':\n",
" score_df['score'] = model.predict(score_df[feature_columns])\n",
"score_df['score_ranks'] = score_df.groupby('trade_date')['score'].rank(ascending=True)\n",
"\n",
"score_df = score_df.groupby('trade_date', group_keys=False).apply(\n",
" lambda x: x[x['score'] >= x['score'].quantile(0.90)] # 计算90%分位数作为阈值，筛选分数>=阈值的行\n",
").reset_index(drop=True) # drop=True 避免添加旧索引列\n",
"# save_df = score_df.groupby('trade_date', group_keys=False).apply(lambda x: x.nlargest(1, 'score')).reset_index()\n",
"save_df = score_df.groupby('trade_date', group_keys=False).apply(lambda x: x.nsmallest(2, 'total_mv')).reset_index()\n",
"save_df = save_df.sort_values(['trade_date', 'score'])\n",
"save_df[['trade_date', 'score', 'ts_code']].to_csv('predictions_test.tsv', index=False)\n"
]
},
{
"cell_type": "code",
"execution_count": 21,
"id": "c1c40917",
"metadata": {},
"outputs": [],
"source": [
"current_date = datetime.datetime.now()\n",
"\n",
"# 2. 格式化日期为字符串，例如 '2025-07-06'\n",
"# 你可以根据需要调整日期格式，例如 '%Y%m%d' 会得到 '20250706'\n",
"date_str = current_date.strftime('%Y-%m-%d')\n",
"\n",
"# 3. 构建包含日期的模型文件名\n",
"# model_filename = f'/mnt/d/PyProject/NewStock/main/train/catboost_model/catboost_model_2025-06-01.cbm'\n",
"\n",
"# model.save_model(model_filename)\n",
"# print(f\"模型已保存到: {model_filename}\")"
]
},
{
"cell_type": "code",
"execution_count": 22,
"id": "e53b209a",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"3244 1345068\n",
" ts_code trade_date turnover_rate\n",
"0 000001.SZ 2022-01-04 0.6025\n",
"1 000001.SZ 2022-01-05 1.0110\n",
"2 000001.SZ 2022-01-06 0.5709\n",
"3 000001.SZ 2022-01-07 0.5806\n",
"4 000001.SZ 2022-01-10 0.4688\n",
"... ... ... ...\n",
"2940680 605599.SH 2025-12-22 0.4745\n",
"2940681 605599.SH 2025-12-23 0.4022\n",
"2940682 605599.SH 2025-12-24 0.5145\n",
"2940683 605599.SH 2025-12-25 0.4805\n",
"2940684 605599.SH 2025-12-26 0.3122\n",
"\n",
"[2940685 rows x 3 columns]\n"
]
}
],
"source": [
"print(len(train_data[train_data['pct_chg'] > 7]), len(train_data))\n",
"print(test_data[['ts_code', 'trade_date', 'turnover_rate']])"
]
},
{
"cell_type": "code",
"execution_count": 23,
"id": "364e821a",
"metadata": {},
"outputs": [],
"source": [
"from sklearn.metrics import accuracy_score, precision_score, recall_score, f1_score, roc_auc_score, roc_curve\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"\n",
"def calculate_binary_classification_metrics(df: pd.DataFrame, score_col: str, label_col: str, future_return_col: str = None, total_mv_col: str = None, n_mv_bins: int = 10, threshold: float = 0.5):\n",
" \"\"\"\n",
" 计算二分类模型的评估指标，可选择计算 score 和 future_return 的相关性，\n",
" 并可选择计算 score 在按总市值 (total_mv) 分为 n 份后的每个分组上的预测性能（ROC AUC）。\n",
"\n",
" Args:\n",
" df: 包含 score (预测概率或置信度), label (真实二分类标签), 可选的 future_return 和 total_mv 的 Pandas DataFrame。\n",
" score_col: 包含模型预测 score 的列名。\n",
" label_col: 包含真实二分类标签 (0 或 1) 的列名。\n",
" future_return_col: (可选) 包含未来收益率的列名，用于计算相关性。\n",
" total_mv_col: (可选) 包含总市值的列名，用于按市值分 n 份分析预测性能。\n",
" n_mv_bins: (可选) 将总市值分为多少份，默认为 5。\n",
" threshold: 将 score 转换为预测类别的阈值，默认为 0.5。\n",
"\n",
" Returns:\n",
" 一个包含以下评估指标的字典：\n",
" - accuracy: 准确率\n",
" - precision: 精确率\n",
" - recall: 召回率\n",
" - f1: F1 分数\n",
" - roc_auc: ROC AUC 值\n",
" - fpr: ROC 曲线的假正率 (False Positive Rate)\n",
" - tpr: ROC 曲线的真正率 (True Positive Rate)\n",
" - thresholds: ROC 曲线的阈值\n",
" - score_return_correlation: (如果 future_return_col 提供) score 和 future_return 的皮尔逊相关系数\n",
" - mv_roc_auc: (如果 total_mv_col 提供) 一个字典，包含按总市值分为 n 份后的每个市值分组对应的 ROC AUC 值\n",
" \"\"\"\n",
" y_true = df[label_col].values\n",
" y_score = df[score_col].values\n",
" y_pred = (y_score >= threshold).astype(int)\n",
"\n",
" metrics = {}\n",
" metrics['accuracy'] = accuracy_score(y_true, y_pred)\n",
" metrics['precision'] = precision_score(y_true, y_pred)\n",
" metrics['recall'] = recall_score(y_true, y_pred)\n",
" metrics['f1'] = f1_score(y_true, y_pred)\n",
" metrics['roc_auc'] = roc_auc_score(y_true, y_score)\n",
" metrics['fpr'], metrics['tpr'], metrics['thresholds'] = roc_curve(y_true, y_score)\n",
"\n",
" if future_return_col in df.columns:\n",
" metrics['score_return_correlation'] = df[score_col].corr(df[future_return_col])\n",
"\n",
" if total_mv_col in df.columns and n_mv_bins > 1:\n",
" metrics['mv_roc_auc'] = {}\n",
" df['mv_quantile'] = pd.cut(df[total_mv_col], bins=n_mv_bins, labels=False, duplicates='drop')\n",
" for i in range(df['mv_quantile'].nunique()):\n",
" mv_group = df[df['mv_quantile'] == i]\n",
" if len(mv_group) > 0 and len(np.unique(mv_group[label_col])) > 1 and len(np.unique(mv_group[score_col])) > 1:\n",
" roc_auc_mv = roc_auc_score(mv_group[label_col], mv_group[score_col])\n",
" lower_bound = df[total_mv_col][df['mv_quantile'] == i].min()\n",
" upper_bound = df[total_mv_col][df['mv_quantile'] == i].max()\n",
" metrics['mv_roc_auc'][f'{lower_bound:.0e}-{upper_bound:.0e}'] = roc_auc_mv\n",
" else:\n",
" lower_bound = df[total_mv_col][df['mv_quantile'] == i].min()\n",
" upper_bound = df[total_mv_col][df['mv_quantile'] == i].max()\n",
" metrics['mv_roc_auc'][f'{lower_bound:.0e}-{upper_bound:.0e}'] = np.nan\n",
" print(f'{lower_bound:.0e}-{upper_bound:.0e}')\n",
" df.drop(columns=['mv_quantile'], inplace=True)\n",
"\n",
" return metrics\n",
"\n",
"def plot_roc_curve(metrics: dict):\n",
" plt.figure(figsize=(8, 6))\n",
" plt.plot(metrics['fpr'], metrics['tpr'], color='darkorange', lw=2, label=f'ROC curve (AUC = {metrics[\"roc_auc\"]:.2f})')\n",
" plt.plot([0, 1], [0, 1], color='navy', lw=2, linestyle='--')\n",
" plt.xlabel('False Positive Rate')\n",
" plt.ylabel('True Positive Rate')\n",
" plt.title('Receiver Operating Characteristic (ROC)')\n",
" plt.legend(loc=\"lower right\")\n",
" plt.grid(True)\n",
" plt.show()\n"
]
},
{
"cell_type": "code",
"execution_count": 24,
"id": "1f6e6336",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"6e+04-9e+04\n",
"9e+04-1e+05\n",
"1e+05-1e+05\n",
"1e+05-2e+05\n",
"2e+05-2e+05\n",
"2e+05-2e+05\n",
"2e+05-2e+05\n",
"2e+05-3e+05\n",
"3e+05-3e+05\n",
"3e+05-3e+05\n",
"二分类评估指标:\n",
"accuracy: 0.6607\n",
"precision: 0.4707\n",
"recall: 0.2058\n",
"f1: 0.2864\n",
"roc_auc: 0.6113\n",
"fpr: (array of length 12329)\n",
"tpr: (array of length 12329)\n",
"thresholds: (array of length 12329)\n",
"score_return_correlation: -0.0256\n",
"mv_roc_auc: {'6e+04-9e+04': 0.5764925373134329, '9e+04-1e+05': 0.5824607329842932, '1e+05-1e+05': 0.5853884314061875, '1e+05-2e+05': 0.5701944637861504, '2e+05-2e+05': 0.611481089969462, '2e+05-3e+05': 0.6152227912178468, '3e+05-3e+05': 0.6330766655507197}\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAArMAAAIjCAYAAAAQgZNYAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAoadJREFUeJzs3XVYVOn7BvB7gKFTRTBQ7O5au0Vdu7Ax19YVe/1au+a6dqwttigWdnf32o2NgkHHMPP+/vDn6DiDMghzmOH+XNde63lO3cxh4OHMOe+RCSEEiIiIiIiMkJnUAYiIiIiIkovNLBEREREZLTazRERERGS02MwSERERkdFiM0tERERERovNLBEREREZLTazRERERGS02MwSERERkdFiM0tERERERovNLJEJ8vT0RJcuXaSOke7UqFEDNWrUkDrGD40fPx4ymQyhoaFSR0lzZDIZxo8fnyLbCgoKgkwmg5+fX4psDwAuXLgAS0tLPH36NMW2mdLatm2LNm3aSB2D0hE2s0R68vPzg0wmU/9nYWGBbNmyoUuXLnj58qXU8dK0qKgo/PXXXyhevDhsbW3h5OSEqlWrYvXq1TCWJ2vfvn0b48ePR1BQkNRRtCiVSqxcuRI1atRAhgwZYGVlBU9PT3Tt2hWXLl2SOl6KWL9+PWbPni11DA2GzDR69Gi0a9cOOXPmVNdq1Kih8TPJxsYGxYsXx+zZs6FSqXRu5927dxg2bBgKFCgAa2trZMiQAV5eXti1a1ei+w4PD8eECRNQokQJ2Nvbw8bGBkWLFsWIESPw6tUr9XIjRozAli1bcP369ZT7wom+QyaM5TcIURrh5+eHrl274s8//0SuXLkQGxuLc+fOwc/PD56enrh58yasra0lzRgXFwczMzPI5XJJc3ztzZs3qF27Nu7cuYO2bduievXqiI2NxZYtW3DixAl4e3tj3bp1MDc3lzrqdwUEBKB169Y4evSo1lnY+Ph4AIClpaXBc8XExKBFixbYt28fqlWrhsaNGyNDhgwICgrCpk2bcP/+fTx79gzZs2fH+PHjMWHCBISEhCBTpkwGz/ozGjVqhJs3b6baHxOxsbGwsLCAhYXFT2cSQiAuLg5yuTxFvq+vXbuGUqVK4cyZM6hYsaK6XqNGDTx69AhTpkwBAISGhmL9+vW4ePEi/vjjD0yaNEljO/fu3UPt2rUREhKCrl27omzZsvj48SPWrVuHa9euYejQoZg+fbrGOo8fP0adOnXw7NkztG7dGlWqVIGlpSX+++8/bNiwARkyZMD9+/fVy1eoUAEFChTA6tWrf/rrJvohQUR6WblypQAgLl68qFEfMWKEACD8/f0lSiatmJgYoVQqE53v5eUlzMzMxI4dO7TmDR06VAAQU6dOTc2IOkVGRuq1/ObNmwUAcfTo0dQJlEz9+vUTAMSsWbO05iUkJIjp06eL58+fCyGEGDdunAAgQkJCUi2PSqUS0dHRKb7dX3/9VeTMmTNFt6lUKkVMTEyy10+NTLoMHDhQ5MiRQ6hUKo169erVRZEiRTRqMTExImfOnMLBwUEkJCSo6/Hx8aJo0aLC1tZWnDt3TmOdhIQE4e3tLQCIjRs3qusKhUKUKFFC2NraipMnT2rlCgsLE3/88YdG7Z9//hF2dnYiIiIi2V8vUVKxmSXSU2LN7K5duwQAMXnyZI36nTt3RMuWLYWLi4uwsrISZcqU0dnQffjwQfz+++8iZ86cwtLSUmTLlk106tRJo+GIjY0VY8eOFXny5BGWlpYie/bsYtiwYSI2NlZjWzlz5hQ+Pj5CCCEuXrwoAAg/Pz+tfe7bt08AEDt37lTXXrx4Ibp27SoyZ84sLC0tReHChcXy5cs11jt69KgAIDZs2CBGjx4tsmbNKmQymfjw4YPO1+zs2bMCgOjWrZvO+QqFQuTLl0+4uLioG6AnT54IAGL69Oli5syZIkeOHMLa2lpUq1ZN3LhxQ2sbSXmdPx+7Y8eOiT59+ghXV1fh7OwshBAiKChI9OnTR+TPn19YW1uLDBkyiFatWoknT55orf/tf58b2+rVq4vq1atrvU7+/v5i4sSJIlu2bMLKykrUqlVLPHjwQOtrmD9/vsiVK5ewtrYW5cqVEydOnNDapi7Pnz8XFhYWom7dut9d7rPPzeyDBw+Ej4+PcHJyEo6OjqJLly4iKipKY9kVK1aImjVrCldXV2FpaSkKFSokFi5cqLXNnDlzil9//VXs27dPlClTRlhZWakb66RuQwgh9uzZI6pVqybs7e2Fg4ODKFu2rFi3bp0Q4tPr++1r/3UTmdT3BwDRr18/sXbtWlG4cGFhYWEhtm3bpp43btw49bLh4eFi0KBB6velq6urqFOnjrh8+fIPM33+Hl65cqXG/u/cuSNat24tMmXKJKytrUX+/Pm1mkFdcuTIIbp06aJV19XMCiFEq1atBADx6tUrdW3Dhg0CgPjzzz917uPjx4/C2dlZFCxYUF3buHGjACAmTZr0w4yfXb9+XQAQW7duTfI6RMmV9M9RiOi7Pn/E6OLioq7dunULlStXRrZs2TBy5EjY2dlh06ZNaNasGbZs2YLmzZsDACIjI1G1alXcuXMH3bp1Q+nSpREaGorAwEC8ePECmTJlgkqlQpMmTXDq1Cn89ttvKFSoEG7cuIFZs2bh/v372L59u85cZcuWRe7cubFp0yb4+PhozPP394eLiwu8vLwAfLoU4JdffoFMJkP//v3h6uqKvXv3onv37ggPD8fvv/+usf5ff/0FS0tLDB06FHFxcYl+vL5z504AQOfOnXXOt7CwQPv27TFhwgScPn0aderUUc9bvXo1IiIi0K9fP8TGxmLOnDmoVasWbty4ATc3N71e58/69u0LV1dXjB07FlFRUQCAixcv4syZM2jbti2yZ8+OoKAg/Pvvv6hRowZu374NW1tbVKtWDQMHDsTcuXPxxx9/oFChQgCg/n9ipk6dCjMzMwwdOhRhYWH4+++/0aFDB5w/f169zL///ov+/fujatWqGDx4MIKCgtCsWTO4uLgge/bs393+3r17kZCQgE6dOn13uW+1adMGuXLlwpQpU3DlyhUsW7YMmTNnxrRp0zRyFSlSBE2aNIGFhQV27tyJvn37QqVSoV+/fhrbu3fvHtq1a4devXqhZ8+eKFCggF7b8PPzQ7du3VCkSBGMGjUKzs7OuHr1Kvbt24f27dtj9OjRCAsLw4sXLzBr1iwAgL29PQDo/f44cuQINm3ahP79+yNTpkzw9PTU+Rr17t0bAQEB6N+/PwoXLox3797h1KlTuHPnDkqXLv3dTLr8999/qFq1KuRyOX777Td4enri0aNH2Llzp9blAF97+fIlnj17htKlSye6zLc+34Dm7Oysrv3ovejk5ISmTZti1apVePjwIfLmzYvAwEAA0Ov7q3DhwrCxscHp06e13n9EKU7qbprI2Hw+O3fo0CEREhIinj9/LgICAoSrq6uwsrJSf5QrhBC1a9cWxYoV0zgzpFKpRKVKlUS+fPnUtbFjxyZ6FuPzR4pr1qwRZmZmWh/zLVq0SAAQp0+fVte+PjMrhBCjRo0ScrlcvH//Xl2Li4sTzs7OGmdLu3fvLrJkySJCQ0M19tG2bVvh5OSkPmv6+Yxj7ty5k/RRcrNmzQSARM/cCiHE1q1bBQAxd+5cIcSXs1o2NjbixYsX6uXOnz8vAIjBgwera0l9nT8fuypVqmh89CqE0Pl1fD6jvHr1anXte5cZJHZmtlChQiIuLk5dnzNnjgCgPsMcFxcnMmbMKMqVKycUCoV6OT8/PwHgh2dmBw8eLACIq1evfne5zz6fmf32THnz5s1FxowZNWq6XhcvLy+RO3dujVrOnDkFALFv3z6t5ZOyjY8fPwoHBwdRoUIFrY/8v/5YPbGP9PV5fwAQZmZm4tatW1rbwTdnZp2cnES/fv20lvtaYpl0nZmtVq2acHBwEE+fPk30a9Tl0KFDWp+ifFa9enVRsGBBERISIkJCQsTdu3fFsGHDBADx66+/aixbsmRJ4eTk9N19zZw5UwAQgYGBQgghSpUq9cN1dMmfP79o0KCB3usR6YujGRAlU506deDq6goPDw+0atUKdnZ2CAwMVJ9Fe//+PY4cOYI2bdogIiICoaGhCA0Nxbt37+Dl5YUHDx6oRz/YsmULSpQoofMMhkwmAwBs3rwZhQoVQsGCBdXbCg0NRa1atQAAR48eTTSrt7c3FAoFtm7dqq4dOHAAHz9+hLe3N4BPN6ts2bIFjRs3hhBCYx9eXl4ICwvDlStXNLbr4+MDGxubH75WERERAAAHB4dEl/k8Lzw8XKPerFkzZMuWTT1dvnx5VKhQAXv27AGg3+v8Wc+ePbVuyPn661AoFHj37h3y5s0LZ2dnra9bX127dtU4a121alUAn26qAYBLly7h3bt36Nmzp8aNRx06dNA405+Yz6/Z915fXXr37q0xXbVqVbx7907jGHz9uoSFhSE0NBTVq1fH48ePERYWprF+rly51Gf5v5aUbRw8eBAREREYOXKk1g2Un98D36Pv+6N69eooXLjwD7fr7OyM8+fPa9ytn1whISE4ceIEunXrhhw5cmjM+9HX+O7dOwBI9Pvh7t27cHV1haurKwoWLIjp06ejSZMmWsOCRURE/PD75Nv3Ynh4uN7fW5+zcvg3MgReZkCUTAsWLED+/PkRFhaGFStW4MSJE7CyslLPf/jwIYQQGDNmDMaMGaNzG2/fvkW2bNnw6NEjtGzZ8rv7e/DgAe7cuQNXV9dEt5WYEiVKoGDBgvD390f37t0BfLrEIFOmTOpf9iEhIfj48SOWLFmCJUuWJGkfuXLl+m7mzz7/IoyIiND4yPNriTW8+fLl01o2f/782LRpEwD9Xufv5Y6JicGUKVOwcuVKvHz5UmOosG+bNn1927h8bkg+fPgAAOoxQ/PmzauxnIWFRaIff3/N0dERwJfXMCVyfd7m6dOnMW7cOJw9exbR0dEay4eFhcHJyUk9ndj3Q1K28ejRIwBA0aJF9foaPtP3/ZHU792///4bPj4+8PDwQJkyZdCwYUN07twZuXPn1jvj5z9ekvs1Akh0CDtPT08sXboUKpUKjx49wqRJkxASEqL1h4GDg8MPG8xv34uOjo7q7PpmTcofIkQ/i80sUTKVL18eZcuWBfDp7GGVKlXQvn173Lt3D/b29urxHYcOHarzbBWg3bx8j0qlQrFixTBz5kyd8z08PL67vre3NyZNmoTQ0FA4ODggMDAQ7dq1U58J/Jy3Y8eOWtfWfla8eHGN6aSclQU+XVO6fft2/Pfff6hWrZrOZf777z8ASNLZsq8l53XWlXvAgAFYuXIlfv/9d1SsWBFOTk6QyWRo27ZtomN1JlViwzIl1pjoq2DBggCAGzduoGTJkkle70e5Hj16hNq1a6NgwYKYOXMmPDw8YGlpiT179mDWrFlar4uu11XfbSSXvu+PpH7vtmnTBlWrVsW2bdtw4MABTJ8+HdOmTcPWrVvRoEGDn86dVBkzZgTw5Q+gb9nZ2Wlca165cmWULl0af/zxB+bOnauuFypUCNeuXcOzZ8+0/pj57Nv3YsGCBXH16lU8f/78hz9nvvbhwwedf4wSpTQ2s0QpwNzcHFOmTEHNmjUxf/58jBw5Un3mRi6Xa/yS0SVPnjy4efPmD5e5fv06ateunayzHd7e3pgwYQK2bNkCNzc3hIeHo23btur5rq6ucHBwgFKp/GFefTVq1AhTpkzB6tWrdTazSqUS69evh4uLCypXrqwx78GDB1rL379/X33GUp/X+XsCAgLg4+ODGTNmqGuxsbH4+PGjxnKpcabp8wD4Dx8+RM2aNdX1hIQEBAUFaf0R8a0GDRrA3Nwca9eu1fsmsO/ZuXMn4uLiEBgYqNH4fO+SluRuI0+ePACAmzdvfvePvMRe/599f3xPlixZ0LdvX/Tt2xdv375F6dKlMWnSJHUzm9T9ff5e/dF7XZfPf7A8efIkScsXL14cHTt2xOLFizF06FD1a9+oUSNs2LABq1evxv/+9z+t9cLDw7Fjxw4ULFhQfRwaN26MDRs2YO3atRg1alSS9p+QkIDnz5+jSZMmSVqe6GfwmlmiFFKjRg2UL18es2fPRmxsLDJnzowaNWpg8eLFeP36tdbyISEh6n+3bNkS169fx7Zt27SW+3yWrE2bNnj58iWWLl2qtUxMTIz6rvzEFCpUCMWKFYO/vz/8/f2RJUsWjcbS3NwcLVu2xJYtW3T+sv06r74qVaqEOnXqYOXKlTqfMDR69Gjcv38fw4cP1zpjtn37do1rXi9cuIDz58+rGwl9XufvMTc31zpTOm/ePCiVSo2anZ0dAGg1uT+jbNmyyJgxI5YuXYqEhAR1fd26dYmeifuah4cHevbsiQMHDmDevHla81UqFWbMmIEXL17olevzmdtvL7lYuXJlim+jXr16cHBwwJQpUxAbG6sx7+t17ezsdF728bPvD12USqXWvjJnzoysWbMiLi7uh5m+5erqimrVqmHFihV49uyZxrwfnaXPli0bPDw89HqS2/Dhw6FQKDTOVrdq1QqFCxfG1KlTtbalUqnQp08ffPjwAePGjdNYp1ixYpg0aRLOnj2rtZ+IiAiMHj1ao3b79m3ExsaiUqVKSc5LlFw8M0uUgoYNG4bWrVvDz88PvXv3xoIFC1ClShUUK1YMPXv2RO7cufHmzRucPXsWL168UD/ucdiwYeonS3Xr1g1lypTB+/fvERgYiEWLFqFEiRLo1KkTNm3ahN69e+Po0aOoXLkylEol7t69i02bNmH//v3qyx4S4+3tjbFjx8La2hrdu3eHmZnm37NTp07F0aNHUaFCBfTs2ROFCxfG+/fvceXKFRw6dAjv379P9muzevVq1K5dG02bNkX79u1RtWpVxMXFYevWrTh27Bi8vb0xbNgwrfXy5s2LKlWqoE+fPoiLi8Ps2bORMWNGDB8+XL1MUl/n72nUqBHWrFkDJycnFC5cGGfPnsWhQ4fUH+9+VrJkSZibm2PatGkICwuDlZUVatWqhcyZMyf7tbG0tMT48eMxYMAA1KpVC23atEFQUBD8/PyQJ0+eJJ35mzFjBh49eoSBAwdi69ataNSoEVxcXPDs2TNs3rwZd+/e1TgTnxT16tWDpaUlGjdujF69eiEyMhJLly5F5syZdf7h8DPbcHR0xKxZs9CjRw+UK1cO7du3h4uLC65fv47o6GisWrUKAFCmTBn4+/vD19cX5cqVg729PRo3bpwi749vRUREIHv27GjVqpX6Ea6HDh3CxYsXNc7gJ5ZJl7lz56JKlSooXbo0fvvtN+TKlQtBQUHYvXs3rl279t08TZs2xbZt25J8LWrhwoXRsGFDLFu2DGPGjEHGjBlhaWmJgIAA1K5dG1WqVNF4Atj69etx5coVDBkyRON7RS6XY+vWrahTpw6qVauGNm3aoHLlypDL5bh165b6U5WvhxY7ePAgbG1tUbdu3R/mJPpphh9Agci4JfbQBCE+PUkoT548Ik+ePOqhnx49eiQ6d+4s3N3dhVwuF9myZRONGjUSAQEBGuu+e/dO9O/fX2TLlk094LuPj4/GMFnx8fFi2rRpokiRIsLKykq4uLiIMmXKiAkTJoiwsDD1ct8OzfXZgwcP1AO7nzp1SufX9+bNG9GvXz/h4eEh5HK5cHd3F7Vr1xZLlixRL/N5yKnNmzfr9dpFRESI8ePHiyJFiggbGxvh4OAgKleuLPz8/LSGJvr6oQkzZswQHh4ewsrKSlStWlVcv35da9tJeZ2/d+w+fPggunbtKjJlyiTs7e2Fl5eXuHv3rs7XcunSpSJ37tzC3Nw8SQ9N+PZ1Smww/blz54qcOXMKKysrUb58eXH69GlRpkwZUb9+/SS8up+e4LRs2TJRtWpV4eTkJORyuciZM6fo2rWrxrBdiT0B7PPr8/WDIgIDA0Xx4sWFtbW18PT0FNOmTRMrVqzQWu7zQxN0Seo2Pi9bqVIlYWNjIxwdHUX58uXFhg0b1PMjIyNF+/bthbOzs9ZDE5L6/sD/PzRBF3w1NFdcXJwYNmyYKFGihHBwcBB2dnaiRIkSWg98SCxTYsf55s2bonnz5sLZ2VlYW1uLAgUKiDFjxujM87UrV64IAFrDjyX20AQhhDh27JjWcGNCCPH27Vvh6+sr8ubNK6ysrISzs7OoU6eOejguXT58+CDGjh0rihUrJmxtbYW1tbUoWrSoGDVqlHj9+rXGshUqVBAdO3b84ddElBJkQqTQHQhERCkoKCgIuXLlwvTp0zF06FCp40hCpVLB1dUVLVq00PnxOaU/tWvXRtasWbFmzRqpoyTq2rVrKF26NK5cuaLXDYlEycVrZomI0oDY2Fit6yZXr16N9+/fo0aNGtKEojRn8uTJ8Pf3Vw/nlhZNnToVrVq1YiNLBsNrZomI0oBz585h8ODBaN26NTJmzIgrV65g+fLlKFq0KFq3bi11PEojKlSogPj4eKljfNfGjRuljkDpDJtZIqI0wNPTEx4eHpg7dy7ev3+PDBkyoHPnzpg6darG08OIiEgTr5klIiIiIqPFa2aJiIiIyGixmSUiIiIio5XurplVqVR49eoVHBwcUuWxlERERET0c4QQiIiIQNasWbUe8POtdNfMvnr1Ch4eHlLHICIiIqIfeP78ObJnz/7dZdJdM+vg4ADg04vj6OiY6vtTKBQ4cOAA6tWrB7lcnur7o5THY2j8eAyNH4+hcePxM36GPobh4eHw8PBQ923fk+6a2c+XFjg6OhqsmbW1tYWjoyPfwEaKx9D48RgaPx5D48bjZ/ykOoZJuSSUN4ARERERkdFiM0tERERERovNLBEREREZLTazRERERGS02MwSERERkdFiM0tERERERovNLBEREREZLTazRERERGS02MwSERERkdFiM0tERERERovNLBEREREZLTazRERERGS02MwSERERkdFiM0tERERERkvSZvbEiRNo3LgxsmbNCplMhu3bt/9wnWPHjqF06dKwsrJC3rx54efnl+o5iYiIiChtkrSZjYqKQokSJbBgwYIkLf/kyRP8+uuvqFmzJq5du4bff/8dPXr0wP79+1M5KRERERGlRRZS7rxBgwZo0KBBkpdftGgRcuXKhRkzZgAAChUqhFOnTmHWrFnw8vJKrZhERERE6U9CLPDuNnDxb8jv+aOc3S+QPYwFCnlLnUyDpM2svs6ePYs6depo1Ly8vPD7778nuk5cXBzi4uLU0+Hh4QAAhUIBhUKRKjm/9nkfhtgXpQ4eQ+PHY2j8eAyNG4+fERECZtcXwPyELwDgYWgG9ApohKWtXZAb5xD/9gYUeVukegx9vleMqpkNDg6Gm5ubRs3NzQ3h4eGIiYmBjY2N1jpTpkzBhAkTtOoHDhyAra1tqmX91sGDBw22L0odPIbGj8fQ+PEYGjcev7TLJfYeqr0YoVHbdK0Iemxugog4K7Rd2wqn+q3AxVeWCN2zJ9XzREdHJ3lZo2pmk2PUqFHw9fVVT4eHh8PDwwP16tWDo6Njqu9foVDg4MGDqFu3LuRyearvj1Iej6Hx4zE0fjyGxo3HLw2KeQezh1tgdmMJZKH/ac5SWGDwjvpYfK6suvZBZEOAwyw0a9rdIMfw8yfpSWFUzay7uzvevHmjUXvz5g0cHR11npUFACsrK1hZWWnV5XK5Qd9Qht4fpTweQ+PHY2j8eAyNG4+fhD4+Bo4OBF6fB2JCE13s3tuMaLOmNf577a6utW9fDPPm1cPJk4cNdgz12YdRNbMVK1bEnm9ObR88eBAVK1aUKBERERFRGnZnA7CnfZIWXXe9PHptaYSoaBUAwNraAvPnN0C3bqWQkJCQmil/iqTNbGRkJB4+fKiefvLkCa5du4YMGTIgR44cGDVqFF6+fInVq1cDAHr37o358+dj+PDh6NatG44cOYJNmzZh9+7dUn0JRERERGmDEEB4EHBtIfBwG/Dx0feXd8oNeNZDdN4eGDjxJZavuQrgUyNbsGAmbN7cGkWLZk712D9L0mb20qVLqFmzpnr687WtPj4+8PPzw+vXr/Hs2TP1/Fy5cmH37t0YPHgw5syZg+zZs2PZsmUclouIiIjSp5j3wJlxwIvjQOiN7y9rbgmUHABUGg9Y2qvL548+wfLlu9TTPj4lsGBBQ9jZWaZS6JQlaTNbo0YNCCESna/r6V41atTA1atXUzEVERERURomBPD8GLC5VtKWz98GqLMQsMmoc3bNmrkwYkRlzJt3AQsXNoSPT8mUSmoQRnXNLBEREVG6o1ICN5YBd9YBZuafGtnEmFt9unwgbzOgRB/A0UNrkZgYBaytLSCTydS1v/6qie7dSyFfPt0Nb1rGZpaIiIgoLYp5D6wtDYQ//fGy5UcCv4wF5LpHd/rsxo03aNMmAAMGlEffvuXUdbnc3CgbWYDNLBEREVHaEv4MWJrzx8vlqA0U6wEU8Aa+OsuqixACy5ZdwcCB+xAbm4DBg/ejYsXsKFUqSwqFlg6bWSIiIqK04mEgsKOp7nnZqwNVpwIZCwNWSX/wU0REHHr12oUNG26qa4UKZYK9vXHc4PUjbGaJiIiIpBb+HFiaQ/e8yn8BFUb/8OyrLlevvkabNgF4+PC9uta3b1nMmOEFa2vTaANN46sgIiIiMjZvrgL+VQFFlO75pQYAteYma9NCCPz77yX4+u5HXJwSAODoaIVlyxqjdesiyU2cJrGZJSIiIjKU6FDgwhTg8szvL9cnBLDNlKxdhIXFokePnQgIuK2ulSmTBf7+rZAnT4ZkbTMtYzNLRERElFoS4oAHW4Cr84HXZ7+/rNweqO8H5G/5U7sUArh06ZV6euDA8vj777qwsjLNts80vyoiIiIiKb04BZwYCrw+/+NlWx4APOum2K6dna3h798KjRtvwOLFjdCsWcEU23ZaxGaWiIiIKKUEHQS2NwKU8YkvY50RqDgGKNYTkNv+9C4/fIhBXJwS7u5fHlFbvnw2PHkyCLa28p/eflrHZpaIiIjoZ9zdCJz+H/DxUeLL5KwLVJkEuJdLfJlkOHfuBdq2DYCnpzMOHeoMCwsz9bz00MgCbGaJiIiI9KOIBq7MAa4vAiKeJb6cY06g7mIgZ71kDav1PSqVwMyZZzFq1GEkJKjw9GkYpk07hdGjq6XofowBm1kiIiKipIiPBBZkAFSK7y/nUQOoOg3IUj5VYoSGRqNLl+3YvfuBula5sgc6dy6RKvtL69jMEhERESVGCODVWWBj5R8v63MTyJS6Y7ieOvUM7dptwYsX4erayJGV8eefNSGXm6fqvtMqNrNERERE34r9+Ok62GsLEl+maPdP18HauaV6HJVKYNq0Uxgz5iiUSgEAyJTJFmvWNEf9+nlTff9pGZtZIiIios8iXwO7vIGXJxNfpvxIoOoUg0WKj1eiSZMN2L//yw1m1avnxPr1LZE1q4PBcqRVbGaJiIiIbq0C9nVJfH7WSkCdRYBrMYNF+szS0hy5cjkD+HQf2f/+Vw1jx1bXGLkgPWMzS0REROmTUgFsbwwE7U98mcp/ARX+AGTSNo6zZtXHkycfMXRoJdSpk1vSLGkNm1kiIiJKXx7tBLY3+f4yTXcAeX+wTCoJDo7Ef/+9Qb16edQ1a2sL7NvXUZI8aR2bWSIiIkofLkwDTo5MfL6lA9DpGuAs3ZnPQ4ceo2PHrYiMjMelS7+hYMFMkmUxFmxmiYiIyHQlxAF31gIHeiS+TPUZQJnBKf5gA30kJKgwYcIxTJp0EuLTYAX4/fd9PBubBGxmiYiIyLQoFcA9f2Bvp+8v1+UWkLGwYTJ9x8uX4WjffitOnHiqrtWvnxerVzeTLpQRYTNLREREpuHNVWB/FyDkv+8v99tzwCG7QSL9yL59D9Gp0zaEhkYDAMzNZZg0qRaGDasMMzPpzhQbEzazREREZNyCLwHryn1/mdK/AyX7AS5p4wEDCoUSY8YcxbRpp9W17NkdsXFjS1SunEPCZMaHzSwREREZHyE+XQt7dgLw8ZHuZTy9gJpzgAwFDJstCdq334qAgNvq6UaN8sPPrykyZrSVMJVxYjNLRERExkMIYH/XTw85SEyTLUC+FobLlAx9+5bF1q13YGYmw9SpteHrWxEyCW9AM2ZsZomIiMg4RL4GFmfVPc/SEWhzDHArZdBIyVWzZi7MmVMfZctmxS+/pI3rd40Vm1kiIiJKu4QAXp4C/Kvpnp+9GlB3KZAhv2Fz6SEo6CMWLbqEyZNra9zU1b9/eQlTmQ42s0RERJQ2XfgbODlC9zz7rEDPZ4CZuWEz6Wnbtjvo1i0QHz/GImNGGwwbVlnqSCZH2gcNExEREX0t8jWwuiQwQ5Z4I9twPdDrZZpuZOPiEjBw4F60aLEJHz/GAgCWL7+KuLgEiZOZHp6ZJSIiImk9DAQuTgNenQUgdC/jkAPoeAmwdTVotOR49Og9vL0DcPnya3WtdevCWLq0Mays2HqlNL6iREREZHhCBdxZ/+OndGWpCLQ9mabPwn5t8+Zb6NFjJ8LD4wAAVlbmmDXLC717l+VoBamEzSwREREZ1vVFwKE+ic+3zwpU+gso2gWQGccVkbGxCfD13Y9//72kruXLlwGbNrVGyZLuEiYzfWxmiYiIKNU5xT2GhV8BIPyJ7gUyFgHqLAKyVQaM8AzmpEknNBrZ9u2LYdGiX+HgYCVhqvSBzSwRERGlDpUSODsB8nN/oUZiy+RpAtScDTjlMlyuVDB8eGVs2nQbz56FYd68BujevRQvKzAQNrNERESUsoQAjvkCV2YnvoyVE9DtIWCbyWCxUpODgxUCAloDAIoVc5M4TfrCZpaIiIhShhDAhSnAqdG6Z2cqAVnNWUCOmgYOlrLu3AlBr167sHp1c3h6OqvrbGKlwWaWiIiIft69TcAub52zlGWGYNfHSmj4a2PI5XIDB0tZq1ZdQ9++exAdrYC3dwBOnuwKS0vjGGnBVLGZJSIiouR7fR5Y/0vi832VUCUogT17DJcpFURFxaNfvz1Yteq6uhYdrUBISBSyZXOUMBmxmSUiIiL9RbwAVhUF4sK05+VqANRf9dUDDpQGjZbSbtx4gzZtAnD3bqi61qNHKcyZ0wC2tsZ9ptkUsJklIiKipPvRmdiB0YDcxnB5UpEQAsuXX8WAAXsRG/vpMbT29pZYvLgR2rcvJnE6+ozNLBEREf1YdCiwsxXw4rju+R0vAW5lDJspFUVExKF3791Yv/6GulaihBs2bWqN/PkzSpiMvsVmloiIiBInVMDB3sCNpdrzMhUFqv8DeHoZPlcqO3v2hUYj27t3GcyaVR/W1myd0hoeESIiItKmUgJb6gHPjuie33wXkPtXw2YyoHr18mDIkIpYsuQyli1rgjZtikgdiRLBZpaIiIg0HR8GXPpH9zyvFUCRLkb5yNnviYqKh62tXOOpXZMn10a/fuWQK5eLhMnoR8ykDkBERERpxINtwIqCuhvZgu0BXyVQtKvJNbKXLr1C8eKLsGTJZY26paU5G1kjwDOzRERE6d3VBcDJkYAiUnte0W5A3SWAmek9GEAIgXnzLmDo0ANQKFQYNGgffvklO0qUcJc6GumBzSwREVF69eEhsCKf7nkVxwPlhpvMMFvf+vAhBt27B2LbtrvqWokS7nByspYwFSUHm1kiIqL0JiEOWJgRUERpz/OsDzTeBFg6GD6XgZw//wLe3gF4+vTLAx+GDKmIyZNr89G0RojNLBERUXoRegsIbAl8uKc9zzkP0PmGyZ6JBT5dVjBz5lmMHHkYCQkqAECGDDbw82uKxo0LSJyOkovNLBERkalTKoCVBYGwx7rnN/IHCrQxbCYDe/8+Bj4+27Fr1311rXJlD2zY0BIeHk4SJqOfxWaWiIjIlF2ZCxwdpHterflAqX6GzSOh//57o/73yJGV8eefNSGX87ICY8dmloiIyBSpEoD5LrpHKDDxBx7okiGDDfz9W6FFC3+sWNEU9evnlToSpRA2s0RERKbm3V3Ar5B2PfevQLOdJjdOrC4hIVFQqQTc3OzVtV9+yY7HjwfxkbQmhg9NICIiMiXPj+luZHsHfzojmw4a2RMnnqJkycVo124LlEqVxjw2sqaHzSwREZEpCL4IzJABm2pq1ksNBIYIwM5NmlwGpFSqMHHiCdSsuQqvXkXg6NEg/PPPGaljUSrjnydERETGShkPnJsInPtL9/wmW4F8zQ2bSSLBwZHo2HErDh9+oq7VqpULPj4lpQtFBsFmloiIyBj5VwdenEh8fsdLgFsZw+WR0OHDj9Ghw1a8efPpIRBmZjKMH18df/xRFebm/BDa1LGZJSIiMiaPdgLbm2jX5XZApuJA5b+AnLUNn0sCSqUKf/55HH/9dQJCfKplyWKP9etbokYNT0mzkeGwmSUiIjIG4c+ApTl1zys3HKg8ETCXGzaThGJjE1C//locP/5UXatXLw/WrGmOzJntJExGhsZmloiIKC27vwXY2Ur3PCtnoN87QJb+Pkq3trZA/vwZcfz4U5ibyzBxYi0MH14ZZmamP1oDaWIzS0RElBZFvgYWZ9U9z8IG6HoPcPQwbKY0Zs6c+nj5MgKjRlVBlSo5pI5DEmEzS0RElJYIFbClPvD0oO75tRcAJfsaNlMa8Px5GO7cCUW9ennUNRsbOXbvbi9hKkoL2MwSERGlBREvgb0dPz30QJeeQYBjItfMmrjdu++jc+ftiI9X4vLl35A/f0apI1Eakv4usiEiIkpLjg/79LCDJdl1N7JNtn566EE6bGQVCiWGDj2ARo024P37GERGxmPYsETOWFO6xTOzREREhiYEcHUecHRQ4svUmAmUHpQub+4CgKCgj2jbNgDnz79U15o1K4gVK3QMS0bpGptZIiIiQxECeLwb2NcZiP2gPV9uB9ScAxTrbvhsacj27XfRtesOfPwYCwCQy83wzz/1MGBAechkHK2ANLGZJSIiMpRVRYF3t7XreZoCTbcB6bxRi4tLwIgRhzBnznl1LXduF/j7t0LZsomM7EDpHptZIiKi1HbXH9jdVrv+y/+ACv8DLKwMnykNatVqM3btuv/VdGEsW9YYTk7WEqaitI7NLBERUWp5cQrwr6p7Xt93gE0Gw+ZJ437/vQJ2774PS0tzzJrlhd69y/KyAvohNrNEREQpTaX89NSuh9t1z/89DjC3NGgkY1C7dm7Mm9cAlSvnQMmS7lLHISPBZpaIiCilqBKA22uA/d2055lbAb1fA9Yuhs+VBj148A5Ll17BtGl1NM6+9utXXsJUZIzYzBIREf0sVQKwrREQtF/HTBnQ8TLgVsrgsdKqDRtu4LffdiEyMh5Zsthj8OCKUkciIyb54HULFiyAp6cnrK2tUaFCBVy4cOG7y8+ePRsFChSAjY0NPDw8MHjwYMTGxhooLRER0TfOTQJmyXU3sgW8gSEqNrL/LyZGgZ49A9G+/VZERsYDAPz8rkOhUEqcjIyZpGdm/f394evri0WLFqFChQqYPXs2vLy8cO/ePWTOnFlr+fXr12PkyJFYsWIFKlWqhPv376NLly6QyWSYOXOmBF8BERGlWy/PABsr655XbgRQ1hew1f5dll49fx6LSpX8cOtWiLrWuXMJLFjQEHK5uYTJyNhJ2szOnDkTPXv2RNeuXQEAixYtwu7du7FixQqMHDlSa/kzZ86gcuXKaN++PQDA09MT7dq1w/nz57WWJSIiShV31gN7OuieV33GpyaWNKxZcwNDh95HXJwKAGBrK8eCBQ3RpUtJaYORSZCsmY2Pj8fly5cxatQodc3MzAx16tTB2bNnda5TqVIlrF27FhcuXED58uXx+PFj7NmzB506dUp0P3FxcYiLi1NPh4eHAwAUCgUUCkUKfTWJ+7wPQ+yLUgePofHjMTR+aeUYWqzMB1nEU626yFQcCd5nAXM5wO8ztaioeAwadACrV/+nrhUunAnr1zdH4cKukh9PSjpDvwf12Y9kzWxoaCiUSiXc3Nw06m5ubrh7967Oddq3b4/Q0FBUqVIFQggkJCSgd+/e+OOPPxLdz5QpUzBhwgSt+oEDB2Bra/tzX4QeDh48aLB9UergMTR+PIbGT6pjaKMIQe2nfSBDgkZdIbPBcY+ZiLLMAuzn99e3Vq16hW3b3qqn69TJgJ49syIo6CKCgqTLRclnqPdgdHR0kpc1qtEMjh07hsmTJ2PhwoWoUKECHj58iEGDBuGvv/7CmDFjdK4zatQo+Pp++cgnPDwcHh4eqFevHhwdHVM9s0KhwMGDB1G3bl3I5fJU3x+lPB5D48djaPykPIay50dgsa2nRk3I7ZDQ+z0gk6G6QdMYlypV4nDjxgq8fh2J337LgokTvfkeNFKGfg9+/iQ9KSRrZjNlygRzc3O8efNGo/7mzRu4u+seKHnMmDHo1KkTevToAQAoVqwYoqKi8Ntvv2H06NEwM9MenMHKygpWVtqPCZTL5QZ9Qxl6f5TyeAyNH4+h8TPoMUzs6V2lBkBWay74naRNCKExZmzGjHJs3eoNmUzg0aPzfA+aAEMdQ332IdnQXJaWlihTpgwOHz6srqlUKhw+fBgVK+oeby46OlqrYTU3/3QHpBAi9cISEVH6cX0RsNBVdyNbcTxQa67BIxmD69eDUanSCjx7FqZRL1bMDQUKZJQoFaUHkl5m4OvrCx8fH5QtWxbly5fH7NmzERUVpR7doHPnzsiWLRumTJkCAGjcuDFmzpyJUqVKqS8zGDNmDBo3bqxuaomIiJLl+mLgUO/E53e5BWQsbLg8RkIIgcWLL+P33/chLk6Jdu224NgxHw63RQYjaTPr7e2NkJAQjB07FsHBwShZsiT27dunvins2bNnGmdi//e//0Emk+F///sfXr58CVdXVzRu3BiTJk2S6ksgIiJj9+ocsCGRJ1BlrQQ02QLY6b78Lb0LC4vFb7/twqZNt9S12NgEvH8fAzc3ewmTUXoi+Q1g/fv3R//+/XXOO3bsmMa0hYUFxo0bh3HjxhkgGRERmbRXZ4EDPYB3t3XP7xMC2GYybCYjcvnyK3h7B+DRow/q2oAB5TF9el1YWUneXlA6wu82IiJKP4QAHgUCO5rpnu9ZH2jkD1il/mg3xkoIgfnzL2Do0IOIj//0GFpnZ2usWNEEzZsXkjgdpUdsZomIKH24Mg84OlD3POc8wK8bAPdyhs1kZD58iEH37oHYtu3LePDly2eDv38reHo6SxeM0jU2s0REZNrCnwFLcyY+v+MVwK2U4fIYsTNnnms0skOGVMTkybVhacmbvUg6kg3NRURElOpurdLdyBbrCXS9CwwRbGT18Ouv+TFoUAVkyGCDwMC2+OefemxkSXI8M0tERKYn5j2wUMfYpo45gR6PARnP5SRFREQc7O0tNR6E8PffdTF0aCVkz87riilt4LuZiIhMh0oJrCmtu5GtNR/oGcRGNonOnHmOIkUWYsWKqxp1S0tzNrKUpvDMLBERmYY3V4C1ZXTP6/cesHYxbB4jpVIJTJ9+GqNHH4FSKTBgwF5UqJAdRYtmljoakU5sZomIyLglxALL8wCRr7TnVZsOlBtq+ExGKiQkCp07b8e+fQ/VtbJls8LFxVrCVETfx2aWiIiMkyoB2FwbeHFCe175kUDVKYbPZMROnHiKdu224NWrCACATAaMHl0V48bVgIUFL82gtIvNLBERGZ+HgcCOprrndbgIuJc1bB4jplSqMGXKKYwbdwwqlQAAZM5sh3XrWqBOndwSpyP6MTazRERkPFRKYFEOICpYe16TbUC+ZgaPZMzevo1Chw5bcejQY3WtVq1cWLu2ObJkcZAwGVHS8XMDIiIyChljbkI+30a7kS0z5NN4sWxk9WZuLsPdu6EAADMzGSZMqIEDBzqykSWjwjOzRESU5snub0aVl//TntH7NWDnbvhAJiJjRlts2NAS7dtvwerVzVGjhqfUkYj0xmaWiIjSruBLwIaKsFAlaNaL9QDqLZUmkxF79SoCFhZmyJzZTl2rUiUHHjwYACsrtgRknHiZARERpU2vLwDryn0ateBrXivYyCbDgQOPULLkInTsuFV9o9dnbGTJmLGZJSKitCUhDtjbGVhfQaP81qYkFD1eAkW7ShTMOCUkqPDHH4fh5bUWISHROHjwMWbPPid1LKIUwz/FiIgobQi+CJwZBzzZqzUrwWs1zj5yRENbVwmCGa8XL8LRrt0WnDr1TF1r2DAfOncuIWEqopTFZpaIiKQjBHCgB3BzReLL+NyEcMoPPNpjuFwmYPfu+/Dx2Y5372IAABYWZpgypTZ8fSvCzEwmcTqilMNmloiIpBEXDsx30j3P3BIo4wtUmQTIzACFwrDZjJhCocQffxzGP/+cVddy5HDCxo0tUbGih4TJiFIHm1kiIjK88Oda18QCADKXBhpvApzzGD6TCYiOVqB27dU4d+6Futa0aQGsWNEUGTLYSJiMKPWwmSUiIsNJiAPmWGvXc9YFmu/6dEaWks3WVo5ChTLh3LkXkMvNMH16XQwcWAEyGS8rINPFZpaIiAzj+iLgUB/teuWJwC+jDZ/HRM2f3xAhIdEYO7YaypXLJnUcolTHZpaIiFJX9FvgXzfd89qfB7KUN2weE/L48Qc8ePAOXl551TVbWzl27mwnYSoiw+I4s0RElHrOT9bdyFaeCAwRbGR/QkDAbZQqtRitW2/Gw4fvpY5DJBmemSUiopT34SGwIp923c4d6PUK4DWcyRYbm4AhQ/Zj4cJL6tqoUYexeXNrCVMRSYfNLBERpRylAlhVDPhwT3tetelAuaGGz2RCHjx4B2/vAFy9GqyutW1bFIsXN5IwFZG02MwSEdHPEwLY1wW4vVp7nnUGoMcTwMrR4LFMycaNN9Gz505ERsYDAKytLTB3bn306FGaoxVQusZmloiIfo5QAVvqA08Pas9rfQTIUdPwmUxITIwCv/++D0uWXFHXChTIiE2bWqN48URurCNKR9jMEhFR8l37FzjcV7tepCtQbwlgxl8zP6tJk404dOixerpTp+JYuPBX2NtzTF4igM0sERElx4uTgH817XqRLoDXCt7glYKGDq2IQ4cew8bGAgsX/oouXUpKHYkoTWEzS0RESadKAGbJdc8rNxyoNs2wedIBL6+8mD+/AWrWzIXChV2ljkOU5rCZJSKipHl/H1hZQPe8QTGAhY7H1JJebt16i5Urr2H69LoaN3X168fxeIkSw2aWiIi+L/YjsLcT8HiX9rw+bwDbzAaPZGqEEFi58hr699+DmJgE5MjhhIEDK0gdi8go8AlgRESUuHubgQUu2o1s2WGfnuDFRvanRUbGo3Pn7ejePRAxMQkAgDVr/oNSqZI4GZFx4JlZIiLSpogC5trrntdgDVC4o2HzmKjr14PRpk0A7t9/p6716lUGs2Z5wdyc55uIkoLNLBERaUpsuK22p4BslQ2fxwQJIbBkyWUMGrQPcXFKAICDgyWWLGmMtm2LSpyOyLiwmSUiok+e7AO2NtA9jzd4pZjw8Dj89ttO+PvfUtdKl84Cf/9WyJs3g4TJiIwTm1kiovQu5h2wMJPueb+MASr/adg8Jm7s2KMajWz//uXwzz/1YGXFX8lEycF3DhFRevZwB7Cjme55v8cB5nzKVEqbMKEGdu68j3fvorF8eRO0bFlY6khERo3NLBFReiRUwIbKwOtz2vO63Qdc8hk+k4kSQmiMGevkZI1t27zh4GCJXLlcJExGZBp4qyQRUXpzZx0w01yzkZXbA60Pfxpui41sirlw4SXKl1+GFy/CNerFi7uxkSVKITwzS0SUXggBzEzkHEb/D4AZfyWkFCEEZs8+hxEjDkGhUKFduy04etQHFhY8h0SU0viuIiJKD57s093IFvAGfFVsZFPQ+/cxaNbMH76+B6BQfHrwgVKpwsePsRInIzJN/OlFRGTKQm4Aq4vrnjdYwSY2hZ09+xze3gF4/vzLZQXDh1fCxIm1IJebS5iMyHTxpxgRkana1xW45addz9UAaL4b+OqmJPo5KpXAP/+cwR9/HIZSKQAAGTPaYPXq5mjYkNcgE6UmNrNERKbm9hpgb2fd83q9AuyzGDaPiQsJiYKPz3bs3ftQXatSJQc2bGiJ7NkdJUxGlD6wmSUiMhWKGGCure55TXcAeZsYNk86cebMc3UjK5MBf/xRFePH1+DNXkQGwncaEZEpODFSdyPrnAfoH8ZGNhU1bVoQ/fuXQ+bMdti/vyMmTqzFRpbIgHhmlojI2F1dAFycpl0fEA5YOhg+j4kLC4uFk5O1Ru2ff+ph9OhqcHe3lygVUfrFPx2JiIzZMV/gSH/NWs25nx5+wEY2xR09+gQFCy6An981jbqVlQUbWSKJ8MwsEZExSmzILQ63lSqUShUmTjyBP/88AZVKoF+/PShfPhsKF3aVOhpRusefeERExiQhDphjrXvegAg2sqng9esIdOiwFUePBqlrlSt7IFOmRG62IyKD4mUGRETG4nB/3Y2sa/FPZ2Qt+TF3Sjt48BFKllysbmTNzGSYNKkW9u3riMyZ7aQNR0QAeGaWiCjte3YE2Fxb97zujwDn3IbNkw4kJKgwfvwxTJ58EuLTMxCQLZsDNmxoiapVc0objog0sJklIkqrVEpgWS4g4rn2vCpTgAojDZ8pHXj9OgLe3gE4efKZutagQV6sXt2clxYQpUFsZomI0qLHe4Btv2rXLWyAgZGAjFeJpRYLCzM8evQBAGBuLsOUKbUxZEglmJnx8b9EaRF/GhIRpTUXp+tuZH97AQyKZiObylxd7bBhQ0vkyuWMkye7YtiwymxkidIwnpklIkpLzk0CTv9Ps5ajNtDqAJvYVPLsWRhsbCzg6vrlhq5q1XLi3r3+kMvNJUxGREnxU81sbGwsrK0TGSKGiIj0s6oYEHrzy7SlI9AzCLB2kSySqQsMvIcuXbajQoXs2L27vcYZWDayRMZB7z/zVSoV/vrrL2TLlg329vZ4/PgxAGDMmDFYvnx5igckIkoXluXWbGQBoOtdNrKpJD5eicGD96Fp04348CEW+/Y9xMKFF6WORUTJoHczO3HiRPj5+eHvv/+GpaWlul60aFEsW7YsRcMREZk8lRJYmgsIe6JZ7/MWsM8iTSYT9+TJB1SpsgKzZ59X11q2LISOHXU8UY2I0jy9m9nVq1djyZIl6NChA8zNv3wEU6JECdy9ezdFwxERmbTQW8AsCyA8SLM+IByw5WNSU8PWrXdQqtRiXLz4CgBgaWmO+fMbYPPm1nB25mVzRMZI72tmX758ibx582rVVSoVFApFioQiIjJ59zYDu9po1rJVBdqekCaPiYuNTcCwYQcwf/6XSwny5HHBpk2tUbo0z4ATGTO9m9nChQvj5MmTyJlT8wkoAQEBKFWqVIoFIyIySUIA68oDby5p1vM2B5pulSaTiYuIiEP16n64ejVYXfP2LoIlSxrD0dFKwmRElBL0bmbHjh0LHx8fvHz5EiqVClu3bsW9e/ewevVq7Nq1KzUyEhGZhqADwBYv7XqTLUC+FobPk044OFihWDE3XL0aDCsrc8yd2wA9e5aGTMaxY4lMgd7NbNOmTbFz5078+eefsLOzw9ixY1G6dGns3LkTdevWTY2MRETG78ND3Y2szw0gU1HD50lnFi5siLCwWPz5Z00UL+4mdRwiSkHJGme2atWqOHjwYEpnISIyTZdnA8cGa9YyFQN8/pMkjqm7dy8UT5+GoV69POqanZ0ltm9vK2EqIkoteo9mkDt3brx7906r/vHjR+TOnTtFQhERmYw7G7Qb2V/GsJFNJWvX/ocyZZagTZvNePz4g9RxiMgA9D4zGxQUBKVSqVWPi4vDy5cvUyQUEZHRU8YDs3XcXNRgDVC4o+HzmLjoaAX699+DlSuvqWvjxh3DmjXNpQtFRAaR5GY2MDBQ/e/9+/fDyclJPa1UKnH48GF4enqmaDgiIqMU9uTTE72+1ekakLmEweOYulu33qJNmwDcvh2irnXtWhLz5jWQMBURGUqSm9lmzZoBAGQyGXx8fDTmyeVyeHp6YsaMGSkajojI6ITcAFbreJJUzyDAMad2nZJNCAE/v2vo128PYmISAAB2dnL8+++v6NSJfzQQpRdJbmZVKhUAIFeuXLh48SIyZcqUaqGIiIzS8+PAphqatRJ9gDoLJYljyiIj49G3726sWfPl2uNixTJj06bWKFiQv5+I0hO9r5l98uTJjxciIkpPFDHAmlLAh3ua9doLgZJ9pMlkwoQQaNhwHU6efKau9epVBrNmecHGRi5hMiKSgt6jGQBAVFQU9uzZg0WLFmHu3Lka/+lrwYIF8PT0hLW1NSpUqIALFy58d/mPHz+iX79+yJIlC6ysrJA/f37s2bMnOV8GEdHPuzANmGur3ci22MNGNpXIZDKMHFkFAODgYIkNG1pi0aJGbGSJ0im9z8xevXoVDRs2RHR0NKKiopAhQwaEhobC1tYWmTNnxsCBA5O8LX9/f/j6+mLRokWoUKECZs+eDS8vL9y7dw+ZM2fWWj4+Ph5169ZF5syZERAQgGzZsuHp06dwdnbW98sgIvo5Hx8By/PqntdiD5CLNx+lpoYN82H+/Abw8sqLvHkzSB2HiCSk95nZwYMHo3Hjxvjw4QNsbGxw7tw5PH36FGXKlME///yj17ZmzpyJnj17omvXrihcuDAWLVoEW1tbrFixQufyK1aswPv377F9+3ZUrlwZnp6eqF69OkqU4IX+RGRAr87pbmTLDQeGCDayKezq1dcYMeIwhBAa9X79yrORJSL9z8xeu3YNixcvhpmZGczNzREXF4fcuXPj77//ho+PD1q0SNrzxePj43H58mWMGjVKXTMzM0OdOnVw9uxZnesEBgaiYsWK6NevH3bs2AFXV1e0b98eI0aMgLm5uc514uLiEBcXp54ODw8HACgUCigUiqR+2cn2eR+G2BelDh5D45eSx1D28iQsttTWqIkMhZDQ9gJgYQXw+yTFCCGwaNFlDBt2GPHxSkRFZUe9enx9jRF/jho/Qx9DffajdzMrl8thZvbphG7mzJnx7NkzFCpUCE5OTnj+/HmStxMaGgqlUgk3N81nZLu5ueHu3bs613n8+DGOHDmCDh06YM+ePXj48CH69u0LhUKBcePG6VxnypQpmDBhglb9wIEDsLW1TXLen8XH/xo/HkPj9zPH0EbxFuWCp8Ml7oFG/b5LS9zJ0Ak4cPhn49FXIiMTsGDBc5w9G6aunTz5Afv3H4CZmUzCZPQz+HPU+BnqGEZHRyd5Wb2b2VKlSuHixYvIly8fqlevjrFjxyI0NBRr1qxB0aJF9d2cXlQqFTJnzowlS5bA3NwcZcqUwcuXLzF9+vREm9lRo0bB19dXPR0eHg4PDw/Uq1cPjo6OqZoX+PSXxcGDB1G3bl3I5bw5wRjxGBq/nzqGQsB8Z3OYPdW80VRY2CCh/SXkcs6HXCmYlYBLl15h8ODtePLkSyPbr18Z1KihgJdXPb4PjRB/jho/Qx/Dz5+kJ4XezezkyZMREREBAJg0aRI6d+6MPn36IF++fFi+fHmSt5MpUyaYm5vjzZs3GvU3b97A3d1d5zpZsmSBXC7XuKSgUKFCCA4ORnx8PCwtLbXWsbKygpWV9iMl5XK5Qd9Qht4fpTweQ+On9zG8swHY017nLFmna5BnyJ9CyQj4dFnBnDnnMXz4QSgUn8Y2d3a2hp9fUzRsmAd79uzh+9DI8fgZP0MdQ332oXczW7ZsWfW/M2fOjH379um7CQCApaUlypQpg8OHD6ufLqZSqXD48GH0799f5zqVK1fG+vXroVKp1Jc63L9/H1myZNHZyBIRJYtQATN1X4ePhuuAgu0AGT/qTknv38ega9cdCAz8MsTZL79kx8aNLZEzpzOvtSSiRCVrnFldrly5gkaNGum1jq+vL5YuXYpVq1bhzp076NOnD6KiotC1a1cAQOfOnTVuEOvTpw/ev3+PQYMG4f79+9i9ezcmT56Mfv36pdSXQUTp3ZO9uhvZyn8BviqgUHs2sqlg9OjDGo3s8OGVcOJEF+TM6SxdKCIyCnqdmd2/fz8OHjwIS0tL9OjRA7lz58bdu3cxcuRI7Ny5E15eXnrt3NvbGyEhIRg7diyCg4NRsmRJ7Nu3T31T2LNnz9RnYAHAw8MD+/fvx+DBg1G8eHFky5YNgwYNwogRI/TaLxGRlui3wL9uuucNjALkhrthND2aPLk29u17hIiIOKxe3RwNG+aTOhIRGYkkN7PLly9Hz549kSFDBnz48AHLli3DzJkzMWDAAHh7e+PmzZsoVKiQ3gH69++f6GUFx44d06pVrFgR586d03s/RESJ+vgYWKXjBtYiPoDXSp6JTQVCCMi+el1dXGywfbs3Mma0RfbsqX9zLhGZjiRfZjBnzhxMmzYNoaGh2LRpE0JDQ7Fw4ULcuHEDixYtSlYjS0QkuftbgOV5gISYL7WMRYBeL4H6fmxkU8HJk09RpswSvHoVoVEvUcKdjSwR6S3JzeyjR4/QunVrAECLFi1gYWGB6dOnI3v27KkWjogoVb0+D+xspVlrcwzochOwzypJJFOmUglMnnwSNWuuwtWrwWjffguUSpXUsYjIyCX5MoOYmBj1QwZkMhmsrKyQJUuWVAtGRJSqPjwA1v+iWWsWCHhUlyaPiXv7NgqdOm3DgQOP1DWZTIbw8Di4uNhImIyIjJ1eN4AtW7YM9vb2AICEhAT4+fkhU6ZMGssMHDgw5dIREaWGiJfAim/GiG17CshWWZo8Ju7o0Sdo334rgoMjAXy6cmPs2OoYM6YazM1TbFAdIkqnktzM5siRA0uXLlVPu7u7Y82aNRrLyGQyNrNElLadnwycGq1Z81rJRjYVKJUqTJx4An/+eQIqlQAAuLvbY926FqhVi89NI6KUkeRmNigoKBVjEBGlPvlcHQ9XydsMKNrF0FFM3uvXEejYcRuOHHmirtWpkxtr1zaHm5u9hMmIyNTo/QQwIiKjE/cRTR82067X9/s0/BaluDNnnqsbWTMzGf78swZGjaoKMzOODkFEKYvNLBGZtqg3kC921673eQPYZjZ8nnSiZcvC6N27DAID72PDhpaoVi2n1JGIyETxynsiMl1hT4BF3zSyBbyBIYKNbAr78CFGqzZrVn1cu9aLjSwRpSo2s0RkmhJigWW5NUrK8n8AjTZKFMh07d37APnzz8fatf9p1K2tLeDqaidRKiJKL9jMEpHpiX4LzHPQKF10GwrVL+OlyWOiFAolRow4iIYN1yM0NBq9e+/C3buhUscionQmWc3so0eP8L///Q/t2rXD27dvAQB79+7FrVu3UjQcEZHenh0B/nUDVAnqkqpwV7xyqCJhKNPz7FkYatRYhb//PqOu1aqVC66uthKmIqL0SO9m9vjx4yhWrBjOnz+PrVu3IjLy0yDY169fx7hx41I8IBFRkl2ZB2yurVmrMgXKOoulyWOiAgPvoWTJRThz5jkAwMLCDDNn1sOOHW2RMSObWSIyLL2b2ZEjR2LixIk4ePAgLC2/jNlYq1YtnDt3LkXDEREl2eH+wNFvHtritQKoMFKaPCYoPl4JX9/9aNp0Iz58iAUAeHo64/Tpbhg8uCJkMg67RUSGp/fQXDdu3MD69eu16pkzZ0ZoKK+VIiIJ3NkAXFugWev5DHD0kCaPCXr2LAytW2/GhQsv1bUWLQph+fImcHa2ljAZEaV3ep+ZdXZ2xuvXr7XqV69eRbZs2VIkFBFRkl2ZC+xpr1nr94GNbAqzsjLHs2dhAABLS3PMm9cAAQGt2cgSkeT0bmbbtm2LESNGIDg4GDKZDCqVCqdPn8bQoUPRuXPn1MhIRKQt6g0wQwYcHaRZ7/kUsHaWJJIpc3Ozx/r1LZA/f0acOdMN/fuX52UFRJQm6N3MTp48GQULFoSHhwciIyNRuHBhVKtWDZUqVcL//ve/1MhIRKRJEaX9MAQAaHMMcMxh8Dim6NGj9wgNjdao1ayZC7du9UWZMlklSkVEpE3va2YtLS2xdOlSjBkzBjdv3kRkZCRKlSqFfPnypUY+IiJNwZeAdeW06wPCAUsH7TrpbdOmW+jRIxDVquVEYGA7mJl9OQNrYcHhyYkobdG7mT116hSqVKmCHDlyIEcOngEhIgN6d1u7kS3cCWiwWpo8JiYmRgFf3/1YtOgyAGD37gdYuvQyevUqK3EyIqLE6f0ndq1atZArVy788ccfuH37dmpkIiLS9jAQ8CuiWas8kY1sCrl3LxS//LJc3cgCQIcOxdC+fTEJUxER/ZjezeyrV68wZMgQHD9+HEWLFkXJkiUxffp0vHjxIjXyERF9utFrR1PNWssDwC+jpcljYtat+w9lyizBf/+9AQDY2Fhg+fImWLOmORwcrCROR0T0fXo3s5kyZUL//v1x+vRpPHr0CK1bt8aqVavg6emJWrVqpUZGIkrPVhXXrjULBDzrGj6LiYmOVqBHj0B07LgNUVEKAEChQplw4UJPdOtWiqMVEJFR0Pua2a/lypULI0eORIkSJTBmzBgcP348pXIRUXonBDDXHkjQvKMevYMBOzdpMpmQjx9jUaXKCty6FaKudelSEvPnN4CdneV31iQiSluSfVvq6dOn0bdvX2TJkgXt27dH0aJFsXv37pTMRkTp2b9u2o3sEMFGNoU4OVmhRIlPw5vZ2sqxalUzrFzZlI0sERkdvc/Mjho1Chs3bsSrV69Qt25dzJkzB02bNoWtrW1q5COi9EaogJnm2vVBsYbPYsJkMhkWLfoVsbEJmDSpFgoWzCR1JCKiZNG7mT1x4gSGDRuGNm3aIFMm/vAjohSUEAcscNGsWToC/T8CvH7zp9y48QavX0eiXr086pqDgxW2bGkjYSoiop+ndzN7+vTp1MhBRATMsdactrAGBoRJk8VECCGwbNkVDBy4D9bWFrh6tRc8PZ2ljkVElGKS1MwGBgaiQYMGkMvlCAwM/O6yTZo0SZFgRJTO7O6gOZ2rAdBijzRZTERERBx69dqFDRtuAgBiYxPw11/HsXx50x+sSURkPJLUzDZr1gzBwcHInDkzmjVrluhyMpkMSqUypbIRUXqgiAHmfnPNvaUDG9mfdPXqa7RpE4CHD9+ra337lsWMGV4SpiIiSnlJamZVKpXOfxMR/ZTHe4Btv2rX+73XrlGSCCHw77+X4Ou7H3Fxn04uODpaYdmyxmjdusgP1iYiMj56D821evVqxMXFadXj4+OxejUfK0lESaBKANaW093IDggHzH5qCOx0KywsFm3aBKBfvz3qRrZs2ay4erUXG1kiMll6N7Ndu3ZFWJj2DRkRERHo2rVrioQiIhO3oTLw5pJmrYD3p3FkLR2kyWTkhBCoW3cNAgJuq2uDBlXAqVNdkTu3y3fWJCIybno3s0IInY84fPHiBZycnFIkFBGZsJOjgOALmrVuD4BGG6XJYyJkMhnGjKkGAHB2tsa2bd6YPbs+rKx4lpuITFuSf8qVKvXpOd0ymQy1a9eGhcWXVZVKJZ48eYL69eunSkgiMhEXpgEXpn6ZdsoN9HgkXR4T07hxASxY0BANG+bj8FtElG4kuZn9PIrBtWvX4OXlBXt7e/U8S0tLeHp6omXLlikekIhMgBDATB0fBHV/aPgsJuLcuRfYtOkWZsyop/FpWd++5SRMRURkeEluZseNGwcA8PT0hLe3N6ytrX+wBhERgPhIYJ6O62A7XORTvZJBpRKYMeMM/vjjCBISVChQICN69SordSwiIsnofc2sj48PG1kiSpoXp3Q3sn1DAXc2YPoKDY1GkyYbMHz4ISQkfBomMSDgDoQQEicjIpJOks7MZsiQAffv30emTJng4uKi8wawz96/5/iQRATgyV5ga0PNmksBoNtdafIYuVOnnqFduy148SJcXRs1qgr+/LPmd38mExGZuiQ1s7NmzYKDg4P63/zBSUTftb8HcHO5Zq3ecqBYN2nyGDGVSmDatFMYM+YolMpPZ2BdXW2xZk1zeHnllTgdEZH0ktTM+vj4qP/dpUuX1MpCRKZgdSkg5JpmrcUeIFcDSeIYs7dvo9Cp0zYcOPBlxIfq1XNi/fqWyJqV4/ESEQHJuGb2ypUruHHjhnp6x44daNasGf744w/Ex8enaDgiMjI7W2s3stX+ZiObTH/8cVjdyMpkwNix1XDoUGc2skREX9G7me3Vqxfu378PAHj8+DG8vb1ha2uLzZs3Y/jw4SkekIiMgBDAqmLA/QDN+qBYoNwwaTKZgL//roscOZzg5maHgwc7YcKEmrCw0PvHNhGRSdP70TD3799HyZIlAQCbN29G9erVsX79epw+fRpt27bF7NmzUzgiEaV5SzyAyJeatcEKwIxPn9KHSiVgZvblnoQMGWwQGNgWbm72cHe3/86aRETpV7IeZ6tSfRoS5tChQ2jY8NPdyh4eHggNDU3ZdESUtimigRkyNrIp4NChxyhVajGCgyM16iVKuLORJSL6Dr2b2bJly2LixIlYs2YNjh8/jl9//RUA8OTJE7i5uaV4QCJKg5QKYGkuYK6d9rzBCWxk9ZCQoMKYMUdQr94a/PffG3TosBVKpUrqWERERkPv3zizZ89Ghw4dsH37dowePRp5834aGiYgIACVKlVK8YBElAYty6V9NhYAejwGzMwNn8dIvXwZjvbtt+LEiafqmqWlOaKiFHB0tJIwGRGR8dC7mS1evLjGaAafTZ8+Hebm/CVGZNLe3QX8CmnXy/gCVacC5nLDZzJS+/Y9RKdO2xAaGg0AMDeXYdKkWhg2rLLGdbNERPR9yf4s8PLly7hz5w4AoHDhwihdunSKhSKiNOjo78CVOZq1DAWBrnckiWOsFAolxow5imnTTqtr2bM7YuPGlqhcOYeEyYiIjJPezezbt2/h7e2N48ePw9nZGQDw8eNH1KxZExs3boSrq2tKZyQiKQkBzNRxeb19NqDLbcPnMWLPn4ehbdstOHPmubrWqFF++Pk1RcaMthImIyIyXnrfADZgwABERkbi1q1beP/+Pd6/f4+bN28iPDwcAwcOTI2MRCQVodLdyDYLBHq9+DSSPyXZmTPP1Y2shYUZZsyoh8DAtmxkiYh+gt5nZvft24dDhw6hUKEv180VLlwYCxYsQL169VI0HBFJ6MMDYEV+7XqfEMA2k+HzmABv76I4fPgJDhx4BH//VqhQIbvUkYiIjJ7ezaxKpYJcrn2Th1wuV48/S0RG7sVJwL+aZs0xJ9DjCc/G6uHdu2its65z5tRHbGwCXFxsJEpFRGRa9L7MoFatWhg0aBBevXqlrr18+RKDBw9G7dq1UzQcERmYEMCZ8dqNbOlBbGT1tHXrHeTJMxcbNmiO/mJjI2cjS0SUgvRuZufPn4/w8HB4enoiT548yJMnD3LlyoXw8HDMmzcvNTISkSEkxAEzzYGzEzTrFccBNWezkU2iuLgEDBiwBy1bbkJYWBx++20XHjx4J3UsIiKTpfdlBh4eHrhy5QoOHz6sHpqrUKFCqFOnToqHIyIDEQKYY61d73IHyFjQ8HmM1KNH7+HtHYDLl1+raw0b5kPmzDqelEZERClCr2bW398fgYGBiI+PR+3atTFgwIDUykVEhiIEMN9Zs2YmBwZEABZ8ClVSbdp0Cz16BCIiIh4AYGVljtmz66NXrzKQ8aw2EVGqSXIz+++//6Jfv37Ily8fbGxssHXrVjx69AjTp09PzXxElJqEALbUB+LDv9SsnIH+HySLZGxiYxMwePA+LFp0WV3Lly8DNm1qjZIl3SVMRkSUPiT5mtn58+dj3LhxuHfvHq5du4ZVq1Zh4cKFqZmNiFJTdMinMWSfHvhSM5OzkdXD48cf8MsvyzQa2fbti+Hy5d/YyBIRGUiSm9nHjx/Dx8dHPd2+fXskJCTg9evX31mLiNKkGyuAfzNr1wfHGz6LEbO1leP160gAgLW1BZYta4y1a5vDwYGXZxARGUqSm9m4uDjY2X25icHMzAyWlpaIiYlJlWBElEr2dAQOdNeuD04wfBYj5+5uj3XrWqBIEVdcvNgT3buX5vWxREQGptcNYGPGjIGt7ZcBwOPj4zFp0iQ4OTmpazNnzky5dESUsq4vBu6s06y12APkaiBNHiNz504I3NzskSHDl3Fi69TJjWvXesPCQu+RDomIKAUkuZmtVq0a7t27p1GrVKkSHj9+rJ7mGQmiNOyuP3Cot2bttxeAQzZp8hgZP79r6NdvD+rUyY3t2701ft6xkSUikk6Sm9ljx46lYgwiSlVHBwNXZmvWut5lI5sEkZHx6NdvD1avvg4ACAy8Bz+/a+jatZTEyYiICEjGQxOIyMjs6wrc8tOsNQsEMhSQJI4xuXHjDdq0CcDdu6HqWo8epeDtXVTCVERE9DU2s0SmSpUArCkNhN7QrA+MBOR8ItX3CCGwfPlVDBiwF7Gxn26Ms7e3xOLFjdC+fTGJ0xER0dfYzBKZoshXwGIdlxD43GAj+wMREXHo3Xs31q//8kdAiRJu2LSpNfLnzyhhMiIi0oXNLJGpOT8ZODVauz4oBrCwNnweI/LuXTQqVlyOBw/eq2t9+5bFjBlesLbmj0siorSIP52JTMm6CkDwBc1apqKfzsjSD2XIYIPSpbPgwYP3cHS0wrJljdG6dRGpYxER0Xckq5k9efIkFi9ejEePHiEgIADZsmXDmjVrkCtXLlSpUiWlMxLRj4Q/BZZ6atfbHAM8qhs6jdGSyWRYsqQxhAAmT66FPHkySB2JiIh+QO/BEbds2QIvLy/Y2Njg6tWriIuLAwCEhYVh8uTJKR6QiH5AEaO7kR0Uy0b2By5deoUDBx5p1BwdreDv34qNLBGRkdC7mZ04cSIWLVqEpUuXQi6Xq+uVK1fGlStXUjQcEf1AXDgw11a73u8DYGFl+DxGQgiBOXPOoVKl5WjbNgDPnoVJHYmIiJJJ72b23r17qFatmlbdyckJHz9+TIlMRJQUV+YB8500a9X/AYYIwNpZkkjG4P37GDRv7o/ff98PhUKFDx9iMW3aKaljERFRMundzLq7u+Phw4da9VOnTiF37tzJCrFgwQJ4enrC2toaFSpUwIULF368EoCNGzdCJpOhWbNmydovkVESAliRHzg6ULNebgRQdog0mYzE+fMvUarUYuzY8eXR3EOGVMSsWfUlTEVERD9D72a2Z8+eGDRoEM6fPw+ZTIZXr15h3bp1GDp0KPr06aN3AH9/f/j6+mLcuHG4cuUKSpQoAS8vL7x9+/a76wUFBWHo0KGoWrWq3vskMlrKeGCmGfDhgWa9/iqg2lRpMhkBlUpg+/a3qFlzjfqSggwZbLBzZzv88089WFqaS5yQiIiSS+/RDEaOHAmVSoXatWsjOjoa1apVg5WVFYYOHYoBAwboHWDmzJno2bMnunbtCgBYtGgRdu/ejRUrVmDkyJE611EqlejQoQMmTJiAkydP8vIGSh+U8cBsHdfB9n4N2LkbPo+RCA2Nho/PNuzZ80pdq1zZAxs2tISHh9N31iQiImOgdzMrk8kwevRoDBs2DA8fPkRkZCQKFy4Me3t7vXceHx+Py5cvY9SoUeqamZkZ6tSpg7Nnzya63p9//onMmTOje/fuOHny5Hf3ERcXpx5xAQDCw8MBAAqFAgqFQu/M+vq8D0Psi1JHWjmGFn6FIPtqWkCGhAExgMwM4PeXTiqVQM2afrh5M0RdGz68EsaNqwq53FzyY0pJl1beh5Q8PH7Gz9DHUJ/9JPuhCZaWlihcuHByVwcAhIaGQqlUws3NTaPu5uaGu3fv6lzn1KlTWL58Oa5du5akfUyZMgUTJkzQqh84cAC2tjruAk8lBw8eNNi+KHVIdgyFQK1n/eGgeKlRDsy7Ddi7T5pMRqRRIzvcvBkCR0dzDB6cE6VKRePgwf1Sx6Jk4s9S48bjZ/wMdQyjo6OTvKzezWzNmjUhk8kSnX/kyBF9N5lkERER6NSpE5YuXYpMmTIlaZ1Ro0bB19dXPR0eHg4PDw/Uq1cPjo6OqRVVTaFQ4ODBg6hbt67GUGZkPKQ+hvK5lhrTqtxNoGwUgIYGT2KcGjYE3N0vwNHxNby9G/J9aKSkfh/Sz+HxM36GPoafP0lPCr2b2ZIlS2pMKxQKXLt2DTdv3oSPj49e28qUKRPMzc3x5s0bjfqbN2/g7q59DeCjR48QFBSExo0bq2sqlQoAYGFhgXv37iFPnjwa61hZWcHKSvs6Q7lcbtA3lKH3RylPkmN4bpJWyazZNpjJ9L53M104fjwIO3bcw4wZ9TT+6O7btzz27NnD96EJ4DE0bjx+xs9Qx1CffejdzM6aNUtnffz48YiMjNRrW5aWlihTpgwOHz6sHl5LpVLh8OHD6N+/v9byBQsWxI0bms+Y/9///oeIiAjMmTMHHh4eeu2fKE17dRY4/T/Nmq8K+M4nI+mVUqnCpEknMWHCcahUAkWKuKJ799JSxyIiIgNI9jWz3+rYsSPKly+Pf/75R6/1fH194ePjg7Jly6J8+fKYPXs2oqKi1KMbdO7cGdmyZcOUKVNgbW2NokWLaqzv7OwMAFp1IqP2/BiwqaZmbWAkG1kdgoMj0aHDVhw58kRd2779Hrp1K/XdS6KIiMg0pFgze/bsWVhbW+u9nre3N0JCQjB27FgEBwejZMmS2Ldvn/qmsGfPnsHMjB+pUjqREAvMsdGutzsLyO0MnyeNO3ToMTp23Io3b6IAAGZmMowfXx1//FGVjSwRUTqhdzPbokULjWkhBF6/fo1Lly5hzJgxyQrRv39/nZcVAMCxY8e+u66fn1+y9kmU5jw7CmyupV3veAlwK2P4PGlYQoIKEyYcw6RJJyHEp1qWLPbYsKElqlf3lDQbEREZlt7NrJOT5iDjZmZmKFCgAP7880/Uq1cvxYIRpSuHBwDX5mvXewcDdm7a9XTs5ctwtG+/FSdOPFXXvLzyYPXq5sicmWeviYjSG72aWaVSia5du6JYsWJwcXFJrUxE6cvNldqNbN0lQPGe0uRJ40aNOqxuZM3NZZg4sRaGD68MMzNeVkBElB7p1cyam5ujXr16uHPnDptZopSwvRnwaIdmrfN/gGsxSeIYg5kzvXDkyBPIZDJs2NASVarkkDoSERFJSO/LDIoWLYrHjx8jV65cqZGHKH1QKYFZOt5+A6MBuY4bwNIxlUponHXNlMkWu3e3R/bsjsiY0XBP8SMiorRJ72ECJk6ciKFDh2LXrl14/fo1wsPDNf4joiRYU0q79nscG9lv7Np1HyVKLMKbN5pjWJco4c5GloiIAOjRzP7555+IiopCw4YNcf36dTRp0gTZs2eHi4sLXFxc4OzszEsPiH5EiE+XFoRqPvwDg2IBc0udq6RH8fFKDBmyH40bb8DNm2/RqdM2qFRC6lhERJQGJfkygwkTJqB37944evRoauYhMl1CADN1/P04hE3a14KCPsLbOwAXLrxU1+zsLBETo4CdHRt+IiLSlORmVvz/YI7Vq1dPtTBEJkupAGbraMS6PzR8ljRs27Y76NYtEB8/xgIA5HIz/PNPPQwYUJ4PQSAiIp30ugGMv0yIkiH4ErCunGbNTA78HgvI+HQ7AIiLS8CwYQcxb94FdS13bhf4+7dC2bJZJUxGRERpnV7NbP78+X/Y0L5///6nAhGZlEc7ge1NNGsZCwNdbkmTJw169Og9vL0DcPnya3WtdevCWLq0MZyc9H9ENhERpS96NbMTJkzQegIYESXi0kzg+BDNWpGuQP0V0uRJo86de6FuZK2szDFrlhd69y7LT4KIiChJ9Gpm27Zti8yZM6dWFiLTcXwYcOkfzVrdpUDxHtLkScM6dCiOw4ef4NSpZ9i0qTVKlnSXOhIRERmRJDezPEtClET7ugK3/DRrzXYCeRpJEietefs2Cpkz22nU5s9vCKVSBQcHK4lSERGRsUry3SefRzMgou845qvdyPYPYyP7/9avv4E8eeZi0ybNa4ZtbeVsZImIKFmS3MyqVCpeYkD0PSdHAZdnadZ+jwesHKXJk4ZERyvQs2cgOnTYisjIePToEYhHj3izKBER/Ty9rpklokTsaA483K5Z6/8RMJdLkSZNuXMnBG3aBODmzbfqWosWheDubi9hKiIiMhVsZol+1rbGwONdmrUeTwArjvyxatU19O27B9HRCgCfLidYuLAhfHxKShuMiIhMBptZop8xQ8eNkf0/pvtGNioqHn377sHq1dfVtSJFXLFpU2sULuwqYTIiIjI1bGaJkiMhDpijY0D/316k+0b23r1QNGvmj7t3Q9W1Hj1KYc6cBrC15WUXRESUstjMEiWHrka23wfA2tngUdIaBwcrvHsXDQCwt7fE4sWN0L59MYlTERGRqeKD4Yn0IQSwu7123VfFRvb/Zc3qgDVrmqNUKXdcvvwbG1kiIkpVPDNLlFQqJTDTUrs+JH2PwXz9ejBy5HCCi4uNuubllRd16uSGuTn/XiYiotTF3zRESSATSlhsKKs9o9crw4dJI4QQ+Pffi6hQYRm6dQvUerAKG1kiIjIE/rYhSoImj1pC9k7zqVX4PQ6wzyJNIImFhcXC2zsAffvuQVycEtu338W6dTekjkVEROkQLzMg+h5FDORzbTVrVacB5YdLkycNuHTpFby9A/D48Qd1bcCA8mjdurCEqYiIKL1iM0uUmBenAP+qmjXnPOm2kRVCYN68Cxg69AAUChUAwNnZGitWNEHz5oUkTkdEROkVm1kiXS7+A5wYplFSFewAs1/XShRIWh8+xKB790Bs23ZXXStfPhv8/VvB09NZumBERJTusZkl+pZQaTWy951bIle9lenyIvM3byJRocIyPH0apq4NGVIRkyfXhqWluYTJiIiI2MwSaUqIBebYaJYabsSd+9bIJVEkqWXObIdy5bLh6dMwZMhgAz+/pmjcuIDUsYiIiACwmSX6IuQGsLq4Zq3SBIi8LYD7e6TJlAbIZDIsW9YYcrkZpk6tgxw50vfjeomIKG1hM0sEAM+PAZtqatYcPYFfxgAJCRIEks7p088QHa1A3bp51DUnJ2usX99SwlRERES6pcdLAIk03duk3ciWHQb0fALIZNJkkoBKJTB16ilUr+6Hdu224MWLcKkjERER/RCbWUrf1pQGdnlr1mrMAqr/LU0eiYSEROHXX9dj1KjDUCoF3r2LwcyZZ6WORURE9EO8zIDSp/gIYF0F4P0dzXrn64Brcd3rmKjjx4PQvv1WvHoVAeDTyejRo6ti3Lga0gYjIiJKAjazlP68uQqsLa1d7/cBsHY2eBypKJUqTJ58EuPHH4dKJQAAbm52WLu2BerUyS1xOiIioqRhM0vpy7u72o1s5lJAx8vp6vrY4OBIdOy4FYcPP1HXatXKhXXrWsDd3V7CZERERPphM0vpx4cHgN83j131WgEU7SpNHokolSrUrLkKd++GAgDMzGQYN646Ro+uCnNzXkZPRETGhb+5KP1YkV9z+pcx6a6RBQBzczNMnPhp9IYsWexx+HBnjB1bnY0sEREZJZ6ZJdMnBDDzm0atwWqgcCdp8qQBLVsWxqJFv6J580LInNlO6jhERETJxlMxZPoWZ9WcdsmfrhrZ/fsfwtd3v1a9V6+ybGSJiMjo8cwsmbYt9YGoYM1at3vSZDGwhAQVxow5gqlTTwMASpRwg49PSWlDERERpTCemSXT9d8SIOibM5JDhDRZDOz58zDUqOGnbmQBYM+ehxImIiIiSh08M0umaV8X4NYqzZqvSpIohrZ793107rwd79/HAAAsLMwwdWpt+PpWlDgZERFRymMzS6ZFCGBXG+B+gGZ9sMLkx5FVKJQYNeowZsz48hjanDmdsHFjK/zyS3YJkxEREaUeNrNkOiJfAYuzadd7PgPMTPtbPSjoI9q2DcD58y/VtWbNCmLFiiZwcbGRMBkREVHqMu3f8JR+xLwDlufRrg+KBSysDJ/HwEaNOqxuZOVyM/zzTz0MGFAeMhM/G01ERMRmlozfuzuAX2HNmrUL0OetyZ+R/Wzu3Po4ceIprK0t4O/fCmXLZv3xSkRERCYgffymJ9MV9QbYWFmz1mgTUKC1NHkMRKlUaTyxy9XVDnv3dkDOnE5wcrKWMBkREZFhcWguMl6HBwCL3IHYD19qNWaafCO7efMtFC++CCEhURr14sXd2MgSEVG6w2aWjNN/y4Br8zVrXe4AZQZLk8cAYmMT0LfvbrRpE4Dbt0PQufN2qFTpY9xcIiKixPAyAzI+8ZHAwZ6atYZrgYwFpcljAA8evEObNgG4du3L08xcXKwRF5cAGxu5hMmIiIikxWaWjM88B83prveADPmlyWIAGzbcwG+/7UJkZDwAwNraAvPmNUD37qU4WgEREaV7bGbJeITeBFYV06zVXmCyjWxMjAKDBu3D0qVX1LWCBTNh06ZWKFbMTcJkREREaQebWTIOj3cD2xpp10v2NXwWA7h7NxStW2/GzZtv1TUfnxJYsKAh7OwsJUxGRESUtrCZpbTvzgZgT3vNml0WoNdL3cubgPPnX6gbWVtbORYubAgfn5LShiIiIkqD2MxS2iaEdiNb3w8o4iNJHEPx8SmJI0eCcOXKa/j7t0Lhwq5SRyIiIkqT2MxS2hUfAcxz1Kw13wXk/lWaPKkoODgS7u72GrWFCxtCJpPB1pajFRARESWG48xS2vRkn3YjW264yTWyQggsX34FuXPPwZYttzXm2dlZspElIiL6ATazlPZEBQNbG3xTlAHVpkkSJ7VERMShU6dt6NFjJ2JiEtC9eyCCgj5KHYuIiMio8DIDSluEABZl0aw13gzkbyVNnlRy/Xow2rQJwP3779S1du2Kal1qQERERN/HZpbSDkUMsLa0Zq3ZTiCPjiG5jJQQAosXX8bvv+9DXJwSAODgYIlly5qgTZsiEqcjIiIyPmxmKW0Ifw4szaFZ86xvUo1sWFgsfvttFzZtuqWulS6dBZs2tUKePBkkTEZERGS82MyS9BQx2o2sa0mgxW5J4qSGmzffomnTjXj8+IO6NmBAeUyfXhdWVnwbEhERJRd/i5L0/AppThftDngtkyZLKnF2tkZYWKz63ytWNEHz5oV+sBYRERH9CEczIGntageEP/0ynaepyTWyAJA9uyNWr26OChWy4erVXmxkiYiIUgjPzJJ0XpwA7m3UrDXdKk2WFHbp0ivky5cBTk7W6lrDhvlQv35emJnJJExGRERkWnhmlgxPpQRmyQH/6pr1gVGAzLi/JYUQmDnzLCpWXI4ePXZCCKExn40sERFRyjLuzoGM07LcgCpBs9b5P0BuK02eFPLuXTSaNNmIIUMOICFBhYCA29i8+faPVyQiIqJk42UGZFjBl4CIZ5q1nkGAY05J4qSUM2eeo23bADx/Hq6ujRhRGc2bF5QwFRERkeljM0uG8/Y6sK6cZm2I0L2skVCpBKZPP43Ro49Aqfz0tWTKZIs1a5qjfv28EqcjIiIyfWxmyTAU0cCakpq1lvsliZJSQkKi0Lnzduzb91Bdq1YtJ9avb4Fs2RwlTEZERJR+sJklw9hQSXO61SEgZ21psqSAFy/CUaHCMrx6FQEAkMmA0aOrYty4GrCw4KXoREREhsLfupT6Lv4DhFz/Ml39H6NuZAEgWzYHVKiQDQDg5maHAwc64a+/arGRJSIiMrA08Zt3wYIF8PT0hLW1NSpUqIALFy4kuuzSpUtRtWpVuLi4wMXFBXXq1Pnu8iSxrQ2BE8M0a2WHSJMlBclkMixf3gSdO5fAtWu9UadObqkjERERpUuSN7P+/v7w9fXFuHHjcOXKFZQoUQJeXl54+/atzuWPHTuGdu3a4ejRozh79iw8PDxQr149vHz50sDJ6YeODQWe7NWsdbklTZaf9N9/EThy5IlGzcXFBqtWNYO7u71EqYiIiEjyZnbmzJno2bMnunbtisKFC2PRokWwtbXFihUrdC6/bt069O3bFyVLlkTBggWxbNkyqFQqHD582MDJ6buCLwKXZ2jW+ocBGQtLkyeZlEoVJkw4gXHjHqFTpx3qa2SJiIgobZD0BrD4+HhcvnwZo0aNUtfMzMxQp04dnD17NknbiI6OhkKhQIYMGXTOj4uLQ1xcnHo6PPzTOKAKhQIKheIn0ifN530YYl9piXxdeY1pRZ8PgJkNYESvw6tXEfDx2YHjxz+NixsSEo05c85h4sQa0gYjvaXX96Ep4TE0bjx+xs/Qx1Cf/UjazIaGhkKpVMLNzU2j7ubmhrt37yZpGyNGjEDWrFlRp04dnfOnTJmCCRMmaNUPHDgAW1vDPXHq4MGDBtuX1Eq9mYMcX00fz/4PPh48Llme5Lh6NRyzZz9DWNinJ5WZmQHt22fBL79EYc+ePRKno+RKT+9DU8VjaNx4/IyfoY5hdHR0kpc16qG5pk6dio0bN+LYsWOwtrbWucyoUaPg6+urng4PD1dfZ+vomPpjgSoUChw8eBB169aFXC5P9f1Jzez6vzB/eFQ9LTIVQ6UWAyVMpJ+EBBXGjz+Bv/++pq5lzWqP/v3dMWhQ83RxDE1RensfmiIeQ+PG42f8DH0MP3+SnhSSNrOZMmWCubk53rx5o1F/8+YN3N3dv7vuP//8g6lTp+LQoUMoXrx4ostZWVnByspKqy6Xyw36hjL0/iQR9gQ4PkijJOt0FXIzc4kC6efFi3C0a7cFp059edxuw4b5sGzZr7hw4Vj6OIYmjsfQ+PEYGjceP+NnqGOozz4kvQHM0tISZcqU0bh56/PNXBUrVkx0vb///ht//fUX9u3bh7JlyxoiKiXFsm+Gp/rtBWAkjaxCoUT16n7qRtbCwgzTp9fFzp3tkCmT4S5HISIiIv1IPpqBr68vli5dilWrVuHOnTvo06cPoqKi0LVrVwBA586dNW4QmzZtGsaMGYMVK1bA09MTwcHBCA4ORmRkpFRfAgGfRi/4ms8NwCGbNFmSQS43x5Qpnx7kkCOHE06e7IqhQyvBzEwmcTIiIiL6HsmvmfX29kZISAjGjh2L4OBglCxZEvv27VPfFPbs2TOYmX3puf/991/Ex8ejVatWGtsZN24cxo8fb8jo9FlCLLCximYtU1FpsvyENm2KICwsFi1bFkaGDDZSxyEiIqIkkLyZBYD+/fujf//+OucdO3ZMYzooKCj1A1HSCQHM+abx63pPmix62LHjLo4ff4qZM7006j17lpEoERERESVHmmhmyYgFttSc9loJZMgvTZYkiI9XYvjwg5gz5zwAoHTpLOjYMfEbCImIiChtk/yaWTJiF/8BHm77Mu1aEijaRao0P/T48QdUrrxC3cgCwKFDjyVMRERERD+LZ2YpeYIOAieGfVWQAZ2vShbnRwICbqN790CEh396GpylpTlmzfJCnz4cDYOIiMiYsZkl/YX8B2ypp1nrGypNlh+IjU3AkCH7sXDhJXUtb94M2LSpFUqVyiJhMiIiIkoJbGZJf6tLaE73fg3YZJAmy3c8ePAO3t4BuHo1WF1r27YoFi9uBEdH7QdpEBERkfFhM0tJlxCrPXJBs0DA7vtPa5PKyJGH1Y2stbUF5s6tjx49SkMm49ixREREpoLNLCWNriG4CrYD8jSWJk8SLFzYEGfOPIeTkxU2bWqN4sXdpI5EREREKYzNLCVNQF3N6SwVgF/XS5MlEQkJKlhYfBmgw83NHvv3d0Tu3C6wt7eUMBkRERGlFg7NRT927V/g2eEv03buQPtz0uXRYc2a6yhW7F+8exetUS9e3I2NLBERkQljM0vfd3kWcLivZq3XK2my6BAVFY9u3Xagc+ftuHs3FD4+26FSCaljERERkYHwMgNK3Ls7wDFfzVq/D0AauYHq1q23aNMmALdvh6hrbm52UCiUsLLitzYREVF6wN/4pNujncD2Jpq1Xi8Ba2dJ4nxNCIGVK6+hf/89iIlJAADY2cmxaFEjPpqWiIgonWEzS9oe7gB2NNOsNdoE2GeVJM7XIiPj0bv3Lqxbd0NdK17cDf7+rVCwYCYJkxEREZEU2MySJiG0G9mqU4ECrSWJ87Xr14PRpk0A7t9/p6716lUGs2Z5wcZGLmEyIiIikgqbWfpCCGDmN/cEtjkGeFSXJM63Ll16pW5kHRwssXRpY3h7F5U4FREREUmJzSx9sbGq5nSOWmmmkQWAbt1K4ciRINy9Gwp//1bImzftPUKXiIiIDIvNLH1ydDDw6rRmrfVh3csayMuX4ciWzVE9LZPJsGRJI1hYmHG0AiIiIgLAcWYJAK4uAK7M1qz5qiSJAnwarWD+/AvIk2cutm+/qzHPzs6SjSwRERGpsZlN716dBY7016x1uy/ZWLIfP8aidevNGDBgL+LilOjadQeePQuTJAsRERGlfTzFlZ4lxAEbKmnWBkQAlvaSxLlw4SW8vQMQFPRRXevatSTc3aXJQ0RERGkfm9n0SqkA5lhr1noGSdLICiEwe/Y5jBhxCArFp8sbXFys4efXDE2aFDB4HiIiIjIebGbTq9mWmtO1FwCOOQ0e4/37GHTtugOBgffUtYoVs2PDhpbImdPZ4HmIiIjIuLCZTY/+W6I5nb81ULKvwWNcvfoaTZtuxPPn4era8OGVMHFiLcjl5gbPQ0RERMaHzWx6I1TAwV6atcabJImSMaMtIiPj///fNli9ujkaNswnSRYiIiIyThzNID2JegPM/OaM56BYabIAyJHDCatWNUO1ajlx7VpvNrJERESkNzaz6YUQwCJ3zZpbGcDCymARzpx5jvDwOI1a48YFcOyYD7Jnd0xkLSIiIqLEsZlNL475ak7nqAV0vGSQXatUApMmnUDVqivx2287IYTQmC+TaExbIiIiMn5sZtODpZ6aT/jK38pgj6p98yYS9euvxf/+dxQqlYC//y3s2HHvxysSERERJQFvADN1WxsC4U81a438DbLrI0eeoEOHrQgOjgTw6aFi48ZVR+PG+Q2yfyIiIjJ9bGZNmRDAk72atT4hgCx1T8grlSr89dcJ/PnncXy+osDd3R7r17dAzZq5UnXfRERElL6wmTVlJ0ZoTg9WAGape8hfv45Ahw5bcfRokLpWt25urF3bApkz26XqvomIiCj9YTNrqsKCgEvTv0wX/y3VG9mgoI+oUGEZ3r6NAgCYmcnw1181MXJkFZiZ8SYvIiIiSnm8AcwURb8Fln3zcX6team+25w5nfDLL9kBANmyOeDYMR/88UdVNrJERESUatjMmpqoYOBfN81a5+uAuWWq71omk2Hlyqbo3r0Url3rjapVc6b6PomIiCh942UGpkSVACzKolmr9CfgWjxVdrdnzwNYW1ugVq0vZ4EzZLDBsmVNUmV/RERERN/imVlTMttac7rOIqDimBTfjUKhxPDhB/Hrr+vRvv0W9dBbRERERIbGZtZUBB0EhPLLdJEuQIleKb6bZ8/CUL26H6ZPPwMAePMmCkuWXE7x/RARERElBS8zMBXbG335t6UDUH9liu8iMPAeunTZjg8fYgEAcrkZ/v67LgYNqpDi+yIiIiJKCjazpmDGN6MF9HmbopuPj1dixIiDmD37vLrm6emMTZtaoVy5bCm6LyIiIiJ9sJk1dlfmaE5X+xuwsNa9bDI8efIB3t4BuHjxlbrWokUhLF/eBM7OKbcfIiIiouRgM2vM3lwBjv6uWSvjm2Kbj49Xolo1P7x4EQ4AsLQ0x8yZ9dC3bznIZBw7loiIiKTHG8CMlSIGWFtGszYwEjAzT7FdWFqa4++/6wAA8uRxwdmz3dGvX3k2skRERJRm8MyssZprqznd6Sogt0vx3bRrVwzR0Qq0bl0Ejo5WKb59IiIiop/BM7PGaHszzekyg4HMJX96s/7+NzFkyH6tevfupdnIEhERUZrEM7PGRhENPNrxZTpjYaDGzJ/aZEyMAr//vg9LllwBAJQrlw1t2xb9qW0SERERGQLPzBqbbx9X2+XWT23u3r1Q/PLLcnUjCwAnTjz9qW0SERERGQrPzBqT/T2A+PAv0w3W/NTm1q79D71770JUlAIAYGNjgQULGqJLl5I/tV0iIiIiQ2EzaywS4oCby79MZ6kIFO6YrE1FRyswYMAerFhxTV0rXNgVmza1QpEimX8yKBEREZHhsJk1BkoFMOebBxS0O52sTd2+HYLWrTfj9u0Qda1bt5KYN68hbG3lP5OSiIiIyODYzBqD2Zaa07UXAMkc63XkyEPqRtbOTo5///0VnTqV+NmERERERJLgDWBp3eXZmtP5WwEl+yZ7c0uWNEbmzHYoViwzLl36jY0sERERGTWemU3LYt4DxwZr1hpv1msTCoUScvmXp4K5u9vj0KFOyJs3A2xseFkBERERGTeemU3LVuTXnO77LsmrCiGwZMllFCv2L96/j9GYV6yYGxtZIiIiMglsZtOq85OB2K+a19oLAJsMSVo1PDwO7dtvRa9eu3Dv3jt07boDQohUCkpEREQkHV5mkBbFfgROjdasJfE62atXX6NNmwA8fPheXfPwcERCgkrjcgMiIiIiU8BmNi3aXFtzum/oD1cRQmDhwovw9T2A+HglAMDJyQrLlzdBy5aFUyMlERERkeTYzKY1V+YBb788WhZ1FwM2Gb+7ysePsejRIxBbttxR18qVy4qNG1shd26X1EpKREREJDk2s2nJq7PA0YFfph08gOK/fXeVixdfwts7AE+efFTXfv+9AqZNqwtLS15WQERERKaNzWxaEfsR2FBJs+Z9/IerXbnyWt3IurhYw8+vGZo0KZDy+YiIiIjSIDazacWqIprTHS4ATrl+uNpvv5XBkSNBePYsDBs3tkTOnM6pk4+IiIgoDWIzmxZsbQhEvvoyXWYI4F5O56LPn4fBw8NJPS2TybBiRRNYWppztAIiIiJKdzjOrNRengae7P0yna0qUOMfrcVUKoHp008jT5652LXrvsY8OztLNrJERESULrGZldqmGprTOq6TDQ2NRuPGGzB8+CEoFCr4+GzHy5fhhslHRERElIbxMgMpnZsEqBK+THe5BchkGoucPPkU7dptwcuXEQA+ze7duwzc3OwNmZSIiIgoTWIzK6XT//vy72I9gIxfHm6gUglMnXoKY8cehVL56VG0rq62WLu2BerVy2PopERERERpEptZqbw8rTldd4n6n2/fRqFjx604ePCxulajhifWr2+BLFkcDJWQiChdE0IgISEBSqVS6ihGT6FQwMLCArGxsXw9jVRqHEO5XA5z85+/54fNrBQUMcDGKl+mLR3VlxecP/8CzZr5Izg4EsCn8tix1TFmTDWYm/MSZyIiQ4iPj8fr168RHR0tdRSTIISAu7s7nj9/Dtk3l9ORcUiNYyiTyZA9e3bY2//cpZNsZqWwpqTmdOMA9T/d3OwRG/vpOlp3d3usW9cCtWr9eLxZIiJKGSqVCk+ePIG5uTmyZs0KS0tLNmA/SaVSITIyEvb29jAz44kZY5TSx1AIgZCQELx48QL58uX7qTO0bGYNbU9H4MNXQ2vVXgB41lVPeno6Y+XKpli48CLWrGnOG72IiAwsPj4eKpUKHh4esLW1lTqOSVCpVIiPj4e1tTWbWSOVGsfQ1dUVQUFBUCgUbGaNypN9X/5tYYNjHxuiTEQcHBys1OVmzQqiadMCPBNARCQhNl1EqSul+hy+Uw3p9hog9h0AIEFphv+93YlatVahT5/dEEJoLMpGloiIiOjH2MwaSux7YG9nAMDLMAfU8huESZNPQQhg3bob2Lv3ocQBiYiIiIwPm1kDMd/RBACw905elJzZGyfvOH2qm8swbVod1K+fV8p4RERE6da9e/fg7u6OiIgIqaOYjLZt22LGjBkG2VeaaGYXLFgAT09PWFtbo0KFCrhw4cJ3l9+8eTMKFiwIa2trFCtWDHv27DFQ0uTJFH0DyleXMGJXHTRc3hGhUXYAAA8PR5w40RXDh1eGmRkvKyAiouTr0qULZDIZZDIZ5HI5cuXKheHDhyM2NlZr2V27dqF69epwcHCAra0typUrBz8/P53b3bJlC2rUqAEnJyfY29ujePHi+PPPP/H+/ftU/ooMZ9SoURgwYAAcHLTHci9YsCCsrKwQHBysNc/T0xOzZ8/Wqo8fPx4lS5bUqAUHB2PAgAHInTs3rKys4OHhgcaNG+Pw4cMp9WXolJyeKS4uDqNHj0bOnDlhZWUFT09PrFixQj3/1q1baNmyJTw9PSGTyXS+Bv/73/8wadIkhIWFpeSXo5Pkzay/vz98fX0xbtw4XLlyBSVKlICXlxfevn2rc/kzZ86gXbt26N69O65evYpmzZqhWbNmuHnzpoGTJ53HrX9Q498u+PvYl7FlGzfOj6tXe6FSJQ8JkxERkSmpX78+Xr9+jcePH2PWrFlYvHgxxo0bp7HMvHnz0LRpU1SuXBnnz5/Hf//9h7Zt26J3794YOnSoxrKjR4+Gt7c3ypUrh7179+LmzZuYMWMGrl+/jjVr1hjs64qPj0+1bT979gy7du1Cly5dtOadOnUKMTExaNWqFVatWpXsfQQFBaFMmTI4cuQIpk+fjhs3bmDfvn2oWbMm+vXr9xPpvy+5PVObNm1w+PBhLF++HPfu3cOGDRtQoEAB9fzo6Gjkzp0bU6dOhbu7u85tFC1aFHny5MHatWtT9GvSSUisfPnyol+/fupppVIpsmbNKqZMmaJz+TZt2ohff/1Vo1ahQgXRq1evJO0vLCxMABBhYWHJD62HByMzCBebEQIYL4DxQi6fIGbOPCNUKpVB9k8/Lz4+Xmzfvl3Ex8dLHYWSicfQ+BnyGMbExIjbt2+LmJiYVN9XSvLx8RFNmzbVqLVo0UKUKlVKPf3s2TMhl8uFr6+v1vpz584VAMS5c+eEEEKcP39eABCzZ8/Wub8PHz4kmuX58+eibdu2wsXFRdja2ooyZcqIgwcPCqVSqTPnoEGDRPXq1dXT1atXF/369RODBg0SGTNmFDVq1BDt2rUTbdq00VgvPj5eZMyYUaxatUoI8amHmDx5svD09BTW1taiePHiYvPmzYnmFEKI6dOni7Jly+qc16VLFzFy5Eixd+9ekT9/fq35OXPmFLNmzdKqjxs3TpQoUUI93aBBA5EtWzYRGRmptez3XseflZyeae/evcLJyUm8e/dOo65UKsWHDx+EUqnUqCf2GgghxIQJE0SVKlUS3df33mv69GuSDs0VHx+Py5cvY9SoUeqamZkZ6tSpg7Nnz+pc5+zZs/D19dWoeXl5Yfv27TqXj4uLQ1xcnHo6PDwcwKfHsikUip/8Cn4sd4YPqJjzOfbczQ9P11is294b5cplRUJCQqrvm1LG5+8TQ3y/UOrgMTR+hjyGCoUCQgioVCqoVCp1XbauPBCt/VFzqrJ1h+jw/UvvPhNCqHMDwM2bN3HmzBnkzJlTXdu8eTMUCgV8fX01vjYA6NmzJ/744w+sX78e5cqVw9q1a2Fvb4/evXtrLQsAjo6OOuuRkZGoXr06smXLhu3bt8Pd3R1XrlyBSqVSZ/w65+fsADRqq1atQu/evXHy5EkAwMOHD+Ht7Y3w8HD1E6P27t2L6OhoNG3aFCqVCpMnT8a6deuwcOFC5MuXDydOnEDHjh2RMWNGVK9eXefrduLECZQpU0bra4mIiMDmzZtx9uxZFCxYEGFhYTh+/DiqVq2q9bp/u+7XX8/79++xb98+TJw4ETY2NlrLJvY6AsC6devQp08fnfM+2717t1amz86ePYvBgwdrbL9evXrYsWNHovvcsWMHypYti2nTpmHt2rWws7ND48aNMWHChO9+vbq2V7ZsWUyaNAkxMTGwsrLSmv/5e0LXOLP6vNclbWZDQ0OhVCrh5uamUXdzc8Pdu3d1rhMcHKxzeV3XsgDAlClT1AfgawcOHDDIYNgNLOzg1247hh5shXq9myAk5Br27LmW6vullHfw4EGpI9BP4jE0foY4hhYWFnB3d0dkZKTGx9uOka9hFv0q1ff/NZVKqE/C/IhCocDu3bvh6OiIhIQExMXFwczMDNOmTVNv4+bNm3B0dISdnZ3O7ebMmRO3b99GeHg47ty5g5w5cyImJgYxMTFJzuzn54eQkBAcOnQILi4uAD5d/gB8ahAVCgUSEhI09h8fH69RS0hIQO7cuTF69Gj1Mq6urrC1tcX69evRtm1bAMDq1atRv3599dOkpkyZgm3btqF8+fIAgBYtWuDYsWNYsGABSpUqpTPvkydPUKxYMa3XY9WqVcidOzc8PDwQFRWF5s2bY/HixShRooR6GZVKhdjYWK114+LioFQqER4ejuvXr0MIgRw5ciT5WH5Wo0YNnDhx4rvLZMmSJdHtBgcHw8HBQWO+o6MjXr9+neg6Dx48wKlTp2Bubo7Vq1fj3bt3GDp0KIKDg7FgwQKtm+QSew0+7ys+Ph4PHjxAjhw5tObHx8cjJiYGJ06c0DrJp8+jpE3+oQmjRo3SOJMbHv5/7d17VFTl+gfw71yYi8igpjiMIHmDzLRE1IOmpsfCLKXU8KQhJaknNVuaFUtN1PKSeXep3VQ8Hk54WZmuUEjKG1QnM1ATglAoW97SUkBBBub5/eFh/xxg0EEZHP1+1po/5p333fvZ+5mRZ173fqcA/v7+eOKJJ2Aymep8/2V/peNg6ndYM3UIPDw86nx/dPtZrVbs3r0bjz/+OHPopphD9+fKHJaUlODkyZNo2LAhDAaD0q5q6Atx8c26qgbmm/5b5eHhgcceewyrV6/G5cuXsWzZMmi1WrzwwgtKn4qf5nW0TY1GA61WC5PJBI1GA41G4/TfyuzsbHTu3BkBAQFKm4igsLAQXl5e8PDwUPZxfVzXt2m1WnTt2rXKviMiIrBt2zaMHTsWly9fxq5du/Cf//wHJpMJx44dw5UrVzBkyBC7MaWlpejcubPD4ygtLYW3t3eV1xMSEjBq1Cil/aWXXkLfvn2xZs0a5UYxtVoNg8FQZaxer1fOXcXEmdFodPpcmkwmtGjRwqkxlVXer9ForPE9UHETYUJCAry9r628pFarERERgUWLFsHHx8duLXxH5wC49gUEgMP3UUlJCYxGI3r37m33WQPgVOFfr8Vs06ZNodFocPbsWbv2s2fPOryg2Gw2O9Vfr9dXO7Xt4eHhmj9qjQNg1Rxz3f6ozjCH7o85dH+uyGF5eTlUKhXUarX9r4BF/lCn+3XkZstnlUqFhg0bIjAwEACwfv16PPzww1i/fj2io6MBAEFBQbh06RLOnDkDi8ViN760tBTHjx9H3759oVarERQUhLS0NJSXlzt1ziuKt+vPXcV/QatUKuW/k69/vWJW7vq2hg0bVvkVthdeeAF9+vTB+fPnsXv3bhiNRgwcOBBqtVqZyUtMTKxSAOr1eoe/6Na0aVNcvHjR7vXMzEx89913+P777xETE6O0l5eXY/PmzRgzZgyAa8VmQUFBlW1funQJ3t7eynlUqVTIyclx+lfl4uPjMW7cuBr77Nq1y+FlBmazGX/88Yfdfs+dOwez2ewwFovFghYtWiiz6gDQoUMHiAhOnTqF5s2bVxlb8Xmp7OLFiwBQ7RjgWr4rVt+o/B5z5j1Xr6sZ6HQ6dOnSxW5ZCpvNhq+++gqhoaHVjgkNDa2yjMXu3bsd9iciIroXqdVqTJs2DTNmzFAuExg6dCg8PDyqXf/zgw8+wOXLl/H8888DAEaMGIGioiKsXr262u1XFCqVderUCRkZGQ6X7mrWrBlOnz5t15aRkXFTx9SjRw/4+/tj06ZNiI+Px3PPPacUPQ8++CD0ej1+++03tG3b1u7h7+945aDOnTsjMzPTrm3t2rXo3bs3Dh8+jIyMDOUxZcoUrF27VukXFBSEQ4cOVdnmjz/+qHypaNKkCcLCwrBq1Spcvny5Sl9H5xEABg8ebLf/6h4hISEOx9emZurZsydOnTqFoqIipa2iEK/8BehGfvrpJ/j5+aFp06ZOjXPaDW8Rq2MJCQmi1+slLi5OMjMzZezYsdKoUSM5c+aMiIhERkZKTEyM0j8tLU20Wq0sWrRIsrKyJDY2Vjw8POTo0aM3tT9Xr2bAu6jdH3Po/phD98fVDG6sulUCrFartGjRQt5//32lbenSpaJWq2XatGmSlZUlubm5snjxYtHr9fL666/bjX/zzTdFo9HIG2+8Id98843k5+dLSkqKDBs2zOEqB1evXpXAwEDp1auXpKamyvHjx2Xz5s2SnJws5eXlkpSUJCqVSjZs2CA5OTkyc+ZMMZlMVVYzeO2116rd/vTp0+XBBx8UrVYrBw4cqPLafffdJ3FxcZKbmyuHDh2SFStWSFxcnMPztmPHDvHx8ZGysjIRufZea9asmaxZs6ZK38zMTAEgP/30k4hcq0nUarW8++67kpmZKUePHpVp06aJVqu1q0uOHz8uZrNZHnzwQdm6davk5ORIZmamLF++XB544AGHsd2qm6mZYmJiJDIyUnleWFgofn5+MmzYMDl27Jjs27dP2rVrJ9HR0cpqBlevXpX09HRJT08XX19fmTp1qqSnp8svv/xit/+oqCgZPXq0w/hu12oG9V7MioisXLlSWrZsKTqdTrp166YsCyJy7Q0dFRVl13/z5s0SGBgoOp1OOnToIImJiTe9Lxaz5Czm0P0xh+6PxeyNVVfMiojMnz9fmjVrZrcs1Pbt26VXr17i6ekpBoNBunTpIuvWrat2u5s2bZLevXuLl5eXeHp6SqdOnWTOnDk1LimVn58vQ4cOFZPJJA0aNJCQkBBJSUlRlnWaOXOmNG/eXLy9vWXy5MkyceLEmy5mKwrKgICAKstc2mw2WbZsmQQFBYmHh4c0a9ZMwsLCZN++fQ5jtVqtYrFYJCkpSUREtm7dKmq1WplUq6x9+/YyefJk5XlycrL07NlTGjdurCwjVt3+Tp06JRMmTJCAgADR6XTSokULGTx4sOzZs8dhbLfDjWqmqKgou3MvIpKVlSX9+/cXo9Eofn5+MmXKFCkqKlKK2by8PAFQ5XH9doqLi8Xb21u+/fZbh7HdrmJWJfK/9SPuEQUFBfD29salS5dccgOY1WrFzp07MXDgQF6r56aYQ/fHHLo/V+awpKQEeXl5aNWqVZWbUqh2bDYbCgoKYDKZnL5u1BVWrVqFHTt2IDk5ub5DuWM5m8M1a9Zg27Zt+PLLLx32qemz5ky9dtevZkBERERUk3HjxuHixYvKigt06zw8PLBy5UqX7IvFLBEREd3TtFqt3Zq2dOtefvlll+3rzpvrJyIiIiK6SSxmiYiIiMhtsZglIiKqxj12fzSRy92uzxiLWSIioutUrJbgzG/DE5HzSktLAUD5Vbja4g1gRERE19FoNGjUqBHOnTsH4NrPs17/W/TkPJvNhtLSUpSUlNyRS3PRjd3uHNpsNvzxxx9o0KABtNpbK0dZzBIREVViNpsBQClo6daICIqLi2E0GvnFwE3VRQ7VajVatmx5y9tjMUtERFSJSqWCr68vfHx8YLVa6zsct2e1WrF//3707t2bP1zipuoihzqd7rbM8rKYJSIickCj0dzy9Xx07TyWlZXBYDCwmHVTd3IOeeEKEREREbktFrNERERE5LZYzBIRERGR27rnrpmtWKC3oKDAJfuzWq24cuUKCgoK7rhrTOjmMIfujzl0f8yhe2P+3J+rc1hRp93MDyvcc8VsYWEhAMDf37+eIyEiIiKimhQWFsLb27vGPiq5x36vz2az4dSpU/Dy8nLJWncFBQXw9/fHyZMnYTKZ6nx/dPsxh+6POXR/zKF7Y/7cn6tzKCIoLCyExWK54fJd99zMrFqthp+fn8v3azKZ+AF2c8yh+2MO3R9z6N6YP/fnyhzeaEa2Am8AIyIiIiK3xWKWiIiIiNwWi9k6ptfrERsbC71eX9+hUC0xh+6POXR/zKF7Y/7c352cw3vuBjAiIiIiuntwZpaIiIiI3BaLWSIiIiJyWyxmiYiIiMhtsZglIiIiIrfFYvY2WLVqFe6//34YDAZ0794d33//fY39t2zZggceeAAGgwEdO3bEzp07XRQpOeJMDj/++GP06tULjRs3RuPGjdG/f/8b5pzqnrOfwwoJCQlQqVR45pln6jZAuiFnc3jx4kVMmDABvr6+0Ov1CAwM5L+n9cjZ/C1btgxBQUEwGo3w9/fH5MmTUVJS4qJoqbL9+/dj0KBBsFgsUKlU+Pzzz284Zu/evQgODoZer0fbtm0RFxdX53FWS+iWJCQkiE6nk3Xr1smxY8dkzJgx0qhRIzl79my1/dPS0kSj0cjChQslMzNTZsyYIR4eHnL06FEXR04VnM3hiBEjZNWqVZKeni5ZWVny4osvire3t/z+++8ujpwqOJvDCnl5edKiRQvp1auXhIeHuyZYqpazObx69aqEhITIwIEDJTU1VfLy8mTv3r2SkZHh4shJxPn8xcfHi16vl/j4eMnLy5Pk5GTx9fWVyZMnuzhyqrBz506ZPn26fPbZZwJAtm3bVmP/EydOSIMGDWTKlCmSmZkpK1euFI1GI0lJSa4J+DosZm9Rt27dZMKECcrz8vJysVgsMn/+/Gr7R0REyFNPPWXX1r17dxk3blydxkmOOZvDysrKysTLy0s2bNhQVyHSDdQmh2VlZdKjRw/55JNPJCoqisVsPXM2h2vWrJHWrVtLaWmpq0KkGjibvwkTJki/fv3s2qZMmSI9e/as0zjp5txMMfvmm29Khw4d7NqGDx8uYWFhdRhZ9XiZwS0oLS3FoUOH0L9/f6VNrVajf//++Pbbb6sd8+2339r1B4CwsDCH/alu1SaHlV25cgVWqxVNmjSpqzCpBrXN4Zw5c+Dj44Po6GhXhEk1qE0Od+zYgdDQUEyYMAHNmzfHQw89hHnz5qG8vNxVYdP/1CZ/PXr0wKFDh5RLEU6cOIGdO3di4MCBLomZbt2dVM9oXb7Hu8j58+dRXl6O5s2b27U3b94cP//8c7Vjzpw5U23/M2fO1Fmc5FhtcljZW2+9BYvFUuVDTa5RmxympqZi7dq1yMjIcEGEdCO1yeGJEyfw9ddfY+TIkdi5cydyc3Mxfvx4WK1WxMbGuiJs+p/a5G/EiBE4f/48Hn30UYgIysrK8M9//hPTpk1zRch0GziqZwoKClBcXAyj0eiyWDgzS3QLFixYgISEBGzbtg0Gg6G+w6GbUFhYiMjISHz88cdo2rRpfYdDtWSz2eDj44OPPvoIXbp0wfDhwzF9+nR88MEH9R0a3YS9e/di3rx5WL16NX788Ud89tlnSExMxDvvvFPfoZEb4szsLWjatCk0Gg3Onj1r13727FmYzeZqx5jNZqf6U92qTQ4rLFq0CAsWLEBKSgo6depUl2FSDZzN4fHjx5Gfn49BgwYpbTabDQCg1WqRnZ2NNm3a1G3QZKc2n0NfX194eHhAo9Eobe3bt8eZM2dQWloKnU5XpzHT/6tN/t5++21ERkbi5ZdfBgB07NgRly9fxtixYzF9+nSo1Zxru9M5qmdMJpNLZ2UBzszeEp1Ohy5duuCrr75S2mw2G7766iuEhoZWOyY0NNSuPwDs3r3bYX+qW7XJIQAsXLgQ77zzDpKSkhASEuKKUMkBZ3P4wAMP4OjRo8jIyFAegwcPRt++fZGRkQF/f39Xhk+o3eewZ8+eyM3NVb6IAEBOTg58fX1ZyLpYbfJ35cqVKgVrxRcTEam7YOm2uaPqGZffcnaXSUhIEL1eL3FxcZKZmSljx46VRo0ayZkzZ0REJDIyUmJiYpT+aWlpotVqZdGiRZKVlSWxsbFcmqueOZvDBQsWiE6nk61bt8rp06eVR2FhYX0dwj3P2RxWxtUM6p+zOfztt9/Ey8tLJk6cKNnZ2fLFF1+Ij4+PvPvuu/V1CPc0Z/MXGxsrXl5e8umnn8qJEyfkyy+/lDZt2khERER9HcI9r7CwUNLT0yU9PV0AyJIlSyQ9PV1+/fVXERGJiYmRyMhIpX/F0lxvvPGGZGVlyapVq7g0lztbuXKltGzZUnQ6nXTr1k2+++475bU+ffpIVFSUXf/NmzdLYGCg6HQ66dChgyQmJro4YqrMmRwGBAQIgCqP2NhY1wdOCmc/h9djMXtncDaH33zzjXTv3l30er20bt1a5s6dK2VlZS6Omio4kz+r1SqzZs2SNm3aiMFgEH9/fxk/frz89ddfrg+cRERkz5491f5tq8hbVFSU9OnTp8qYRx55RHQ6nbRu3VrWr1/v8rhFRFQinM8nIiIiIvfEa2aJiIiIyG2xmCUiIiIit8ViloiIiIjcFotZIiIiInJbLGaJiIiIyG2xmCUiIiIit8ViloiIiIjcFotZIiIiInJbLGaJiADExcWhUaNG9R1GralUKnz++ec19nnxxRfxzDPPuCQeIiJXYTFLRHeNF198ESqVqsojNze3vkNDXFycEo9arYafnx9eeuklnDt37rZs//Tp03jyyScBAPn5+VCpVMjIyLDrs3z5csTFxd2W/Tkya9Ys5Tg1Gg38/f0xduxY/Pnnn05th4U3Ed0sbX0HQER0Ow0YMADr16+3a2vWrFk9RWPPZDIhOzsbNpsNhw8fxksvvYRTp04hOTn5lrdtNptv2Mfb2/uW93MzOnTogJSUFJSXlyMrKwujR4/GpUuXsGnTJpfsn4juLZyZJaK7il6vh9lstntoNBosWbIEHTt2hKenJ/z9/TF+/HgUFRU53M7hw4fRt29feHl5wWQyoUuXLvjhhx+U11NTU9GrVy8YjUb4+/tj0qRJuHz5co2xqVQqmM1mWCwWPPnkk5g0aRJSUlJQXFwMm82GOXPmwM/PD3q9Ho888giSkpKUsaWlpZg4cSJ8fX1hMBgQEBCA+fPn22274jKDVq1aAQA6d+4MlUqFxx57DID9bOdHH30Ei8UCm81mF2N4eDhGjx6tPN++fTuCg4NhMBjQunVrzJ49G2VlZTUep1arhdlsRosWLdC/f38899xz2L17t/J6eXk5oqOj0apVKxiNRgQFBWH58uXK67NmzcKGDRuwfft2ZZZ37969AICTJ08iIiICjRo1QpMmTRAeHo78/Pwa4yGiuxuLWSK6J6jVaqxYsQLHjh3Dhg0b8PXXX+PNN9902H/kyJHw8/PDwYMHcejQIcTExMDDwwMAcPz4cQwYMABDhw7FkSNHsGnTJqSmpmLixIlOxWQ0GmGz2VBWVobly5dj8eLFWLRoEY4cOYKwsDAMHjwYv/zyCwBgxYoV2LFjBzZv3ozs7GzEx8fj/vvvr3a733//PQAgJSUFp0+fxmeffValz3PPPYcLFy5gz549Stuff/6JpKQkjBw5EgBw4MABjBo1Cq+99hoyMzPx4YcfIi4uDnPnzr3pY8zPz0dycjJ0Op3SZrPZ4Ofnhy1btiAzMxMzZ87EtGnTsHnzZgDA1KlTERERgQEDBuD06dM4ffo0evToAavVirCwMHh5eeHAgQNIS0tDw4YNMWDAAJSWlt50TER0lxEiortEVFSUaDQa8fT0VB7Dhg2rtu+WLVvkvvvuU56vX79evL29ledeXl4SFxdX7djo6GgZO3asXduBAwdErVZLcXFxtWMqbz8nJ0cCAwMlJCREREQsFovMnTvXbkzXrl1l/PjxIiLy6quvSr9+/cRms1W7fQCybds2ERHJy8sTAJKenm7XJyoqSsLDw5Xn4eHhMnr0aOX5hx9+KBaLRcrLy0VE5O9//7vMmzfPbhsbN24UX1/famMQEYmNjRW1Wi2enp5iMBgEgACQJUuWOBwjIjJhwgQZOnSow1gr9h0UFGR3Dq5evSpGo1GSk5Nr3D4R3b14zSwR3VX69u2LNWvWKM89PT0BXJulnD9/Pn7++WcUFBSgrKwMJSUluHLlCho0aFBlO1OmTMHLL7+MjRs3Kv9V3qZNGwDXLkE4cuQI4uPjlf4iApvNhry8PLRv377a2C5duoSGDRvCZrOhpKQEjz76KD755BMUFBTg1KlT6Nmzp13/nj174vDhwwCuXSLw+OOPIygoCAMGDMDTTz+NJ5544pbO1ciRIzFmzBisXr0aer0e8fHx+Mc//gG1Wq0cZ1pamt1MbHl5eY3nDQCCgoKwY8cOlJSU4N///jcyMjLw6quv2vVZtWoV1q1bh99++w3FxcUoLS3FI488UmO8hw8fRm5uLry8vOzaS0pKcPz48VqcASK6G7CYJaK7iqenJ9q2bWvXlp+fj6effhqvvPIK5s6diyZNmiA1NRXR0dEoLS2ttiibNWsWRowYgcTEROzatQuxsbFISEjAs88+i6KiIowbNw6TJk2qMq5ly5YOY/Py8sKPP/4ItVoNX19fGI1GAEBBQcENjys4OBh5eXnYtWsXUlJSEBERgf79+2Pr1q03HOvIoEGDICJITExE165dceDAASxdulR5vaioCLNnz8aQIUOqjDUYDA63q9PplBwsWLAATz31FGbPno133nkHAJCQkICpU6di8eLFCA0NhZeXF95//33897//rTHeoqIidOnSxe5LRIU75SY/InI9FrNEdNc7dOgQbDYbFi9erMw6VlyfWZPAwEAEBgZi8uTJeP7557F+/Xo8++yzCA4ORmZmZpWi+UbUanW1Y0wmEywWC9LS0tCnTx+lPS0tDd26dbPrN3z4cAwfPhzDhg3DgAED8Oeff6JJkyZ226u4PrW8vLzGeAwGA4YMGYL4+Hjk5uYiKCgIwcHByuvBwcHIzs52+jgrmzFjBvr164dXXnlFOc4ePXpg/PjxSp/KM6s6na5K/MHBwdi0aRN8fHxgMpluKSYiunvwBjAiuuu1bdsWVqsVK1euxIkTJ7Bx40Z88MEHDvsXFxdj4sSJ2Lt3L3799VekpaXh4MGDyuUDb731Fr755htMnDgRGRkZ+OWXX7B9+3anbwC73htvvIH33nsPmzZtQnZ2NmJiYpCRkYHXXnsNALBkyRJ8+umn+Pnnn5GTk4MtW7bAbDZX+0MPPj4+MBqNSEpKwtmzZ3Hp0iWH+x05ciQSExOxbt065cavCjNnzsS//vUvzJ49G8eOHUNWVhYSEhIwY8YMp44tNDQUnTp1wrx58wAA7dq1ww8//IDk5GTk5OTg7bffxsGDB+3G3H///Thy5Aiys7Nx/vx5WK1WjBw5Ek2bNkV4eDgOHDiAvLw87N27F5MmTcLvv//uVExEdPdgMUtEd72HH34YS5YswXvvvYeHHnoI8fHxdstaVabRaHDhwgWMGjUKgYGBiIiIwJNPPonZs2cDADp16oR9+/YhJycHvXr1QufOnTFz5kxYLJZaxzhp0iRMmTIFr7/+Ojp27IikpCTs2LED7dq1A3DtEoWFCxciJCQEXbt2RX5+Pnbu3KnMNF9Pq9VixYoV+PDDD2GxWBAeHu5wv/369UOTJk2QnZ2NESNG2L0WFhaGL774Al9++SW6du2Kv/3tb1i6dCkCAgKcPr7Jkyfjk08+wcmTJzFu3DgMGTIEw4cPR/fu3XHhwgW7WVoAGDNmDIKCghASEoJmzZohLS0NDRo0wP79+9GyZUsMGTIE7du3R3R0NEpKSjhTS3QPU4mI1HcQRERERES1wZlZIiIiInJbLGaJiIiIyG2xmCUiIiIit8ViloiIiIjcFotZIiIiInJbLGaJiIiIyG2xmCUiIiIit8ViloiIiIjcFotZIiIiInJbLGaJiIiIyG2xmCUiIiIit/V/+f05AxVclewAAAAASUVORK5CYII=",
"text/plain": [
"<Figure size 800x600 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"\n",
"for sd in [\n",
" # save_df1, \n",
" # save_df2, \n",
" score_df\n",
" ]:\n",
" # 计算二分类指标\n",
" evaluation_metrics = calculate_binary_classification_metrics(sd, score_col='score', label_col='label', threshold=0.6,\n",
" future_return_col='future_return', total_mv_col='total_mv')\n",
"\n",
" # 打印指标\n",
" print(\"二分类评估指标:\")\n",
" for metric, value in evaluation_metrics.items():\n",
" if isinstance(value, (float, int)):\n",
" print(f\"{metric}: {value:.4f}\")\n",
" elif isinstance(value, (list, tuple, np.ndarray)):\n",
" print(f\"{metric}: (array of length {len(value)})\")\n",
" else:\n",
" print(f\"{metric}: {value}\")\n",
"\n",
" # 绘制 ROC 曲线\n",
" plot_roc_curve(evaluation_metrics)"
]
},
{
"cell_type": "code",
"execution_count": 25,
"id": "7e9023cc",
"metadata": {},
"outputs": [],
"source": [
"def analyze_factors(\n",
" df: pd.DataFrame,\n",
" feature_columns: list[str],\n",
" target_column: str = 'target', # 假设目标列默认为 'target'\n",
" trade_date_col: str = 'trade_date', # 假设日期列默认为 'trade_date'\n",
" mcap_col: str = 'total_mv', # 新增: 市值列名称\n",
" mcap_bins: int = 5 # 新增: 市值分位数的数量 (例如 5 表示五分位数)\n",
") -> pd.DataFrame:\n",
" \"\"\"\n",
" 分析DataFrame中指定特征列的各种指标，包括基本统计、相关性、日间IC、ICIR以及在不同市值分位数上的IC。\n",
"\n",
" Args:\n",
" df (pd.DataFrame): 包含日期、目标列、特征列和市值列的DataFrame。\n",
" 需要包含 trade_date_col, target_column, feature_columns 和 mcap_col 中的所有列。\n",
" feature_columns (list[str]): 需要分析的特征列名称列表。\n",
" target_column (str): 目标变量列的名称。\n",
" trade_date_col (str): 交易日期列的名称。\n",
" mcap_col (str): 市值列的名称。\n",
" mcap_bins (int): 市值分位数的数量 (例如 5 表示五分位数)。\n",
"\n",
" Returns:\n",
" pd.DataFrame: 包含各个因子分析指标的汇总DataFrame。\n",
" 同时打印因子在不同市值分位数上的平均IC表格。\n",
" 如果输入数据或列有问题，可能返回空或包含NaN的DataFrame。\n",
" \"\"\"\n",
"\n",
" # --- 数据校验 ---\n",
" required_cols = [trade_date_col, target_column, mcap_col] + feature_columns\n",
" if not all(col in df.columns for col in required_cols):\n",
" missing = [col for col in required_cols if col not in df.columns]\n",
" print(f\"错误: 输入DataFrame缺少必需的列: {missing}\")\n",
" return pd.DataFrame() # 返回空DataFrame\n",
"\n",
" # 确保日期列是 datetime 类型\n",
" df = df.copy() # 在副本上操作\n",
" df[trade_date_col] = pd.to_datetime(df[trade_date_col], errors='coerce')\n",
" df.dropna(subset=[trade_date_col], inplace=True) # 移除日期转换失败的行\n",
"\n",
" # 过滤掉那些在 feature_columns, target_column, mcap_col 上有 NaN 的行，以确保后续计算是在完整数据上\n",
" # 直接在 df 副本上进行清洗\n",
" initial_rows_before_clean = len(df)\n",
" df.dropna(subset=feature_columns + [target_column, mcap_col], inplace=True)\n",
" rows_dropped_clean = initial_rows_before_clean - len(df)\n",
" if rows_dropped_clean > 0:\n",
" print(f\"警告: 移除了 {rows_dropped_clean} 行，因为其特征、目标或市值列存在空值。\")\n",
"\n",
" if df.empty:\n",
" print(\"错误: 清理缺失值后数据为空，无法进行因子分析。\")\n",
" return pd.DataFrame() # 返回空DataFrame\n",
"\n",
"\n",
" print(f\"开始分析 {len(feature_columns)} 个因子指标...\")\n",
"\n",
" # --- 1. 基本因子统计量 ---\n",
" basic_stats = df[feature_columns].describe().T\n",
"\n",
" print(\"\\n--- 基本因子统计量 ---\")\n",
" print(basic_stats)\n",
"\n",
" # --- 2. 因子与目标变量的整体相关性 ---\n",
" overall_correlation = {}\n",
" for feature in feature_columns:\n",
" # 在清理后的 df 上计算相关性\n",
" if df[[feature, target_column]].dropna().shape[0] > 1: # 确保至少有两个有效数据点\n",
" overall_correlation[feature] = {\n",
" 'Pearson_Correlation_with_Target': df[feature].corr(df[target_column], method='pearson'),\n",
" 'Spearman_Correlation_with_Target': df[feature].corr(df[target_column], method='spearman')\n",
" }\n",
" else:\n",
" overall_correlation[feature] = {\n",
" 'Pearson_Correlation_with_Target': np.nan,\n",
" 'Spearman_Correlation_with_Target': np.nan\n",
" }\n",
" overall_corr_df = pd.DataFrame.from_dict(overall_correlation, orient='index')\n",
"\n",
" print(\"\\n--- 因子与目标变量的整体相关性 ---\")\n",
" print(overall_corr_df)\n",
"\n",
" # --- 3. 因子之间的相关性矩阵 ---\n",
" # 在清理后的 df 上计算相关性\n",
" factor_correlation_matrix = df[feature_columns].corr(method='spearman') # 改回 Spearman\n",
"\n",
" print(\"\\n--- 因子之间的相关性矩阵 (Spearman) ---\") # 修正打印信息\n",
" print(factor_correlation_matrix)\n",
"\n",
" # --- 4. 日间 IC 和 ICIR ---\n",
" print(\"\\n--- 计算日间 IC (Spearman 相关性) 和 ICIR ---\")\n",
"\n",
" # 直接在清理后的 df 上计算每日 IC\n",
" if df.empty: # 理论上上面已经检查过，这里再检查一次更安全\n",
" daily_ic_series = pd.Series(dtype=float) # 空 Series\n",
" ic_stats = pd.DataFrame({\n",
" 'Mean_IC (Spearman)': np.nan, 'Std_Dev_IC': np.nan, 'ICIR': np.nan\n",
" }, index=feature_columns)\n",
" else:\n",
" daily_ic_series = df.groupby(trade_date_col).apply(\n",
" lambda day_group: {\n",
" feature: day_group[feature].corr(day_group[target_column], method='spearman')\n",
" for feature in feature_columns if day_group.shape[0] > 1 # 确保每日数据点多于1才能计算相关性\n",
" }\n",
" ).apply(pd.Series) # 将字典结果转换为 DataFrame\n",
"\n",
" # 计算 IC 的统计量\n",
" if not daily_ic_series.empty:\n",
" ic_mean = daily_ic_series.mean()\n",
" ic_std = daily_ic_series.std()\n",
" # 避免除以零\n",
" ic_ir = ic_mean / ic_std.replace(0, np.nan) # 使用 replace 0 为 NaN\n",
"\n",
" ic_stats = pd.DataFrame({\n",
" 'Mean_IC (Spearman)': ic_mean,\n",
" 'Std_Dev_IC': ic_std,\n",
" 'ICIR': ic_ir\n",
" })\n",
" print(\"\\n--- 日间 IC 和 ICIR (Spearman) ---\")\n",
" print(ic_stats)\n",
" else:\n",
" ic_stats = pd.DataFrame({\n",
" 'Mean_IC (Spearman)': np.nan, 'Std_Dev_IC': np.nan, 'ICIR': np.nan\n",
" }, index=feature_columns)\n",
"\n",
"\n",
" # --- 5. 因子在不同市值分位数上的平均 IC ---\n",
" print(f\"\\n--- 计算因子在 {mcap_bins} 个市值分位数上的平均 IC (Spearman) ---\")\n",
"\n",
" # 在清理后的 df 上计算每日市值分位数，直接添加到 df 中\n",
" # 使用 transform() 和 qcut() 在每个日期分组内计算分位数\n",
" # labels=False 返回整数 0 to mcap_bins-1\n",
" # duplicates='drop' 处理在某些日期股票数量少于 bins 导致分位数边缘重复的情况，会返回 NaN\n",
" # 添加一个临时列来存储分位数\n",
" mcap_bin_col_name = f'_mcap_bin_{mcap_bins}'\n",
" df[mcap_bin_col_name] = df.groupby(trade_date_col)[mcap_col].transform(\n",
" lambda x: pd.qcut(x, q=mcap_bins, labels=False, duplicates='drop') if len(x) >= mcap_bins else np.nan # 确保股票数量足够进行分位数划分\n",
" )\n",
"\n",
" # 过滤掉无法划分分位数 (NaN) 的行，进行分位数 IC 计算\n",
" # 创建一个临时 DataFrame df_binned_analysis\n",
" df_binned_analysis = df.dropna(subset=[mcap_bin_col_name]).copy()\n",
"\n",
" if df_binned_analysis.empty:\n",
" print(\"错误: 划分市值分位数后数据为空，无法计算分位数上的 IC。\")\n",
" avg_ic_by_bin = pd.DataFrame(index=range(mcap_bins), columns=feature_columns) # Placeholder\n",
" else:\n",
" # 按日期和市值分位数分组，计算每个分组内的因子与目标变量的截面相关性 (分位数IC)\n",
" binned_ic_by_day = df_binned_analysis.groupby([trade_date_col, mcap_bin_col_name]).apply(\n",
" lambda group: {\n",
" feature: group[feature].corr(group[target_column], method='spearman')\n",
" for feature in feature_columns if group.shape[0] > 1 # 确保分位数组内数据点多于1\n",
" }\n",
" ).apply(pd.Series) # 将嵌套结果转为 DataFrame\n",
"\n",
" # 对每个分位数组的每日 IC 求平均\n",
" # unstack(level=mcap_bin_col_name) 将 mcap_bin 作为列\n",
" # mean(axis=0) 对日期索引求平均\n",
" avg_ic_by_bin = binned_ic_by_day.unstack(level=mcap_bin_col_name).mean(axis=0).unstack()\n",
"\n",
" # 重命名索引和列，使表格更清晰\n",
" if not avg_ic_by_bin.empty:\n",
" # Index name will be the original column name used for grouping ('_mcap_bin_X')\n",
" # Rename the index name explicitly\n",
" avg_ic_by_bin.index.name = 'MarketCap_Bin'\n",
" avg_ic_by_bin.columns.name = 'Feature'\n",
" # 可以根据需要对分位数 bin 索引进行排序 (虽然 pd.qcut labels=False usually sorts)\n",
" avg_ic_by_bin = avg_ic_by_bin.sort_index()\n",
"\n",
" print(avg_ic_by_bin)\n",
"\n",
"\n",
" # --- 6. 汇总所有指标 ---\n",
" # 将基本统计、整体相关性、IC/ICIR 合并到一个 DataFrame\n",
" # 注意：合并时需要根据索引进行对齐 (因子名称)\n",
" summary_df = basic_stats\n",
" summary_df = summary_df.merge(overall_corr_df, left_index=True, right_index=True, how='left')\n",
" summary_df = summary_df.merge(ic_stats, left_index=True, right_index=True, how='left')\n",
"\n",
" # print(\"\\n--- 因子分析汇总报告 ---\")\n",
" # print(summary_df)\n",
"\n",
" # --- 清理临时列 'mcap_bin' ---\n",
" # 修正：在函数结束时从我们一直在操作的 df 副本中删除临时列\n",
" if mcap_bin_col_name in df.columns:\n",
" df.drop(columns=[mcap_bin_col_name], inplace=True)\n",
"\n",
"\n",
" return summary_df # 主要返回汇总报告，分位数IC单独打印\n",
"\n",
"# # 运行分析函数\n",
"# factor_analysis_report = analyze_factors(test_data.copy(), feature_columns, 'future_return')\n",
"\n",
"# print(\"\\n--- 最终汇总报告 DataFrame ---\")\n",
"# print(factor_analysis_report)"
]
},
{
"cell_type": "code",
"execution_count": 26,
"id": "a0000d75",
"metadata": {},
"outputs": [
{
"ename": "ModuleNotFoundError",
"evalue": "No module named 'seaborn'",
"output_type": "error",
"traceback": [
"\u001b[31m---------------------------------------------------------------------------\u001b[39m",
"\u001b[31mModuleNotFoundError\u001b[39m Traceback (most recent call last)",
"\u001b[36mCell\u001b[39m\u001b[36m \u001b[39m\u001b[32mIn[26]\u001b[39m\u001b[32m, line 2\u001b[39m\n\u001b[32m 1\u001b[39m \u001b[38;5;28;01mimport\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34;01mmatplotlib\u001b[39;00m\u001b[34;01m.\u001b[39;00m\u001b[34;01mpyplot\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;28;01mas\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34;01mplt\u001b[39;00m\n\u001b[32m----> \u001b[39m\u001b[32m2\u001b[39m \u001b[38;5;28;01mimport\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34;01mseaborn\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;28;01mas\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34;01msns\u001b[39;00m\n\u001b[32m 3\u001b[39m \u001b[38;5;28;01mfrom\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34;01mscipy\u001b[39;00m\u001b[34;01m.\u001b[39;00m\u001b[34;01mstats\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;28;01mimport\u001b[39;00m spearmanr\n\u001b[32m 4\u001b[39m \u001b[38;5;28;01mfrom\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[34;01mtqdm\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;28;01mimport\u001b[39;00m tqdm \u001b[38;5;66;03m# 用于显示进度条 (可选)\u001b[39;00m\n",
"\u001b[31mModuleNotFoundError\u001b[39m: No module named 'seaborn'"
]
}
],
"source": [
"import matplotlib.pyplot as plt\n",
"import seaborn as sns\n",
"from scipy.stats import spearmanr\n",
"from tqdm import tqdm # 用于显示进度条 (可选)\n",
"\n",
"# 设置 Matplotlib/Seaborn 样式 (可选)\n",
"sns.set_theme(style=\"whitegrid\")\n",
"plt.rcParams['font.sans-serif'] = ['SimHei'] # 或者其他支持中文的字体\n",
"plt.rcParams['axes.unicode_minus'] = False # 解决负号显示问题\n",
"\n",
"def analyze_score_performance_2d(score_df: pd.DataFrame,\n",
" score_col: str = 'score',\n",
" label_col: str = 'label',\n",
" condition1_col: str = 'circ_mv',\n",
" condition2_col: str = 'future_return',\n",
" n_bins: int = 100,\n",
" min_samples_per_bin: int = 30): # 每个格子最少样本数\n",
" \"\"\"\n",
" 分析 score 在两个条件下 (如市值、未来收益) 的二维分箱表现。\n",
"\n",
" Args:\n",
" score_df (pd.DataFrame): 包含分数、标签和条件列的 DataFrame。\n",
" score_col (str): 预测分数所在的列名。\n",
" label_col (str): 目标标签所在的列名 (应为数值或可排序类别)。\n",
" condition1_col (str): 第一个条件列名 (例如 'circ_mv')。\n",
" condition2_col (str): 第二个条件列名 (例如 'future_return')。\n",
" n_bins (int): 每个条件划分的箱数 (分位数数量)。\n",
" min_samples_per_bin (int): 计算指标所需的最小样本数，小于此数目的格子结果将被屏蔽。\n",
"\n",
" Returns:\n",
" tuple: 包含 (performance_pivot, count_pivot, fig)\n",
" performance_pivot: 以二维分箱为索引/列的 Spearman 相关系数矩阵。\n",
" count_pivot: 每个二维分箱的样本数量矩阵。\n",
" fig: 生成的热力图 Matplotlib Figure 对象。\n",
" \"\"\"\n",
" print(f\"开始分析 '{score_col}' 在 '{condition1_col}' 和 '{condition2_col}' 下的表现...\")\n",
"\n",
" required_cols = [score_col, label_col, condition1_col, condition2_col]\n",
" if not all(col in score_df.columns for col in required_cols):\n",
" missing = [col for col in required_cols if col not in score_df.columns]\n",
" raise ValueError(f\"输入 DataFrame 缺少必需列: {missing}\")\n",
"\n",
" # --- 1. 数据准备和清洗 ---\n",
" print(\"准备数据，处理 NaN 值...\")\n",
" # 只保留需要的列，并移除包含 NaN 的行，避免影响分箱和计算\n",
" analysis_df = score_df[required_cols].dropna().copy()\n",
" n_original = len(score_df)\n",
" n_after_drop = len(analysis_df)\n",
" print(f\"原始数据 {n_original} 行，移除 NaN 后剩余 {n_after_drop} 行用于分析。\")\n",
"\n",
" if n_after_drop < min_samples_per_bin * n_bins: # 检查数据量是否过少\n",
" print(f\"警告: 清理 NaN 后数据量 ({n_after_drop}) 可能不足以支持 {n_bins}x{n_bins} 的精细分箱分析。\")\n",
" if n_after_drop < min_samples_per_bin:\n",
" print(\"错误: 有效数据过少，无法进行分析。\")\n",
" return None, None, None\n",
"\n",
" # --- 2. 二维分箱 ---\n",
" print(f\"对 '{condition1_col}' 和 '{condition2_col}' 进行 {n_bins} 分位数分箱...\")\n",
" bin1_col = f'{condition1_col}_bin'\n",
" bin2_col = f'{condition2_col}_bin'\n",
"\n",
" try:\n",
" # 使用 qcut 进行分位数分箱，labels=False 返回 0 到 n_bins-1 的整数标签\n",
" # duplicates='drop' 会丢弃导致边界不唯一的重复值所在的箱子，可能导致某些箱号缺失\n",
" # 对于可视化，这通常可以接受，但如果需要严格的等分，需先 rank\n",
" analysis_df[bin1_col] = pd.qcut(analysis_df[condition1_col], q=n_bins, labels=False, duplicates='drop')\n",
" analysis_df[bin2_col] = pd.qcut(analysis_df[condition2_col], q=n_bins, labels=False, duplicates='drop')\n",
" except Exception as e:\n",
" print(f\"错误: 分箱失败，请检查数据分布或减少 n_bins。错误信息: {e}\")\n",
" # 可以尝试先 rank 再 qcut\n",
" # analysis_df[bin1_col] = pd.qcut(analysis_df[condition1_col].rank(method='first'), q=n_bins, labels=False, duplicates='raise')\n",
" # analysis_df[bin2_col] = pd.qcut(analysis_df[condition2_col].rank(method='first'), q=n_bins, labels=False, duplicates='raise')\n",
" return None, None, None\n",
"\n",
" # --- 3. 分组计算表现指标 (Spearman Rank IC) ---\n",
" print(\"按二维分箱分组计算 Spearman Rank IC...\")\n",
"\n",
" def safe_spearmanr(x, y):\n",
" \"\"\"安全计算 Spearman 相关性，处理数据量过少的情况\"\"\"\n",
" if len(x) < max(2, min_samples_per_bin): # 要求至少有 min_samples_per_bin 个点才计算\n",
" return np.nan\n",
" corr, p_value = spearmanr(x, y)\n",
" return corr if not np.isnan(corr) else np.nan # 确保返回 NaN 而不是 None 或其他\n",
"\n",
" # 按两个分箱列分组\n",
" grouped = analysis_df.groupby([bin1_col, bin2_col])\n",
"\n",
" # 计算每个格子的 Spearman 相关系数\n",
" # apply 可能较慢，但计算相关性通常需要 apply\n",
" performance_series = grouped.apply(lambda sub: safe_spearmanr(sub[score_col], sub[label_col]))\n",
"\n",
" # 计算每个格子的样本数量\n",
" count_series = grouped.size()\n",
"\n",
" # --- 4. 结果整理成 Pivot Table (用于绘图) ---\n",
" print(\"整理结果用于绘图...\")\n",
" try:\n",
" # 将 performance_series 转换成二维矩阵\n",
" # index 为 condition1_bin, columns 为 condition2_bin\n",
" performance_pivot = performance_series.unstack(level=0) # level=0 对应第一个 groupby key (bin1_col)\n",
" count_pivot = count_series.unstack(level=0)\n",
"\n",
" # 可选：按列和索引排序，确保顺序正确\n",
" performance_pivot = performance_pivot.sort_index(axis=0).sort_index(axis=1)\n",
" count_pivot = count_pivot.sort_index(axis=0).sort_index(axis=1)\n",
" \n",
" print(performance_pivot)\n",
"\n",
" except Exception as e:\n",
" print(f\"错误: 无法将结果转换为二维矩阵，可能因为分箱不均匀或数据问题: {e}\")\n",
" return None, None, None\n",
"\n",
" # --- 5. 可视化：绘制热力图 ---\n",
" print(\"生成热力图...\")\n",
" fig, ax = plt.subplots(figsize=(16, 12)) # 调整图像大小\n",
"\n",
" # 使用 count_pivot 创建一个 mask，屏蔽掉样本量过小的格子\n",
" mask = count_pivot < min_samples_per_bin\n",
"\n",
" # 绘制热力图\n",
" sns.heatmap(performance_pivot,\n",
" annot=False, # 100x100 个格子加注释会太密集\n",
" fmt=\".2f\",\n",
" cmap=\"viridis\", # 选择颜色映射, 'viridis', 'coolwarm', 'RdYlGn' 等都不错\n",
" linewidths=.5,\n",
" linecolor='lightgray',\n",
" # mask=mask, # 应用 mask\n",
" ax=ax,\n",
" cbar_kws={'label': f'Spearman Rank IC ({score_col} vs {label_col})'}) # 颜色条标签\n",
"\n",
" # 设置标题和轴标签\n",
" ax.set_title(f'{score_col} 表现分析 (Rank IC vs {label_col})\\n基于 {condition1_col} 和 {condition2_col} {n_bins}x{n_bins} 分箱', fontsize=16)\n",
" ax.set_xlabel(f'{condition1_col} 分位数 (0 -> 高)', fontsize=12)\n",
" ax.set_ylabel(f'{condition2_col} 分位数 (0 -> 高)', fontsize=12)\n",
"\n",
" # 可选：调整刻度标签，避免显示所有 100 个刻度\n",
" if n_bins > 20:\n",
" tick_interval = n_bins // 10 # 大约显示 10 个刻度\n",
" ax.set_xticks(np.arange(0, n_bins, tick_interval) + 0.5)\n",
" ax.set_yticks(np.arange(0, n_bins, tick_interval) + 0.5)\n",
" ax.set_xticklabels(np.arange(0, n_bins, tick_interval))\n",
" ax.set_yticklabels(np.arange(0, n_bins, tick_interval))\n",
"\n",
" plt.xticks(rotation=45, ha='right')\n",
" plt.yticks(rotation=0)\n",
" plt.tight_layout() # 调整布局\n",
"\n",
" print(\"分析完成。\")\n",
" return performance_pivot, count_pivot, fig\n",
"\n",
"# --- 如何使用 ---\n",
"# 假设你的包含预测结果和所需列的 DataFrame 是 final_predictions_df\n",
"# 确保它包含 'score', 'label', 'circ_mv', 'future_return'\n",
"\n",
"# # 示例调用 (你需要有实际的 score_df)\n",
"try:\n",
" # 确保数据类型正确\n",
" cols_to_numeric = ['score', 'label', 'circ_mv', 'future_return']\n",
" for col in cols_to_numeric:\n",
" if col in score_df.columns:\n",
" score_df[col] = pd.to_numeric(score_df[col], errors='coerce')\n",
"\n",
" # 调用分析函数\n",
" performance_matrix, count_matrix, heatmap_figure = analyze_score_performance_2d(\n",
" score_df,\n",
" n_bins=100, # 你要求的100分箱\n",
" min_samples_per_bin=50 # 每个格子至少需要50个样本才显示IC，可以调整\n",
" )\n",
"\n",
" # 显示图像\n",
" if heatmap_figure:\n",
" plt.show()\n",
"\n",
" # 可以查看具体的 performance_matrix 和 count_matrix\n",
" # print(\"\\nPerformance Matrix (Spearman IC):\")\n",
" # print(performance_matrix)\n",
" # print(\"\\nCount Matrix:\")\n",
" # print(count_matrix)\n",
"\n",
"except ValueError as ve:\n",
" print(f\"数据错误: {ve}\")\n",
"except Exception as e:\n",
" print(f\"发生未知错误: {e}\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a436dba4",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Empty DataFrame\n",
"Columns: [ts_code, trade_date, is_st]\n",
"Index: []\n"
]
}
],
"source": [
"print(df[(df['ts_code'] == '600242.SH') & (df['trade_date'] >= '2023-06-01')][['ts_code', 'trade_date', 'is_st']])"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "stock",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.11"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Reference in New Issue View Git Blame Copy Permalink