""" 因果与风险动力学模块 Causal & Risk Dynamics Layer (FIN 552, FIN 553, ECON 610) """ import pandas as pd import numpy as np from statsmodels.tsa.api import VAR from statsmodels.tsa.stattools import grangercausalitytests, adfuller from scipy.linalg import cholesky, eig from scipy.stats import norm, t, chi2, jarque_bera from scipy.optimize import minimize from sklearn.covariance import LedoitWolf, EmpiricalCovariance from sklearn.ensemble import RandomForestRegressor from sklearn.metrics import mutual_info_score import networkx as nx from typing import Dict, List, Optional, Tuple, Union, Callable import warnings warnings.filterwarnings('ignore') try: from nolitsa import dimension, lyapunov, delay NOLITSA_AVAILABLE = True except ImportError: NOLITSA_AVAILABLE = False print("Warning: nolitsa not available. Some chaos analysis features will be disabled.") class CausalDiscovery: """因果关系发现器""" def __init__(self): self.causal_graph = None self.var_results = None self.granger_results = {} def discover_causal_structure(self, data: pd.DataFrame, maxlags: int = 5, method: str = 'var') -> Dict[str, any]: """ 发现变量间的因果结构 Args: data: 输入数据 maxlags: 最大滞后阶数 method: 方法选择 ('var', 'granger', 'mutual_info', 'random_forest') Returns: Dict: 因果结构分析结果 """ results = {} if method == 'var': results = self._var_causal_analysis(data, maxlags) elif method == 'granger': results = self._granger_causality_analysis(data, maxlags) elif method == 'mutual_info': results = self._mutual_info_causality(data, maxlags) elif method == 'random_forest': results = self._random_forest_causality(data, maxlags) else: raise ValueError(f"Unknown method: {method}") # 构建因果图 self.causal_graph = self._build_causal_graph(results) return results def _var_causal_analysis(self, data: pd.DataFrame, maxlags: int) -> Dict[str, any]: """VAR模型因果分析""" # 数据预处理 data_clean = data.dropna() # 检查平稳性 stationarity_results = {} for col in data_clean.columns: adf_result = adfuller(data_clean[col]) stationarity_results[col] = { 'adf_statistic': adf_result[0], 'p_value': adf_result[1], 'is_stationary': adf_result[1] < 0.05 } # 拟合VAR模型 try: model = VAR(data_clean) var_results = model.fit(maxlags=maxlags, ic='aic') self.var_results = var_results # 提取因果关系 causal_matrix = self._extract_var_causality(var_results) # 脉冲响应分析 irf = var_results.irf(periods=10) # 方差分解 fevd = var_results.fevd(periods=10) return { 'type': 'var', 'model_results': var_results, 'causal_matrix': causal_matrix, 'stationarity': stationarity_results, 'irf': irf, 'fevd': fevd, 'aic': var_results.aic, 'bic': var_results.bic, 'optimal_lags': var_results.k_ar } except Exception as e: print(f"VAR analysis failed: {str(e)}") return {'type': 'var', 'error': str(e)} def _granger_causality_analysis(self, data: pd.DataFrame, maxlags: int) -> Dict[str, any]: """格兰杰因果性分析""" results = {'type': 'granger', 'pairwise_results': {}} variables = data.columns.tolist() causal_matrix = pd.DataFrame(0, index=variables, columns=variables) for i, cause in enumerate(variables): for j, effect in enumerate(variables): if i != j: try: # 构造双变量数据 test_data = data[[effect, cause]].dropna() if len(test_data) > maxlags * 3: # 确保有足够数据 gc_test = grangercausalitytests(test_data, maxlags, verbose=False) # 提取最优滞后的p值 p_values = [gc_test[lag][0]['ssr_ftest'][1] for lag in range(1, maxlags + 1)] min_p_value = min(p_values) optimal_lag = p_values.index(min_p_value) + 1 results['pairwise_results'][f"{cause}->{effect}"] = { 'p_value': min_p_value, 'optimal_lag': optimal_lag, 'is_causal': min_p_value < 0.05 } causal_matrix.loc[cause, effect] = 1 if min_p_value < 0.05 else 0 except Exception as e: print(f"Granger test failed for {cause}->{effect}: {str(e)}") results['causal_matrix'] = causal_matrix return results def _mutual_info_causality(self, data: pd.DataFrame, maxlags: int) -> Dict[str, any]: """基于互信息的因果分析""" variables = data.columns.tolist() causal_matrix = pd.DataFrame(0.0, index=variables, columns=variables) for cause in variables: for effect in variables: if cause != effect: # 计算不同滞后的互信息 mi_scores = [] for lag in range(1, maxlags + 1): cause_lagged = data[cause].shift(lag).dropna() effect_current = data[effect].iloc[lag:].dropna() if len(cause_lagged) == len(effect_current) and len(cause_lagged) > 0: # 离散化连续变量 cause_discrete = pd.cut(cause_lagged, bins=10, labels=False) effect_discrete = pd.cut(effect_current, bins=10, labels=False) # 计算互信息 mi = mutual_info_score(cause_discrete, effect_discrete) mi_scores.append(mi) if mi_scores: causal_matrix.loc[cause, effect] = max(mi_scores) return { 'type': 'mutual_info', 'causal_matrix': causal_matrix, 'threshold': causal_matrix.values.mean() + causal_matrix.values.std() } def _random_forest_causality(self, data: pd.DataFrame, maxlags: int) -> Dict[str, any]: """基于随机森林的因果分析""" variables = data.columns.tolist() causal_matrix = pd.DataFrame(0.0, index=variables, columns=variables) for target in variables: # 构造特征矩阵 feature_data = [] feature_names = [] for var in variables: if var != target: for lag in range(1, maxlags + 1): feature_data.append(data[var].shift(lag)) feature_names.append(f"{var}_lag_{lag}") if feature_data: X = pd.concat(feature_data, axis=1).dropna() y = data[target].iloc[maxlags:].dropna() if len(X) == len(y) and len(X) > 50: # 确保有足够数据 rf = RandomForestRegressor(n_estimators=100, random_state=42) rf.fit(X, y) # 提取特征重要性 importances = rf.feature_importances_ # 聚合每个变量的重要性 for var in variables: if var != target: var_importances = [importances[i] for i, name in enumerate(feature_names) if name.startswith(f"{var}_lag_")] if var_importances: causal_matrix.loc[var, target] = sum(var_importances) return { 'type': 'random_forest', 'causal_matrix': causal_matrix, 'threshold': causal_matrix.values.mean() + causal_matrix.values.std() } def _extract_var_causality(self, var_results) -> pd.DataFrame: """从VAR结果中提取因果关系""" variables = var_results.names n_vars = len(variables) causal_matrix = pd.DataFrame(0, index=variables, columns=variables) # 使用系数的t统计量判断因果关系 for i, effect in enumerate(variables): for j, cause in enumerate(variables): if i != j: # 检查所有滞后期的系数 causality_strength = 0 for lag in range(var_results.k_ar): coef = var_results.coefs[lag][i, j] # 从cause到effect的系数 stderr = var_results.stderr[lag][i, j] t_stat = abs(coef / (stderr + 1e-8)) if t_stat > 1.96: # 95%置信度 causality_strength += t_stat causal_matrix.loc[cause, effect] = causality_strength return causal_matrix def _build_causal_graph(self, causal_results: Dict) -> nx.DiGraph: """构建因果关系图""" G = nx.DiGraph() if 'causal_matrix' in causal_results: causal_matrix = causal_results['causal_matrix'] threshold = causal_results.get('threshold', causal_matrix.values.mean()) # 添加节点 G.add_nodes_from(causal_matrix.index) # 添加边 for cause in causal_matrix.index: for effect in causal_matrix.columns: if cause != effect and causal_matrix.loc[cause, effect] > threshold: G.add_edge(cause, effect, weight=causal_matrix.loc[cause, effect]) return G class GradientDynamicsSimulator: """梯度动力学模拟器""" def __init__(self): self.trajectory_history = [] def simulate_gradient_dynamics(self, b0: np.ndarray, potential_func: Callable, eta: float = 0.1, sigma: float = 0.05, dt: float = 0.01, steps: int = 1000, method: str = 'euler') -> Dict[str, any]: """ 模拟梯度动力学: db_t = -η ∇U(b_t) dt + σ dW_t Args: b0: 初始状态 potential_func: 势函数 eta: 学习率/衰减参数 sigma: 噪声强度 dt: 时间步长 steps: 模拟步数 method: 数值方法 ('euler', 'heun', 'rk4') Returns: Dict: 模拟结果 """ trajectory = [b0.copy()] current_b = b0.copy() potential_values = [potential_func(b0)] for step in range(steps): if method == 'euler': new_b = self._euler_step(current_b, potential_func, eta, sigma, dt) elif method == 'heun': new_b = self._heun_step(current_b, potential_func, eta, sigma, dt) elif method == 'rk4': new_b = self._rk4_step(current_b, potential_func, eta, sigma, dt) else: raise ValueError(f"Unknown method: {method}") trajectory.append(new_b.copy()) potential_values.append(potential_func(new_b)) current_b = new_b # 检查收敛 if step > 100 and np.linalg.norm(trajectory[-1] - trajectory[-2]) < 1e-6: print(f"Converged at step {step}") break trajectory = np.array(trajectory) self.trajectory_history.append(trajectory) return { 'trajectory': trajectory, 'potential_values': np.array(potential_values), 'final_state': trajectory[-1], 'convergence_time': len(trajectory), 'path_length': self._calculate_path_length(trajectory) } def _euler_step(self, b: np.ndarray, potential_func: Callable, eta: float, sigma: float, dt: float) -> np.ndarray: """欧拉方法步进""" grad_U = self._numerical_gradient(potential_func, b) drift = -eta * grad_U diffusion = sigma * np.random.normal(0, 1, size=b.shape) return b + drift * dt + diffusion * np.sqrt(dt) def _heun_step(self, b: np.ndarray, potential_func: Callable, eta: float, sigma: float, dt: float) -> np.ndarray: """Heun方法步进""" # 预测步 grad_U1 = self._numerical_gradient(potential_func, b) drift1 = -eta * grad_U1 noise = np.random.normal(0, 1, size=b.shape) b_pred = b + drift1 * dt + sigma * noise * np.sqrt(dt) # 修正步 grad_U2 = self._numerical_gradient(potential_func, b_pred) drift2 = -eta * grad_U2 drift_avg = (drift1 + drift2) / 2 return b + drift_avg * dt + sigma * noise * np.sqrt(dt) def _rk4_step(self, b: np.ndarray, potential_func: Callable, eta: float, sigma: float, dt: float) -> np.ndarray: """四阶Runge-Kutta方法""" noise = np.random.normal(0, 1, size=b.shape) k1 = -eta * self._numerical_gradient(potential_func, b) k2 = -eta * self._numerical_gradient(potential_func, b + 0.5 * dt * k1) k3 = -eta * self._numerical_gradient(potential_func, b + 0.5 * dt * k2) k4 = -eta * self._numerical_gradient(potential_func, b + dt * k3) drift = (k1 + 2*k2 + 2*k3 + k4) / 6 diffusion = sigma * noise * np.sqrt(dt) return b + drift * dt + diffusion def _numerical_gradient(self, func: Callable, x: np.ndarray, eps: float = 1e-6) -> np.ndarray: """数值计算梯度""" grad = np.zeros_like(x) for i in range(len(x)): x_plus = x.copy() x_plus[i] += eps x_minus = x.copy() x_minus[i] -= eps grad[i] = (func(x_plus) - func(x_minus)) / (2 * eps) return grad def _calculate_path_length(self, trajectory: np.ndarray) -> float: """计算轨迹长度""" diffs = np.diff(trajectory, axis=0) distances = np.linalg.norm(diffs, axis=1) return np.sum(distances) def create_vrp_potential(self, vix_target: float = 0.2, rv_target: float = 0.15) -> Callable: """创建VRP势函数""" def vrp_potential(b): """U(b) = (VIX² - RV)² / 2 + regularization""" if len(b) < 2: return 0.5 * (b[0] - vix_target)**2 vix_sq = b[0] rv = b[1] vrp = vix_sq - rv # 主要项：VRP平方 main_term = 0.5 * vrp**2 # 正则化项：防止极值 reg_term = 0.01 * (vix_sq**2 + rv**2) # 约束项：确保非负 constraint_term = 0.1 * (max(0, -vix_sq)**2 + max(0, -rv)**2) return main_term + reg_term + constraint_term return vrp_potential def analyze_attractor_dynamics(self, trajectory: np.ndarray) -> Dict[str, any]: """分析吸引子动力学""" results = {} # 计算平衡点 final_points = trajectory[-100:] # 最后100个点 equilibrium = np.mean(final_points, axis=0) equilibrium_stability = np.std(final_points, axis=0) results['equilibrium'] = equilibrium results['equilibrium_stability'] = equilibrium_stability # 计算回归时间 distances_to_eq = np.linalg.norm(trajectory - equilibrium, axis=1) results['return_times'] = self._calculate_return_times(distances_to_eq) # 相空间分析 if trajectory.shape[1] >= 2: results['phase_space_area'] = self._calculate_phase_space_area(trajectory) results['limit_cycle_detection'] = self._detect_limit_cycles(trajectory) return results def _calculate_return_times(self, distances: np.ndarray, threshold: float = None) -> List[int]: """计算回归时间""" if threshold is None: threshold = np.std(distances) return_times = [] last_return = 0 for i, dist in enumerate(distances): if dist < threshold and i - last_return > 10: # 避免重复计数 return_times.append(i - last_return) last_return = i return return_times def _calculate_phase_space_area(self, trajectory: np.ndarray) -> float: """计算相空间面积""" if trajectory.shape[1] < 2: return 0.0 x, y = trajectory[:, 0], trajectory[:, 1] # 使用鞋带公式计算多边形面积 return 0.5 * np.abs(np.dot(x, np.roll(y, 1)) - np.dot(y, np.roll(x, 1))) def _detect_limit_cycles(self, trajectory: np.ndarray, tolerance: float = 0.1) -> Dict[str, any]: """检测极限环""" if trajectory.shape[0] < 100: return {'has_limit_cycle': False} # 简单的极限环检测：检查是否存在周期性 mid_point = len(trajectory) // 2 recent_traj = trajectory[mid_point:] # 计算自相关来检测周期性 if len(recent_traj) > 50: x_centered = recent_traj[:, 0] - np.mean(recent_traj[:, 0]) autocorr = np.correlate(x_centered, x_centered, mode='full') autocorr = autocorr[autocorr.size // 2:] autocorr = autocorr / autocorr[0] # 标准化 # 寻找主要周期 peaks = [] for i in range(10, min(len(autocorr), 200)): if (autocorr[i] > autocorr[i-1] and autocorr[i] > autocorr[i+1] and autocorr[i] > 0.3): peaks.append((i, autocorr[i])) has_cycle = len(peaks) > 0 main_period = peaks[0][0] if peaks else 0 return { 'has_limit_cycle': has_cycle, 'main_period': main_period, 'cycle_strength': peaks[0][1] if peaks else 0 } return {'has_limit_cycle': False} class RiskMetrics: """风险计量工具""" def __init__(self): self.risk_models = {} def calculate_var_es(self, returns: pd.Series, confidence_levels: List[float] = [0.95, 0.99], methods: List[str] = ['historical', 'parametric', 'cornish_fisher']) -> Dict[str, Dict]: """ 计算VaR和ES Args: returns: 收益率序列 confidence_levels: 置信水平 methods: 计算方法 Returns: Dict: VaR和ES结果 """ results = {} for method in methods: method_results = {} for alpha in confidence_levels: var, es = self._calculate_single_var_es(returns, alpha, method) method_results[f'VaR_{int(alpha*100)}'] = var method_results[f'ES_{int(alpha*100)}'] = es results[method] = method_results return results def _calculate_single_var_es(self, returns: pd.Series, alpha: float, method: str) -> Tuple[float, float]: """计算单个VaR和ES""" returns_clean = returns.dropna() if method == 'historical': var = np.percentile(returns_clean, (1 - alpha) * 100) es = returns_clean[returns_clean <= var].mean() elif method == 'parametric': mu = returns_clean.mean() sigma = returns_clean.std() var = norm.ppf(1 - alpha, mu, sigma) es = mu - sigma * norm.pdf(norm.ppf(1 - alpha)) / (1 - alpha) elif method == 't_distribution': # 拟合t分布 df, mu, sigma = t.fit(returns_clean) var = t.ppf(1 - alpha, df, mu, sigma) es = self._calculate_t_es(returns_clean, alpha, df, mu, sigma) elif method == 'cornish_fisher': # Cornish-Fisher展开 mu = returns_clean.mean() sigma = returns_clean.std() skew = returns_clean.skew() kurt = returns_clean.kurtosis() z_alpha = norm.ppf(1 - alpha) cf_adjustment = (z_alpha**2 - 1) * skew / 6 + \ (z_alpha**3 - 3*z_alpha) * kurt / 24 - \ (2*z_alpha**3 - 5*z_alpha) * skew**2 / 36 z_cf = z_alpha + cf_adjustment var = mu + sigma * z_cf es = self._estimate_es_cornish_fisher(returns_clean, var, alpha) else: raise ValueError(f"Unknown method: {method}") return var, es def _calculate_t_es(self, returns: pd.Series, alpha: float, df: float, mu: float, sigma: float) -> float: """计算t分布的期望短缺""" t_val = t.ppf(1 - alpha, df) es = mu - sigma * (t.pdf(t_val, df) / (1 - alpha)) * ((df + t_val**2) / (df - 1)) return es def _estimate_es_cornish_fisher(self, returns: pd.Series, var: float, alpha: float) -> float: """估计Cornish-Fisher展开的ES""" tail_returns = returns[returns <= var] if len(tail_returns) > 0: return tail_returns.mean() else: return var * 1.2 # 简单估计 def calculate_risk_decomposition(self, portfolio_returns: pd.Series, component_returns: pd.DataFrame) -> Dict[str, any]: """计算风险分解""" results = {} # 总体组合风险 portfolio_var = portfolio_returns.var() results['total_variance'] = portfolio_var # 成分风险贡献 component_contributions = {} component_weights = {} for component in component_returns.columns: # 计算相关系数 corr = portfolio_returns.corr(component_returns[component]) component_var = component_returns[component].var() # 风险贡献 = 权重 × 相关系数 × 成分标准差 × 组合标准差 contribution = corr * np.sqrt(component_var * portfolio_var) component_contributions[component] = contribution # 计算隐含权重 weight = contribution / portfolio_var if portfolio_var > 0 else 0 component_weights[component] = weight results['component_contributions'] = component_contributions results['component_weights'] = component_weights results['diversification_ratio'] = sum(component_contributions.values()) / np.sqrt(portfolio_var) return results def perform_stress_testing(self, returns: pd.Series, scenarios: Dict[str, Dict]) -> Dict[str, Dict]: """执行压力测试""" stress_results = {} for scenario_name, scenario_params in scenarios.items(): scenario_results = {} if scenario_params['type'] == 'historical': # 历史情景 start_date = scenario_params['start_date'] end_date = scenario_params['end_date'] historical_returns = returns.loc[start_date:end_date] if isinstance(returns.index, pd.DatetimeIndex) else returns scenario_results['total_return'] = (1 + historical_returns).prod() - 1 scenario_results['max_drawdown'] = self._calculate_max_drawdown(historical_returns) scenario_results['var_95'] = np.percentile(historical_returns, 5) elif scenario_params['type'] == 'monte_carlo': # 蒙特卡罗情景 n_simulations = scenario_params.get('n_simulations', 10000) shock_size = scenario_params.get('shock_size', -0.1) # 生成随机冲击 random_shocks = np.random.normal(shock_size, np.abs(shock_size) * 0.5, n_simulations) scenario_results['simulated_returns'] = random_shocks scenario_results['var_95'] = np.percentile(random_shocks, 5) scenario_results['worst_case'] = np.min(random_shocks) elif scenario_params['type'] == 'parametric': # 参数化冲击 shock_magnitude = scenario_params.get('shock_magnitude', -3) # 3个标准差 mu = returns.mean() sigma = returns.std() shocked_return = mu + shock_magnitude * sigma scenario_results['shocked_return'] = shocked_return scenario_results['probability'] = norm.cdf(shock_magnitude) stress_results[scenario_name] = scenario_results return stress_results def _calculate_max_drawdown(self, returns: pd.Series) -> float: """计算最大回撤""" cumulative = (1 + returns).cumprod() running_max = cumulative.expanding().max() drawdown = (cumulative - running_max) / running_max return drawdown.min() class ChaosAnalysis: """混沌分析工具""" def __init__(self): self.embedding_params = {} def analyze_chaos_indicators(self, series: pd.Series, max_dim: int = 10) -> Dict[str, any]: """分析混沌指标""" if not NOLITSA_AVAILABLE: return {'error': 'nolitsa library not available'} results = {} series_array = series.dropna().values if len(series_array) < 100: return {'error': 'Insufficient data for chaos analysis'} try: # 估计嵌入维数 embedding_dimension = self._estimate_embedding_dimension(series_array, max_dim) results['embedding_dimension'] = embedding_dimension # 估计延迟时间 delay_time = self._estimate_delay_time(series_array) results['delay_time'] = delay_time # 重构相空间 embedded = dimension.embed(series_array, dim=embedding_dimension, tau=delay_time) results['phase_space_dimension'] = embedded.shape # 计算最大Lyapunov指数 lle = self._calculate_lyapunov_exponent(embedded) results['lyapunov_exponent'] = lle results['is_chaotic'] = lle > 0.01 # 计算关联维数 correlation_dimension = self._calculate_correlation_dimension(embedded) results['correlation_dimension'] = correlation_dimension # 计算Hurst指数 hurst_exponent = self._calculate_hurst_exponent(series_array) results['hurst_exponent'] = hurst_exponent results['persistence_type'] = self._interpret_hurst(hurst_exponent) except Exception as e: results['error'] = f"Chaos analysis failed: {str(e)}" return results def _estimate_embedding_dimension(self, series: np.ndarray, max_dim: int) -> int: """估计嵌入维数""" try: # 使用false nearest neighbors方法 dimensions = range(1, max_dim + 1) fnn_percentages = [] for dim in dimensions: embedded = dimension.embed(series, dim=dim, tau=1) if len(embedded) > 50: # 简化的false nearest neighbors计算 fnn_pct = self._calculate_fnn_percentage(embedded, dim) fnn_percentages.append(fnn_pct) else: fnn_percentages.append(100) # 数据不足 # 找到FNN百分比显著下降的点 optimal_dim = 1 for i in range(1, len(fnn_percentages)): if fnn_percentages[i] < 10 and fnn_percentages[i-1] - fnn_percentages[i] > 20: optimal_dim = i + 1 break return min(optimal_dim, max_dim) except: return 3 # 默认值 def _estimate_delay_time(self, series: np.ndarray) -> int: """估计延迟时间""" try: # 使用自相关函数的第一个零点 autocorr = np.correlate(series, series, mode='full') autocorr = autocorr[autocorr.size // 2:] autocorr = autocorr / autocorr[0] # 找到第一个接近零的点 for i in range(1, min(len(autocorr), 50)): if autocorr[i] <= 0.1: return i return 1 # 默认值 except: return 1 def _calculate_fnn_percentage(self, embedded: np.ndarray, dim: int) -> float: """计算虚假最近邻百分比""" if len(embedded) < 100: return 100 # 简化实现 n_points = min(len(embedded), 500) sample_indices = np.random.choice(len(embedded), n_points, replace=False) false_neighbors = 0 for i in sample_indices[:100]: # 限制计算量 point = embedded[i] distances = np.linalg.norm(embedded - point, axis=1) nearest_idx = np.argsort(distances)[1] # 第二小的距离（排除自己） # 简单的虚假邻居判断 if distances[nearest_idx] > np.std(embedded.flatten()): false_neighbors += 1 return (false_neighbors / min(100, len(sample_indices))) * 100 def _calculate_lyapunov_exponent(self, embedded: np.ndarray) -> float: """计算最大Lyapunov指数""" try: # 使用nolitsa计算 lle = lyapunov.mle(embedded, maxt=50, window=10) return float(lle) if not np.isnan(lle) else 0.0 except: # 简化实现 return self._simple_lyapunov_estimate(embedded) def _simple_lyapunov_estimate(self, embedded: np.ndarray) -> float: """简化的Lyapunov指数估计""" if len(embedded) < 50: return 0.0 # 计算相邻点的平均发散速度 divergences = [] for i in range(len(embedded) - 10): initial_separation = np.linalg.norm(embedded[i+1] - embedded[i]) if initial_separation > 1e-10: final_separation = np.linalg.norm(embedded[i+10] - embedded[i+9]) if final_separation > 1e-10: divergence = np.log(final_separation / initial_separation) / 10 if not np.isnan(divergence) and np.isfinite(divergence): divergences.append(divergence) return np.mean(divergences) if divergences else 0.0 def _calculate_correlation_dimension(self, embedded: np.ndarray) -> float: """计算关联维数""" try: n_points = min(len(embedded), 1000) sample_data = embedded[np.random.choice(len(embedded), n_points, replace=False)] # 计算不同半径下的关联积分 radii = np.logspace(-3, 0, 20) correlation_integrals = [] for r in radii: distances = np.linalg.norm(sample_data[:, np.newaxis] - sample_data, axis=2) correlation_sum = np.sum(distances < r) - n_points # 排除对角线 correlation_integral = correlation_sum / (n_points * (n_points - 1)) correlation_integrals.append(correlation_integral + 1e-10) # 使用线性回归估计维数 log_r = np.log(radii) log_c = np.log(correlation_integrals) # 选择线性区域 valid_indices = np.isfinite(log_c) & np.isfinite(log_r) if np.sum(valid_indices) > 5: slope, _ = np.polyfit(log_r[valid_indices], log_c[valid_indices], 1) return max(0, min(slope, 10)) # 限制在合理范围内 return 1.0 except: return 1.0 def _calculate_hurst_exponent(self, series: np.ndarray) -> float: """计算Hurst指数""" try: n = len(series) if n < 50: return 0.5 # R/S分析 max_k = int(n / 4) rs_values = [] time_scales = [] for k in range(10, max_k, max(1, max_k // 20)): # 将序列分割成长度为k的子序列 n_subseries = n // k rs_subseries = [] for i in range(n_subseries): subseries = series[i*k:(i+1)*k] mean_sub = np.mean(subseries) # 计算累积离差 cumulative_deviations = np.cumsum(subseries - mean_sub) # 计算R和S R = np.max(cumulative_deviations) - np.min(cumulative_deviations) S = np.std(subseries) if S > 1e-10: rs_subseries.append(R / S) if rs_subseries: rs_values.append(np.mean(rs_subseries)) time_scales.append(k) if len(rs_values) >= 5: # 对数线性回归 log_rs = np.log(rs_values) log_scales = np.log(time_scales) valid_indices = np.isfinite(log_rs) & np.isfinite(log_scales) if np.sum(valid_indices) >= 5: hurst, _ = np.polyfit(log_scales[valid_indices], log_rs[valid_indices], 1) return max(0, min(1, hurst)) # 限制在[0,1]范围内 return 0.5 except: return 0.5 def _interpret_hurst(self, hurst: float) -> str: """解释Hurst指数""" if hurst < 0.4: return "Mean-reverting (Anti-persistent)" elif hurst > 0.6: return "Trending (Persistent)" else: return "Random walk" # 示例和测试 if __name__ == "__main__": # 创建示例数据 np.random.seed(42) dates = pd.date_range('2023-01-01', periods=1000, freq='D') # 生成模拟的金融时间序列 returns = np.random.normal(0.001, 0.02, 1000) prices = 100 * np.exp(np.cumsum(returns)) # 创建多变量数据 data = pd.DataFrame({ 'Price': prices, 'Volume': np.random.lognormal(10, 0.5, 1000), 'Returns': returns, 'Volatility': pd.Series(returns).rolling(20).std().fillna(0.02) }, index=dates) print("Testing Causal & Risk Dynamics...") # 1. 因果发现测试 causal_discovery = CausalDiscovery() causal_results = causal_discovery.discover_causal_structure(data, method='var') print(f"Causal analysis type: {causal_results.get('type', 'failed')}") if 'causal_matrix' in causal_results: print(f"Causal matrix shape: {causal_results['causal_matrix'].shape}") # 2. 梯度动力学模拟测试 simulator = GradientDynamicsSimulator() vrp_potential = simulator.create_vrp_potential() initial_state = np.array([0.25, 0.15]) # VIX^2, RV dynamics_results = simulator.simulate_gradient_dynamics( initial_state, vrp_potential, steps=500 ) print(f"Dynamics simulation final state: {dynamics_results['final_state']}") print(f"Convergence time: {dynamics_results['convergence_time']}") # 3. 风险计量测试 risk_metrics = RiskMetrics() returns_series = pd.Series(returns, index=dates) var_es_results = risk_metrics.calculate_var_es(returns_series) print("VaR/ES Results:") for method, results in var_es_results.items(): print(f" {method}: VaR_95 = {results['VaR_95']:.4f}, ES_95 = {results['ES_95']:.4f}") # 4. 混沌分析测试 chaos_analyzer = ChaosAnalysis() chaos_results = chaos_analyzer.analyze_chaos_indicators(returns_series) if 'error' not in chaos_results: print("Chaos Analysis Results:") print(f" Lyapunov exponent: {chaos_results.get('lyapunov_exponent', 'N/A'):.4f}") print(f" Hurst exponent: {chaos_results.get('hurst_exponent', 'N/A'):.4f}") print(f" Is chaotic: {chaos_results.get('is_chaotic', False)}") print(f" Persistence type: {chaos_results.get('persistence_type', 'N/A')}") else: print(f"Chaos analysis error: {chaos_results['error']}")