多因子选股实战：从因子构建到组合优化的完整流程

import pandas as pd
import numpy as np
from scipy import stats

def construct_factor(df, factor_name, method='zscore'):
    """
    因子构建函数
    df: 包含因子原始数据的DataFrame
    factor_name: 因子名称
    method: 标准化方法 ('zscore', 'rank', 'normalize')
    """
    factor_raw = df[factor_name]
    
    if method == 'zscore':
        # Z-score标准化
        factor_processed = (factor_raw - factor_raw.mean()) / factor_raw.std()
    elif method == 'rank':
        # 排名标准化
        factor_processed = factor_raw.rank(pct=True) * 2 - 1
    elif method == 'normalize':
        # 0-1归一化
        factor_processed = (factor_raw - factor_raw.min()) / (factor_raw.max() - factor_raw.min())
    
    # 去极值（Winsorize）
    factor_processed = winsorize(factor_processed, limits=0.025)
    
    # 中性化（行业、市值）
    factor_processed = neutralize(factor_processed, df['industry'], df['market_cap'])
    
    return factor_processed

def winsorize(series, limits=0.025):
    """去极值处理"""
    lower = series.quantile(limits)
    upper = series.quantile(1 - limits)
    return series.clip(lower, upper)

def neutralize(factor, industry_dummies, market_cap):
    """因子中性化"""
    # 行业中性化
    X = pd.concat([industry_dummies, market_cap], axis=1)
    model = sm.OLS(factor, X).fit()
    factor_neutral = model.resid
    return factor_neutral

def analyze_factor_decay(factor_data, return_data, max_lag=20):
    """
    因子衰减分析
    """
    ic_series = []
    for lag in range(1, max_lag + 1):
        ic = factor_data.corrwith(return_data.shift(-lag), axis=0)
        ic_series.append(ic.mean())
    
    # 绘制IC衰减曲线
    plt.figure(figsize=(10, 6))
    plt.plot(range(1, max_lag + 1), ic_series, marker='o')
    plt.xlabel('预测周期(天)')
    plt.ylabel('IC值')
    plt.title('因子IC衰减曲线')
    plt.grid(True)
    plt.savefig('factor_decay.png')
    
    # 计算半衰期
    half_life = np.argmin(np.array(ic_series) > ic_series[0] / 2)
    print(f"因子半衰期: {half_life}天")
    
    return ic_series

def ic_ir_weighted_synthesis(factor_matrix, return_data, lookback=252):
    """
    IC_IR动态加权合成
    """
    weights = []
    for factor_name in factor_matrix.columns:
        factor = factor_matrix[factor_name]
        ic_series = [factor.shift(i).corr(return_data) for i in range(1, 21)]
        ic_mean = np.mean(ic_series)
        ic_std = np.std(ic_series)
        ic_ir = ic_mean / ic_std if ic_std > 0 else 0
        weights.append(abs(ic_ir))
    
    weights = np.array(weights) / sum(weights)
    composite_factor = (factor_matrix * weights).sum(axis=1)
    return composite_factor, weights

import cvxpy as cp

def mean_variance_optimization(expected_return, cov_matrix, risk_aversion=1.0):
    """
    均值-方差优化
    expected_return: 预期收益率(N,)
    cov_matrix: 协方差矩阵(N, N)
    risk_aversion: 风险厌恶系数
    """
    n = len(expected_return)
    w = cp.Variable(n)
    
    # 目标函数：最大化效用 = 收益 - 风险厌恶 * 方差
    utility = expected_return @ w - risk_aversion * cp.quad_form(w, cov_matrix)
    
    # 约束条件
    constraints = [
        cp.sum(w) == 1,  # 全额投资
        w >= 0  # 不允许做空
    ]
    
    problem = cp.Problem(cp.Maximize(utility), constraints)
    problem.solve()
    
    return w.value

def risk_parity_optimization(cov_matrix):
    """
    风险平价优化：每个资产贡献相同风险
    """
    n = cov_matrix.shape[0]
    w = cp.Variable(n)
    
    # 风险贡献
    portfolio_vol = cp.sqrt(cp.quad_form(w, cov_matrix))
    risk_contribution = (w @ cov_matrix) / portfolio_vol
    
    # 目标：风险贡献相等
    target_risk = cp.Variable()
    objective = cp.Minimize(cp.sum_squares(risk_contribution - target_risk))
    
    constraints = [
        cp.sum(w) == 1,
        w >= 0
    ]
    
    problem = cp.Problem(objective, constraints)
    problem.solve()
    
    return w.value

def optimization_with_constraints(expected_return, cov_matrix, industry_matrix, 
                                  max_weight=0.05, target_industry_weight=None):
    """
    带行业中性约束的优化
    industry_matrix: 行业哑变量矩阵
    """
    n = len(expected_return)
    w = cp.Variable(n)
    
    # 目标函数
    utility = expected_return @ w - 0.5 * cp.quad_form(w, cov_matrix)
    
    # 约束条件
    constraints = [
        cp.sum(w) == 1,
        w >= 0,
        w <= max_weight,  # 个股权重上限
    ]
    
    # 行业中性约束
    if target_industry_weight is not None:
        industry_exposure = industry_matrix.T @ w
        constraints.append(industry_exposure == target_industry_weight)
    
    problem = cp.Problem(cp.Maximize(utility), constraints)
    problem.solve()
    
    return w.value

# 回测引擎核心代码
class MultiFactorBacktest:
    def __init__(self, factor_data, price_data, benchmark_data):
        self.factor_data = factor_data
        self.price_data = price_data
        self.benchmark_data = benchmark_data
        
    def run_backtest(self, factor_weights, rebalance_freq='M'):
        """
        执行回测
        """
        returns = []
        weights_history = []
        
        for date in self.get_rebalance_dates(rebalance_freq):
            # 1. 因子合成
            composite_factor = (self.factor_data.loc[date] * factor_weights).sum(axis=1)
            
            # 2. 选股（Top 30）
            selected_stocks = composite_factor.nlargest(30).index
            
            # 3. 组合优化
            expected_return = self.predict_return(composite_factor, selected_stocks)
            cov_matrix = self.estimate_covariance(selected_stocks, date)
            optimal_weights = mean_variance_optimization(expected_return, cov_matrix)
            
            # 4. 计算收益
            portfolio_return = (optimal_weights * self.price_data[selected_stocks].pct_change()).sum()
            returns.append(portfolio_return)
            weights_history.append(optimal_weights)
        
        return pd.Series(returns, index=self.get_rebalance_dates(rebalance_freq))

class RiskMonitor:
    def __init__(self, portfolio_weights, factor_exposures):
        self.weights = portfolio_weights
        self.factor_exposures = factor_exposures
        
    def check_risk_limits(self):
        """风险限额检查"""
        risks = {}
        
        # 1. 行业集中度
        industry_concentration = self.calculate_industry_concentration()
        risks['industry'] = industry_concentration
        
        # 2. 因子暴露
        factor_exposure = self.factor_exposures.T @ self.weights
        risks['factor'] = factor_exposure
        
        # 3. 换手率
        turnover = self.calculate_turnover()
        risks['turnover'] = turnover
        
        # 4. VaR
        var_95 = self.calculate_var(confidence=0.95)
        risks['var'] = var_95
        
        return risks
    
    def generate_alert(self, risks, limits):
        """生成风险预警"""
        alerts = []
        for risk_type, value in risks.items():
            if abs(value) > limits[risk_type]:
                alerts.append(f"{risk_type}超限: {value:.4f} > {limits[risk_type]:.4f}")
        return alerts