""" Portfolio Modeling Module for RiskOptimix This module implements a portfolio modeling using the CCC (Constant Conditional Correlation) and DCC (Dynamic Conditional Correlation) GARCH ensemble approaches for portfolio analysis. It is made to work with user-submitted portfolio and provide risk modeling results. The paper that this module is based on and where all mathematical models are described can be found at www.riskoptimix.com/approach For any questions or feedback, please contact riskoptimix@gmail.com. """ # Import used libraries import numpy as np import pandas as pd import yfinance as yf from arch import arch_model from scipy.optimize import minimize from scipy import stats import warnings import json from datetime import datetime import logging from typing import List, Dict, Tuple, Optional import gc import matplotlib.pyplot as plt import os # Filter out all warnings warnings.filterwarnings('ignore') # Set up logging, this is used to log the progress of the portfolio modeling logging.basicConfig(level=logging.INFO) logger = logging.getLogger(__name__) # Set the number of days between rebalancing the individual GARCH models # This here set as a constant as this number is frequently changed when testing the models REBALANCE_FREQUENCY = 50 # Number of days between rebalancing # This is the main class that implements the portfolio modeling pipeline of the CCC-GARCH ensemble model class PortfolioCCCEnsembleModel: """ This class contains the CCC-GARCH ensemble model for portfolio analysis. This class implements the complete modeling pipeline for user-submitted portfolios: 1. Data obtaining for portfolio stocks 2. Ensemble GARCH model fitting for each stock 3. CCC correlation matrix estimation 4. Portfolio optimization and risk assessment 5. Results generation """ def __init__(self, portfolio_data, start_date='2014-01-01', end_date=None, confidence_level=0.05, train_ratio=0.8): """ Initialize the portfolio modeling system Args: portfolio_data: List of (ticker, amount) tuples start_date: Start date for historical data end_date: End date for historical data (None = today) confidence_level: VaR confidence level train_ratio: Training/testing split ratio """ # Set the initial parameters for the portfolio modeling self.portfolio_data = portfolio_data self.start_date = start_date self.end_date = end_date or datetime.now().strftime('%Y-%m-%d') self.confidence_level = confidence_level self.train_ratio = train_ratio # Portfolio composition self.tickers = [ticker for ticker, amount in portfolio_data] self.amounts = [float(amount) for ticker, amount in portfolio_data] self.total_value = sum(self.amounts) self.weights = [amount / self.total_value for amount in self.amounts] # Data storage self.returns = None self.prices = None self.ensemble_models = {} self.conditional_volatilities = {} self.correlation_matrix = None self.standardized_residuals = None # Results storage self.modeling_results = {} self.optimization_results = {} logger.info(f"Initialized portfolio model with {len(self.tickers)} stocks") logger.info(f"Portfolio value: ${self.total_value:,.2f}") logger.info(f"Data period: {self.start_date} to {self.end_date}") def download_portfolio_data(self) -> bool: """ Download historical data for all portfolio stocks Returns: True if the data is downloaded successfully, False otherwise """ # Log the start of the data downloading logger.info("Downloading portfolio data...") try: # Download data for all tickers using yfinance data = yf.download(self.tickers, start=self.start_date, end=self.end_date) if len(self.tickers) == 1: # Handle single stock case - data['Close'] is already a Series if isinstance(data['Close'], pd.Series): prices = data['Close'].to_frame() prices.columns = self.tickers else: # If it's already a DataFrame prices = data['Close'] if prices.shape[1] == 1: prices.columns = self.tickers else: # Multiple stocks - data['Close'] should be a DataFrame prices = data['Close'] if hasattr(prices, 'columns'): prices.columns = self.tickers # Calculate log returns returns = np.log(prices / prices.shift(1)).dropna() # Store data self.prices = prices self.returns = returns logger.info(f"Downloaded {len(returns)} observations") logger.info(f"Date range: {returns.index[0]} to {returns.index[-1]}") # Print basic statistics for ticker in self.tickers: mean_ret = returns[ticker].mean() std_ret = returns[ticker].std() logger.info(f"{ticker}: Mean={mean_ret:.4f}, Std={std_ret:.4f}") return True except Exception as e: logger.error(f"Error downloading data: {str(e)}") return False def create_garch_model_specifications(self) -> List[Dict]: """ Create the GARCH model specifications for the individual stocks in the portfolio. Returns: List of GARCH model specifications """ # Initialize the list of GARCH model specifications models = [] # Define the GARCH model configurations # In the ensemble, we use GARCH, EGARCH and TARCH models with different parameters model_configs = [ ('GARCH', [(1,1), (1,2), (2,1), (2,2)]), ('EGARCH', [(1,1), (1,2), (2,1)]), ('TGARCH', [(1,1), (1,2), (2,1)]) ] # For each model, we use three different distributions: normal, t and ged distributions = ['normal', 't', 'ged'] # We iterate over all the model configurations and distributions for model_type, param_combos in model_configs: for p, q in param_combos: for dist in distributions: # Create the model specification for the current model configuration and distribution model_spec = { 'model_type': model_type, 'p': p, 'q': q, 'distribution': dist, 'name': f"{model_type}({p},{q})_{dist}" } # Set the volatility type for the current model configuration if model_type == 'EGARCH': model_spec['vol_type'] = 'EGARCH' elif model_type == 'TGARCH': model_spec['vol_type'] = 'GARCH' model_spec['power'] = 2.0 else: model_spec['vol_type'] = 'Garch' # Add the model specification to the list of model specifications models.append(model_spec) return models def fit_garch_model(self, returns, model_spec, verbose=False) -> Optional[arch_model]: """ Fit the GARCH model for a given stock and model specification. Args: returns: The returns of the stock model_spec: The model specification verbose: Whether to log the progress Returns: The fitted GARCH model """ # Try to fit the GARCH model try: # If the model type is GARCH, we use the GARCH model if model_spec['model_type'] == 'GARCH': model = arch_model( returns, vol='Garch', p=model_spec['p'], q=model_spec['q'], dist=model_spec['distribution'] ) # If the model type is EGARCH, we use the EGARCH model elif model_spec['model_type'] == 'EGARCH': model = arch_model( returns, vol='EGARCH', p=model_spec['p'], q=model_spec['q'], dist=model_spec['distribution'] ) # If the model type is TGARCH, we use the TGARCH model elif model_spec['model_type'] == 'TGARCH': model = arch_model( returns, vol='GARCH', p=model_spec['p'], q=model_spec['q'], power=model_spec.get('power', 2.0), dist=model_spec['distribution'] ) # Fit the model result = model.fit(disp='off', show_warning=False) if verbose: logger.info(f"Successfully fitted {model_spec['name']}") return result except Exception as e: if verbose: logger.error(f"Failed to fit {model_spec['name']}: {str(e)}") return None def calculate_var(self, returns_forecast, volatility_forecast, confidence_level=0.05, distribution='normal', dist_params=None) -> float: """ Calculate Value at Risk using GARCH forecasts Args: returns_forecast: The forecasted returns volatility_forecast: The forecasted volatility confidence_level: The confidence level for the VaR distribution: The distribution of the returns dist_params: The distribution parameters Returns: The calculated VaR """ # Calculate the VaR using the GARCH forecasts if distribution == 'normal': z_score = stats.norm.ppf(confidence_level) elif distribution == 't': if dist_params is not None and 'nu' in dist_params: nu = dist_params['nu'] z_score = stats.t.ppf(confidence_level, df=nu) else: z_score = stats.norm.ppf(confidence_level) elif distribution == 'ged': if dist_params is not None and 'nu' in dist_params: nu = dist_params['nu'] z_score = stats.gennorm.ppf(confidence_level, beta=nu) else: z_score = stats.norm.ppf(confidence_level) else: z_score = stats.norm.ppf(confidence_level) # Calculate the VaR using the GARCH forecasts and the calculated z-score var = -(returns_forecast + z_score * volatility_forecast) return var def asymmetric_quantile_loss(self, actual_return, var_forecast, confidence_level=0.05) -> float: """ Asymmetric quantile loss function Args: actual_return: The actual return var_forecast: The forecasted VaR confidence_level: The confidence level for the VaR Returns: The asymmetric quantile loss """ # Convert the actual return and the forecasted VaR to floats if hasattr(actual_return, 'item'): actual_return = actual_return.item() if hasattr(var_forecast, 'item'): var_forecast = var_forecast.item() # Calculate the actual loss and the forecasted VaR loss actual_loss = -actual_return var_loss = var_forecast # Check if the actual loss exceeds the forecasted VaR exceedance = actual_loss > var_loss # Calculate the asymmetric quantile loss if exceedance: loss = confidence_level * (actual_loss - var_loss) else: loss = (1 - confidence_level) * (var_loss - actual_loss) return loss def regulatory_loss_function(self, actual_return, var_forecast, confidence_level=0.05) -> float: """ Regulatory-style loss function Args: actual_return: The actual return var_forecast: The forecasted VaR confidence_level: The confidence level for the VaR Returns: The regulatory-style loss """ # Convert the actual return and the forecasted VaR to floats if hasattr(actual_return, 'item'): actual_return = actual_return.item() if hasattr(var_forecast, 'item'): var_forecast = var_forecast.item() # Calculate the actual loss and the forecasted VaR loss actual_loss = -actual_return var_loss = var_forecast # Check if the actual loss exceeds the forecasted VaR exceedance = actual_loss > var_loss # Calculate the base loss base_loss = self.asymmetric_quantile_loss(actual_return, var_forecast, confidence_level) # Calculate the regulatory-style loss if exceedance: excess_magnitude = (actual_loss - var_loss) / var_loss if var_loss != 0 else 0 magnitude_penalty = excess_magnitude ** 2 base_loss += magnitude_penalty return base_loss def christoffersen_test(self, exceedances, confidence_level=0.05) -> Dict: """ Christoffersen conditional coverage test Args: exceedances: The exceedances of the VaR confidence_level: The confidence level for the VaR Returns: The Christoffersen conditional coverage test results """ # Convert the exceedances to an array of integers exceedances = np.array(exceedances, dtype=int) n = len(exceedances) n_exceedances = np.sum(exceedances) # If there are no exceedances or all exceedances, return NaN if n_exceedances == 0 or n_exceedances == n: return { 'lr_uc': np.nan, 'lr_cc': np.nan, 'lr_ind': np.nan, 'p_value_uc': np.nan, 'p_value_cc': np.nan, 'p_value_ind': np.nan, 'status': 'Cannot compute (extreme exceedance rate)' } # Calculate the observed and expected exceedance rates p_obs = n_exceedances / n p_exp = confidence_level # Unconditional coverage test (LR test) lr_uc = -2 * (n_exceedances * np.log(p_exp) + (n - n_exceedances) * np.log(1 - p_exp) - n_exceedances * np.log(p_obs) - (n - n_exceedances) * np.log(1 - p_obs)) # Independence test (LR test) n00 = n01 = n10 = n11 = 0 for i in range(n - 1): if exceedances[i] == 0 and exceedances[i + 1] == 0: n00 += 1 elif exceedances[i] == 0 and exceedances[i + 1] == 1: n01 += 1 elif exceedances[i] == 1 and exceedances[i + 1] == 0: n10 += 1 elif exceedances[i] == 1 and exceedances[i + 1] == 1: n11 += 1 # Calculate the total number of exceedances and non-exceedances n_0 = n00 + n01 n_1 = n10 + n11 if n_0 > 0 and n_1 > 0 and n01 > 0 and n10 > 0: p01 = n01 / n_0 p11 = n11 / n_1 lr_ind = -2 * (n01 * np.log(p_obs) + n11 * np.log(p_obs) + n00 * np.log(1 - p_obs) + n10 * np.log(1 - p_obs) - n01 * np.log(p01) - n11 * np.log(p11) - n00 * np.log(1 - p01) - n10 * np.log(1 - p11)) else: lr_ind = np.nan lr_cc = lr_uc + lr_ind if not np.isnan(lr_ind) else np.nan # Calculate the p-values for the unconditional coverage, independence and conditional coverage tests p_value_uc = 1 - stats.chi2.cdf(lr_uc, df=1) if not np.isnan(lr_uc) else np.nan p_value_ind = 1 - stats.chi2.cdf(lr_ind, df=1) if not np.isnan(lr_ind) else np.nan p_value_cc = 1 - stats.chi2.cdf(lr_cc, df=2) if not np.isnan(lr_cc) else np.nan return { 'lr_uc': lr_uc, 'lr_cc': lr_cc, 'lr_ind': lr_ind, 'p_value_uc': p_value_uc, 'p_value_cc': p_value_cc, 'p_value_ind': p_value_ind, 'exceedance_rate': p_obs, 'status': 'computed' } def rolling_var_backtest_single_model(self, returns, model_spec, min_train_window=500, rebalance_freq=REBALANCE_FREQUENCY) -> pd.DataFrame: """ Compute the rolling window VaR backtest for a single model. After each rebalancing, we fit the GARCH model and compute the VaR. Then, the VaR is compared to the actual returns to compute the loss. Args: returns: The returns of the stock model_spec: The model specification min_train_window: The minimum training window rebalance_freq: The rebalancing frequency Returns: The results of the rolling window VaR backtest """ # Get the number of observations and the initial training window n_obs = len(returns) initial_train_size = max(int(n_obs * self.train_ratio), min_train_window) # Initialize the list of results results = [] # Iterate over the training window and rebalance the model for i in range(initial_train_size, n_obs, rebalance_freq): train_start = 0 train_end = i train_returns = returns.iloc[train_start:train_end] try: # Fit the GARCH model garch_result = self.fit_garch_model(train_returns, model_spec, verbose=False) if garch_result is None: continue # Compute the forecasted volatility and expected return forecast = garch_result.forecast(horizon=1) vol_forecast = float(np.sqrt(forecast.variance.iloc[-1, 0])) expected_return = float(train_returns.mean()) # Compute the distribution parameters dist_params = {} if hasattr(garch_result, 'params'): if 'nu' in garch_result.params.index: dist_params['nu'] = garch_result.params['nu'] # Compute the VaR forecast var_forecast = self.calculate_var(expected_return, vol_forecast, self.confidence_level, model_spec['distribution'], dist_params) # Compute the VaR backtest for j in range(min(rebalance_freq, n_obs - i)): if i + j < n_obs: # Compute the actual return actual_return = float(returns.iloc[i + j]) actual_loss = -actual_return # Compute the exceedance exceedance = actual_loss > var_forecast # Compute the asymmetric quantile loss loss_asym = self.asymmetric_quantile_loss(actual_return, var_forecast, self.confidence_level) # Compute the regulatory-style loss loss_reg = self.regulatory_loss_function(actual_return, var_forecast, self.confidence_level) results.append({ 'date': returns.index[i + j], 'var_forecast': var_forecast, 'actual_return': actual_return, 'exceedance': bool(exceedance), 'loss_asymmetric': loss_asym, 'loss_regulatory': loss_reg }) except Exception as e: # Log the error logger.error(f"Error fitting GARCH model: {e}") continue return pd.DataFrame(results) def select_best_models_by_quality(self, all_model_results, max_models=5) -> Dict: """ Select the best models based on quality scores when no models pass coverage tests. Uses the same multi-criteria approach as ensemble weighting but selects top N models. Args: all_model_results (dict): Dictionary of all model results max_models (int): Maximum number of models to select (default: 5) Returns: dict: Dictionary containing the best models """ if len(all_model_results) == 0: return {} # If we have fewer models than max_models, return all if len(all_model_results) <= max_models: logger.info(f"Only {len(all_model_results)} models available - using all") return all_model_results logger.info(f"Selecting best {max_models} models from {len(all_model_results)} available models") # Calculate quality scores using the same approach as ensemble weighting model_names = list(all_model_results.keys()) alpha, beta, gamma, delta = 0.4, 0.3, 0.2, 0.1 # Same weights as ensemble method uc_pvalues = [] cc_pvalues = [] reg_losses = [] exc_rates = [] asym_losses = [] for model_name in model_names: model_data = all_model_results[model_name] uc_pvalues.append(model_data['christoffersen']['p_value_uc']) cc_pvalues.append(model_data['christoffersen']['p_value_cc']) reg_losses.append(model_data['avg_loss_regulatory']) exc_rates.append(abs(model_data['exceedance_rate'] - self.confidence_level)) asym_losses.append(model_data['avg_loss_asymmetric']) uc_pvalues = np.array(uc_pvalues) cc_pvalues = np.array(cc_pvalues) reg_losses = np.array(reg_losses) exc_rates = np.array(exc_rates) asym_losses = np.array(asym_losses) # Handle NaN values uc_pvalues = np.where(np.isnan(uc_pvalues), 1e-6, uc_pvalues) cc_pvalues = np.where(np.isnan(cc_pvalues), 1e-6, cc_pvalues) # Calculate component scores stat_scores = np.sqrt(np.maximum(uc_pvalues, 1e-6) * np.maximum(cc_pvalues, 1e-6)) stat_weights = stat_scores / np.sum(stat_scores) reg_scores = 1.0 / (reg_losses + 1e-10) reg_weights = reg_scores / np.sum(reg_scores) coverage_scores = 1.0 / (exc_rates + 1e-10) coverage_weights = coverage_scores / np.sum(coverage_scores) asym_scores = 1.0 / (asym_losses + 1e-10) asym_weights = asym_scores / np.sum(asym_scores) # Calculate final quality scores final_scores = (alpha * stat_weights + beta * reg_weights + gamma * coverage_weights + delta * asym_weights) # Get indices of top models top_indices = np.argsort(final_scores)[::-1][:max_models] # Sort descending, take top N # Select best models best_models = {} selected_names = [] for idx in top_indices: model_name = model_names[idx] best_models[model_name] = all_model_results[model_name] selected_names.append(model_name) logger.info(f"Selected models: {selected_names}") logger.info(f"Quality scores: {[f'{final_scores[idx]:.4f}' for idx in top_indices]}") return best_models def calculate_ensemble_weights_multicriteria(self, acceptable_models, alpha=0.4, beta=0.3, gamma=0.2, delta=0.1) -> Tuple[Dict, Dict]: """ Calculate the ensemble weights for a given set of models using a multi-criteria approach. Args: acceptable_models (dict): Dictionary of acceptable models alpha (float): Weight for the statistical component beta (float): Weight for the regulatory component gamma (float): Weight for the coverage component delta (float): Weight for the asymmetric component Returns: Tuple[Dict, Dict]: Tuple containing the ensemble weights and the weight breakdown """ # Initialize the lists to store the model data model_names = list(acceptable_models.keys()) uc_pvalues = [] cc_pvalues = [] reg_losses = [] exc_rates = [] asym_losses = [] # For each model, append the Christoffersen p-values, regulatory losses, exceedance rates, and asymmetric losses for model_name in model_names: model_data = acceptable_models[model_name] uc_pvalues.append(model_data['christoffersen']['p_value_uc']) cc_pvalues.append(model_data['christoffersen']['p_value_cc']) reg_losses.append(model_data['avg_loss_regulatory']) exc_rates.append(abs(model_data['exceedance_rate'] - self.confidence_level)) asym_losses.append(model_data['avg_loss_asymmetric']) uc_pvalues = np.array(uc_pvalues) cc_pvalues = np.array(cc_pvalues) reg_losses = np.array(reg_losses) exc_rates = np.array(exc_rates) asym_losses = np.array(asym_losses) uc_pvalues = np.where(np.isnan(uc_pvalues), 1e-6, uc_pvalues) cc_pvalues = np.where(np.isnan(cc_pvalues), 1e-6, cc_pvalues) # Calculate component weights stat_scores = np.sqrt(np.maximum(uc_pvalues, 1e-6) * np.maximum(cc_pvalues, 1e-6)) stat_weights = stat_scores / np.sum(stat_scores) reg_scores = 1.0 / (reg_losses + 1e-10) reg_weights = reg_scores / np.sum(reg_scores) coverage_scores = 1.0 / (exc_rates + 1e-10) coverage_weights = coverage_scores / np.sum(coverage_scores) asym_scores = 1.0 / (asym_losses + 1e-10) asym_weights = asym_scores / np.sum(asym_scores) # Combine weights final_weights = (alpha * stat_weights + beta * reg_weights + gamma * coverage_weights + delta * asym_weights) final_weights = final_weights / np.sum(final_weights) weight_breakdown = { 'statistical': dict(zip(model_names, stat_weights)), 'regulatory': dict(zip(model_names, reg_weights)), 'coverage': dict(zip(model_names, coverage_weights)), 'asymmetric': dict(zip(model_names, asym_weights)), 'final': dict(zip(model_names, final_weights)) } return dict(zip(model_names, final_weights)), weight_breakdown def fit_ensemble_univariate_model(self, returns, ticker) -> Dict: """ Fit ensemble of GARCH models for a single stock Args: returns (pd.Series): The returns of the stock ticker (str): The ticker symbol of the stock Returns: Dict: Dictionary containing the ensemble information """ logger.info(f"Fitting ensemble models for {ticker}...") # Create all model specifications model_specs = self.create_garch_model_specifications() all_model_results = {} # Test each model specification with faster rebalancing for portfolio analysis for i, model_spec in enumerate(model_specs): try: model_results = self.rolling_var_backtest_single_model( returns, model_spec, min_train_window=500, rebalance_freq=REBALANCE_FREQUENCY ) if len(model_results) > 0: exceedances = model_results['exceedance'].values total_obs = len(model_results) total_exceedances = np.sum(exceedances) exceedance_rate = total_exceedances / total_obs avg_loss_asym = model_results['loss_asymmetric'].mean() avg_loss_reg = model_results['loss_regulatory'].mean() christoffersen_results = self.christoffersen_test(exceedances, self.confidence_level) all_model_results[model_spec['name']] = { 'model_spec': model_spec, 'results_df': model_results, 'total_obs': total_obs, 'total_exceedances': total_exceedances, 'exceedance_rate': exceedance_rate, 'avg_loss_asymmetric': avg_loss_asym, 'avg_loss_regulatory': avg_loss_reg, 'christoffersen': christoffersen_results } # Memory cleanup every 10 models if i % 10 == 0: gc.collect() except Exception as e: logger.error(f"Error fitting model {model_spec['name']}: {e}") continue logger.info(f"Successfully tested {len(all_model_results)} models for {ticker}") # Filter acceptable models acceptable_models = {} for model_name, results in all_model_results.items(): uc_pvalue = results['christoffersen']['p_value_uc'] cc_pvalue = results['christoffersen']['p_value_cc'] # Only add the model if it passes the coverage tests with a significance level of 0.05 if (not np.isnan(uc_pvalue) and not np.isnan(cc_pvalue) and uc_pvalue > 0.05 and cc_pvalue > 0.05): acceptable_models[model_name] = results if len(acceptable_models) == 0: logger.warning(f"No models pass coverage tests for {ticker} - selecting best 5 models based on quality scores") # Calculate quality scores for all models to select the best ones best_models = self.select_best_models_by_quality(all_model_results, max_models=5) acceptable_models = best_models logger.info(f"Using {len(acceptable_models)} models for {ticker} ensemble") # Calculate ensemble weights weight_dict, weight_breakdown = self.calculate_ensemble_weights_multicriteria(acceptable_models) # Fit final models on full data final_fitted_models = {} for model_name in acceptable_models.keys(): model_spec = acceptable_models[model_name]['model_spec'] try: fitted_model = self.fit_garch_model(returns, model_spec, verbose=False) if fitted_model is not None: final_fitted_models[model_name] = fitted_model except Exception as e: continue ensemble_info = { 'acceptable_models': acceptable_models, 'fitted_models': final_fitted_models, 'weights': weight_dict, 'weight_breakdown': weight_breakdown, 'ticker': ticker, 'validation_summary': { 'total_models_tested': len(model_specs), 'successful_models': len(all_model_results), 'acceptable_models': len(acceptable_models), 'final_fitted_models': len(final_fitted_models) } } return ensemble_info def fit_all_ensemble_models(self) -> None: """ Fit ensemble models for all portfolio stocks """ logger.info("Fitting ensemble models for all portfolio stocks...") for ticker in self.tickers: logger.info(f"Processing {ticker}...") ensemble_info = self.fit_ensemble_univariate_model(self.returns[ticker], ticker) self.ensemble_models[ticker] = ensemble_info # Get conditional volatility conditional_vol = self.get_ensemble_conditional_volatility(ensemble_info, self.returns[ticker]) self.conditional_volatilities[ticker] = conditional_vol # Get standardized residuals for correlation estimation self.standardized_residuals = pd.DataFrame(index=self.returns.index) for ticker in self.tickers: std_resid = self.get_ensemble_standardized_residuals(self.ensemble_models[ticker], self.returns[ticker]) self.standardized_residuals[ticker] = std_resid self.standardized_residuals = self.standardized_residuals.dropna() logger.info("Completed ensemble model fitting for all stocks") def get_ensemble_conditional_volatility(self, ensemble_info, returns) -> pd.Series: """ Get weighted conditional volatility from ensemble models Args: ensemble_info (dict): Dictionary containing the ensemble information returns (pd.Series): The returns of the stock Returns: pd.Series: The weighted conditional volatility """ # Initialize the weighted volatility weighted_volatility = pd.Series(0.0, index=returns.index) total_weight = 0.0 # For each model, append the conditional volatility for model_name, fitted_model in ensemble_info['fitted_models'].items(): try: weight = ensemble_info['weights'][model_name] conditional_vol = fitted_model.conditional_volatility aligned_vol = conditional_vol.reindex(returns.index, fill_value=0) weighted_volatility += weight * aligned_vol total_weight += weight except Exception as e: continue if total_weight > 0: weighted_volatility = weighted_volatility / total_weight return weighted_volatility def get_ensemble_standardized_residuals(self, ensemble_info, returns) -> pd.Series: """ Get weighted standardized residuals from ensemble models Args: ensemble_info (dict): Dictionary containing the ensemble information returns (pd.Series): The returns of the stock Returns: pd.Series: The weighted standardized residuals """ # Initialize the weighted residuals weighted_residuals = pd.Series(0.0, index=returns.index) total_weight = 0.0 # For each model, append the standardized residuals for model_name, fitted_model in ensemble_info['fitted_models'].items(): try: weight = ensemble_info['weights'][model_name] residuals = fitted_model.resid conditional_vol = fitted_model.conditional_volatility std_residuals = residuals / conditional_vol aligned_residuals = std_residuals.reindex(returns.index, fill_value=0) weighted_residuals += weight * aligned_residuals total_weight += weight except Exception as e: continue if total_weight > 0: weighted_residuals = weighted_residuals / total_weight return weighted_residuals def estimate_correlation_matrix(self) -> None: """ Estimate constant conditional correlation matrix using the standardized residuals. """ logger.info("Estimating correlation matrix...") if self.standardized_residuals is None: raise ValueError("Must fit ensemble models first") self.correlation_matrix = self.standardized_residuals.corr().values logger.info("Correlation Matrix:") corr_df = pd.DataFrame(self.correlation_matrix, columns=self.tickers, index=self.tickers) logger.info(f"\n{corr_df.round(4)}") def forecast_ensemble_volatility(self, ensemble_info, horizon=1) -> float: """ Generate volatility forecast from ensemble Args: ensemble_info (dict): Dictionary containing the ensemble information returns (pd.Series): The returns of the stock horizon (int): The horizon for the forecast Returns: float: The weighted forecast """ # Initialize the weighted forecast weighted_forecast = 0.0 total_weight = 0.0 # For each model, append the forecasted volatility for model_name, fitted_model in ensemble_info['fitted_models'].items(): try: weight = ensemble_info['weights'][model_name] forecast = fitted_model.forecast(horizon=horizon) vol_forecast = np.sqrt(forecast.variance.iloc[-1, 0]) weighted_forecast += weight * vol_forecast total_weight += weight except Exception as e: continue if total_weight > 0: weighted_forecast = weighted_forecast / total_weight return weighted_forecast def forecast_covariance_matrix(self, horizon=1) -> Tuple[np.ndarray, List[float]]: """ Forecast covariance matrix using ensemble models Args: horizon (int): The horizon for the forecast Returns: Tuple[np.ndarray, List[float]]: The forecasted covariance matrix and the volatilities """ logger.info(f"Forecasting covariance matrix for {horizon} period(s) ahead...") # Initialize the volatilities vol_forecasts = [] # For each ticker, forecast the volatility for ticker in self.tickers: vol_forecast = self.forecast_ensemble_volatility(self.ensemble_models[ticker], horizon) vol_forecasts.append(vol_forecast) logger.info(f"{ticker} forecasted volatility: {vol_forecast:.4f}") # Construct covariance matrix vol_matrix = np.diag(vol_forecasts) forecast_cov = vol_matrix @ self.correlation_matrix @ vol_matrix return forecast_cov, vol_forecasts def optimize_portfolio(self) -> None: """ Portfolio optimization using CCC ensemble model """ logger.info("Optimizing portfolio...") # Get forecasted covariance matrix forecast_cov, vol_forecasts = self.forecast_covariance_matrix() logger.info("Forecasted volatilities:") for i, ticker in enumerate(self.tickers): logger.info(f"{ticker}: {vol_forecasts[i]:.4f}") print("\n") # Minimum variance portfolio n_assets = len(self.tickers) def portfolio_variance(weights, cov_matrix): return weights.T @ cov_matrix @ weights # Constraints and bounds constraints = {'type': 'eq', 'fun': lambda w: np.sum(w) - 1} bounds = [(0, 1) for _ in range(n_assets)] initial_weights = np.array([1/n_assets] * n_assets) # Optimize result = minimize(portfolio_variance, initial_weights, args=(forecast_cov,), method='SLSQP', bounds=bounds, constraints=constraints, options={'maxiter': 1000, 'ftol': 1e-10}) if result.success: optimal_weights = result.x min_variance = result.fun logger.info("Portfolio optimization successful!") logger.info("Optimal weights:") for i, ticker in enumerate(self.tickers): logger.info(f" {ticker}: {optimal_weights[i]:.4f}") logger.info(f"Portfolio Variance: {min_variance:.6f}") logger.info(f"Portfolio Volatility: {np.sqrt(min_variance):.4f}") # Store optimization results self.optimization_results = { 'optimal_weights': optimal_weights, 'portfolio_variance': min_variance, 'portfolio_volatility': np.sqrt(min_variance), 'forecasted_volatilities': vol_forecasts, 'forecasted_covariance': forecast_cov, 'current_weights': self.weights, 'optimization_success': True } else: logger.error("Portfolio optimization failed!") logger.error(f"Reason: {result.message}") # Use current weights as fallback optimal_weights = np.array(self.weights) min_variance = portfolio_variance(optimal_weights, forecast_cov) self.optimization_results = { 'optimal_weights': optimal_weights, 'portfolio_variance': min_variance, 'portfolio_volatility': np.sqrt(min_variance), 'forecasted_volatilities': vol_forecasts, 'forecasted_covariance': forecast_cov, 'current_weights': self.weights, 'optimization_success': False } return self.optimization_results def generate_results(self) -> Dict: """ Generate cmodeling results for the portfolio Returns: Dict: Dictionary containing the results """ logger.info("Generating results...") # Generate forecasts forecast_results = self.volatility_forecast() # Create plots plot_files = self.create_plots() # Portfolio summary portfolio_summary = { 'total_value': self.total_value, 'num_stocks': len(self.tickers), 'tickers': self.tickers, 'amounts': self.amounts, 'weights': self.weights, 'data_period': { 'start_date': self.start_date, 'end_date': self.end_date, 'total_observations': len(self.returns) } } # Ensemble model summary ensemble_summary = {} for ticker in self.tickers: ensemble_info = self.ensemble_models[ticker] ensemble_summary[ticker] = { 'total_models_tested': ensemble_info['validation_summary']['total_models_tested'], 'successful_models': ensemble_info['validation_summary']['successful_models'], 'acceptable_models': ensemble_info['validation_summary']['acceptable_models'], 'final_fitted_models': ensemble_info['validation_summary']['final_fitted_models'], 'top_3_models': dict(sorted(ensemble_info['weights'].items(), key=lambda x: x[1], reverse=True)[:3]), 'avg_conditional_volatility': self.conditional_volatilities[ticker].mean(), 'model_performance_stats': { 'avg_exceedance_rate': np.mean([ model_data['exceedance_rate'] for model_data in ensemble_info['acceptable_models'].values() ]), 'avg_regulatory_loss': np.mean([ model_data['avg_loss_regulatory'] for model_data in ensemble_info['acceptable_models'].values() ]), 'avg_uc_pvalue': np.mean([ model_data['christoffersen']['p_value_uc'] for model_data in ensemble_info['acceptable_models'].values() if not np.isnan(model_data['christoffersen']['p_value_uc']) ]) } } # Risk metrics # Handle correlation calculations for single vs multiple stocks if len(self.tickers) == 1: # Single stock - no correlations to calculate avg_correlation = None correlation_analysis = { 'min_correlation': None, 'max_correlation': None, 'correlation_std': None } else: # Multiple stocks - calculate correlations upper_triangle = self.correlation_matrix[np.triu_indices_from(self.correlation_matrix, k=1)] avg_correlation = np.mean(upper_triangle) correlation_analysis = { 'min_correlation': float(np.min(upper_triangle)), 'max_correlation': float(np.max(upper_triangle)), 'correlation_std': float(np.std(upper_triangle)) } risk_metrics = { 'correlation_matrix': self.correlation_matrix.tolist(), 'average_correlation': avg_correlation, 'portfolio_volatility_forecast': self.optimization_results['portfolio_volatility'], 'individual_volatility_forecasts': dict(zip(self.tickers, self.optimization_results['forecasted_volatilities'])), 'comprehensive_forecasts': forecast_results, 'correlation_analysis': correlation_analysis } # Optimization results optimization_summary = { 'optimization_successful': self.optimization_results['optimization_success'], 'optimal_weights': dict(zip(self.tickers, self.optimization_results['optimal_weights'])), 'current_weights': dict(zip(self.tickers, self.optimization_results['current_weights'])), 'portfolio_variance': self.optimization_results['portfolio_variance'], 'portfolio_volatility': self.optimization_results['portfolio_volatility'], 'weight_changes': { ticker: self.optimization_results['optimal_weights'][i] - self.optimization_results['current_weights'][i] for i, ticker in enumerate(self.tickers) }, 'rebalancing_analysis': { 'significant_changes': [ (ticker, change) for ticker, change in [(ticker, self.optimization_results['optimal_weights'][i] - self.optimization_results['current_weights'][i]) for i, ticker in enumerate(self.tickers)] if abs(change) > 0.05 ], 'total_turnover': sum([ abs(self.optimization_results['optimal_weights'][i] - self.optimization_results['current_weights'][i]) for i in range(len(self.tickers)) ]) / 2 } } # Compile all results self.modeling_results = { 'portfolio_summary': portfolio_summary, 'ensemble_summary': ensemble_summary, 'risk_metrics': risk_metrics, 'optimization_results': optimization_summary, 'plot_files': plot_files, 'model_metadata': { 'confidence_level': self.confidence_level, 'train_ratio': self.train_ratio, 'modeling_approach': 'CCC-GARCH with Ensemble Univariate Models', 'rebalance_frequency': REBALANCE_FREQUENCY, 'generation_timestamp': datetime.now().isoformat(), 'technical_details': { 'garch_specifications': len(self.create_garch_model_specifications()), 'distributions_used': ['normal', 't', 'ged'], 'model_types': ['GARCH', 'EGARCH', 'TGARCH'], 'validation_method': 'Rolling window VaR backtesting with Christoffersen tests', 'ensemble_weighting': 'Multi-criteria: statistical validity (40%), regulatory loss (30%), coverage accuracy (20%), asymmetric loss (10%)', 'correlation_estimation': 'Constant Conditional Correlation from standardized residuals', 'optimization_method': 'Minimum variance with full investment constraint' } } } logger.info("Results generated successfully") return self.modeling_results def run_complete_analysis(self) -> Dict: """ Run the complete portfolio modeling pipeline Returns: Dict: Dictionary containing the results """ logger.info("Starting complete portfolio analysis...") try: # Step 1: Download data if not self.download_portfolio_data(): raise Exception("Failed to download portfolio data") # Step 2: Fit ensemble models self.fit_all_ensemble_models() # Step 3: Estimate correlation matrix self.estimate_correlation_matrix() # Step 4: Optimize portfolio self.optimize_portfolio() # Step 5: Generate results results = self.generate_results() logger.info("Complete portfolio analysis finished successfully!") return results except Exception as e: logger.error(f"Error in portfolio analysis: {str(e)}") raise def save_results_to_file(self, filename=None) -> Optional[str]: """ Save modeling results to JSON file Args: filename (str): The name of the file to save the results to Returns: Optional[str]: The name of the file the results were saved to """ if filename is None: timestamp = datetime.now().strftime('%Y%m%d_%H%M%S') filename = f"portfolio_analysis_{timestamp}.json" try: with open(filename, 'w') as f: json.dump(self.modeling_results, f, indent=2, default=str) logger.info(f"Results saved to {filename}") return filename except Exception as e: logger.error(f"Error saving results: {str(e)}") return None def volatility_forecast(self, horizons=[1, 5, 10, 22], confidence_levels=[0.05, 0.10, 0.25]) -> Dict: """ Generate volatility forecasts with multiple horizons and confidence intervals Args: horizons: List of forecast horizons (in days) confidence_levels: List of confidence levels for intervals Returns: Dictionary with forecast results """ logger.info("Generating volatility forecasts...") forecast_results = {} for ticker in self.tickers: logger.info(f"Forecasting volatility for {ticker}...") ticker_forecasts = { 'horizons': horizons, 'confidence_levels': confidence_levels, 'point_forecasts': [], 'confidence_intervals': {}, 'var_forecasts': {} } ensemble_info = self.ensemble_models[ticker] # Initialize confidence intervals storage for conf_level in confidence_levels: ticker_forecasts['confidence_intervals'][conf_level] = { 'lower': [], 'upper': [] } ticker_forecasts['var_forecasts'][conf_level] = [] # Generate forecasts for each horizon for horizon in horizons: # Collect forecasts from all models in ensemble model_forecasts = [] var_forecasts_by_conf = {conf: [] for conf in confidence_levels} for model_name, fitted_model in ensemble_info['fitted_models'].items(): try: weight = ensemble_info['weights'][model_name] # Get model forecast forecast = fitted_model.forecast(horizon=horizon) vol_forecast = np.sqrt(forecast.variance.iloc[-1, 0]) # Add to collection (weighted) for _ in range(int(weight * 1000)): # Weight by adding multiple copies model_forecasts.append(vol_forecast) # Calculate VaR for each confidence level model_spec = ensemble_info['acceptable_models'][model_name]['model_spec'] expected_return = self.returns[ticker].mean() for conf_level in confidence_levels: var_forecast = self.calculate_var( expected_return, vol_forecast, conf_level, model_spec['distribution'] ) for _ in range(int(weight * 1000)): var_forecasts_by_conf[conf_level].append(var_forecast) except Exception as e: continue # Calculate ensemble statistics if model_forecasts: point_forecast = np.mean(model_forecasts) ticker_forecasts['point_forecasts'].append(point_forecast) # Calculate confidence intervals for volatility for conf_level in confidence_levels: lower_bound = np.percentile(model_forecasts, (conf_level/2) * 100) upper_bound = np.percentile(model_forecasts, (1 - conf_level/2) * 100) ticker_forecasts['confidence_intervals'][conf_level]['lower'].append(lower_bound) ticker_forecasts['confidence_intervals'][conf_level]['upper'].append(upper_bound) # VaR forecast var_point = np.mean(var_forecasts_by_conf[conf_level]) ticker_forecasts['var_forecasts'][conf_level].append(var_point) else: # Fallback to single forecast vol_forecast = self.forecast_ensemble_volatility(ensemble_info, horizon) ticker_forecasts['point_forecasts'].append(vol_forecast) for conf_level in confidence_levels: # Simple confidence intervals (±20% of point forecast) margin = vol_forecast * 0.2 ticker_forecasts['confidence_intervals'][conf_level]['lower'].append(vol_forecast - margin) ticker_forecasts['confidence_intervals'][conf_level]['upper'].append(vol_forecast + margin) # VaR calculation expected_return = self.returns[ticker].mean() var_forecast = self.calculate_var(expected_return, vol_forecast, conf_level) ticker_forecasts['var_forecasts'][conf_level].append(var_forecast) forecast_results[ticker] = ticker_forecasts # Portfolio-level forecasts portfolio_forecasts = { 'horizons': horizons, 'portfolio_volatility': [], 'portfolio_var': {} } for conf_level in confidence_levels: portfolio_forecasts['portfolio_var'][conf_level] = [] for horizon in horizons: # Get individual volatility forecasts for this horizon individual_vols = [forecast_results[ticker]['point_forecasts'][horizons.index(horizon)] for ticker in self.tickers] # Calculate portfolio volatility vol_matrix = np.diag(individual_vols) portfolio_cov = vol_matrix @ self.correlation_matrix @ vol_matrix portfolio_vol = np.sqrt(np.array(self.weights) @ portfolio_cov @ np.array(self.weights)) portfolio_forecasts['portfolio_volatility'].append(portfolio_vol) # Portfolio VaR portfolio_return = np.sum([w * self.returns[ticker].mean() for w, ticker in zip(self.weights, self.tickers)]) for conf_level in confidence_levels: portfolio_var = self.calculate_var(portfolio_return, portfolio_vol, conf_level) portfolio_forecasts['portfolio_var'][conf_level].append(portfolio_var) results = { 'individual_forecasts': forecast_results, 'portfolio_forecasts': portfolio_forecasts, 'metadata': { 'horizons': horizons, 'confidence_levels': confidence_levels, 'forecast_date': datetime.now().isoformat(), 'base_currency': 'USD' } } logger.info("Volatility forecasting completed") return results def create_plots(self, save_directory='plots') -> Dict: """ Create visualization plots Args: save_directory: Directory to save plots Returns: Dictionary with plot file paths """ logger.info("Creating plots...") # Create plots directory if not os.path.exists(save_directory): os.makedirs(save_directory) plot_files = {} # Plot 1: Cumulative Returns Evolution plt.figure(figsize=(12, 8)) cumulative_returns = self.returns.cumsum().apply(np.exp) for ticker in self.tickers: plt.plot(cumulative_returns.index, cumulative_returns[ticker], label=ticker, linewidth=2) plt.title('Portfolio Cumulative Returns Evolution', fontsize=16, fontweight='bold') plt.xlabel('Date', fontsize=12) plt.ylabel('Cumulative Return', fontsize=12) plt.legend(fontsize=12) plt.grid(True, alpha=0.3) plt.tight_layout() plot_path = os.path.join(save_directory, '01_cumulative_returns.png') plt.savefig(plot_path, dpi=300, bbox_inches='tight') plt.close() plot_files['cumulative_returns'] = plot_path # Plot 2: Daily Returns plt.figure(figsize=(14, 8)) for ticker in self.tickers: plt.plot(self.returns.index, self.returns[ticker], label=ticker, alpha=0.7, linewidth=0.8) plt.title('Daily Returns Time Series', fontsize=16, fontweight='bold') plt.xlabel('Date', fontsize=12) plt.ylabel('Daily Return', fontsize=12) plt.legend(fontsize=12) plt.grid(True, alpha=0.3) plt.tight_layout() plot_path = os.path.join(save_directory, '02_daily_returns.png') plt.savefig(plot_path, dpi=300, bbox_inches='tight') plt.close() plot_files['daily_returns'] = plot_path # Plot 3: Ensemble Conditional Volatilities plt.figure(figsize=(14, 8)) for ticker in self.tickers: vol_series = self.conditional_volatilities[ticker] plt.plot(vol_series.index, vol_series, label=f'{ticker} Ensemble Volatility', linewidth=2) plt.title('Ensemble Conditional Volatilities', fontsize=16, fontweight='bold') plt.xlabel('Date', fontsize=12) plt.ylabel('Conditional Volatility', fontsize=12) plt.legend(fontsize=12) plt.grid(True, alpha=0.3) plt.tight_layout() plot_path = os.path.join(save_directory, '03_ensemble_volatilities.png') plt.savefig(plot_path, dpi=300, bbox_inches='tight') plt.close() plot_files['ensemble_volatilities'] = plot_path # Plot 4: Rolling vs Model Volatilities Comparison plt.figure(figsize=(14, 8)) window = 60 for ticker in self.tickers: rolling_vol = self.returns[ticker].rolling(window=window).std() model_vol = self.conditional_volatilities[ticker] plt.plot(rolling_vol.index, rolling_vol, label=f'{ticker} Rolling Vol ({window}d)', alpha=0.7, linewidth=1.5, linestyle='--') plt.plot(model_vol.index, model_vol.rolling(5).mean(), label=f'{ticker} Ensemble Vol', linewidth=2) plt.title('Rolling vs Ensemble Volatilities Comparison', fontsize=16, fontweight='bold') plt.xlabel('Date', fontsize=12) plt.ylabel('Volatility', fontsize=12) plt.legend(fontsize=10) plt.grid(True, alpha=0.3) plt.tight_layout() plot_path = os.path.join(save_directory, '04_volatility_comparison.png') plt.savefig(plot_path, dpi=300, bbox_inches='tight') plt.close() plot_files['volatility_comparison'] = plot_path # Plot 5: Model Weights for each asset for ticker in self.tickers: plt.figure(figsize=(12, 8)) weights = self.ensemble_models[ticker]['weights'] # Get all models and their weights, sort by weight sorted_weights = sorted(weights.items(), key=lambda x: x[1], reverse=True) models, weight_values = zip(*sorted_weights) # Clean model names for display clean_models = [m.replace('_', '\n').replace('(', '\n(') for m in models] colors = plt.cm.viridis(np.linspace(0, 1, len(models))) bars = plt.barh(range(len(models)), weight_values, color=colors) plt.yticks(range(len(models)), clean_models, fontsize=8) plt.xlabel('Weight', fontsize=12) plt.title(f'{ticker} - Ensemble Model Weights', fontsize=16, fontweight='bold') plt.grid(True, alpha=0.3, axis='x') # Add weight values on bars for i, (bar, weight) in enumerate(zip(bars, weight_values)): plt.text(weight + 0.001, i, f'{weight:.3f}', va='center', fontsize=8) plt.tight_layout() plot_path = os.path.join(save_directory, f'05_model_weights_{ticker}.png') plt.savefig(plot_path, dpi=300, bbox_inches='tight') plt.close() plot_files[f'model_weights_{ticker}'] = plot_path # Plot 6: Correlation Matrix Heatmap plt.figure(figsize=(10, 8)) im = plt.imshow(self.correlation_matrix, cmap='RdBu_r', vmin=-1, vmax=1) plt.colorbar(im, label='Correlation Coefficient') # Add text annotations for i in range(len(self.tickers)): for j in range(len(self.tickers)): plt.text(j, i, f'{self.correlation_matrix[i, j]:.3f}', ha='center', va='center', fontsize=14, fontweight='bold') plt.xticks(range(len(self.tickers)), self.tickers, fontsize=12) plt.yticks(range(len(self.tickers)), self.tickers, fontsize=12) plt.title('CCC Correlation Matrix', fontsize=16, fontweight='bold') plt.tight_layout() plot_path = os.path.join(save_directory, '06_correlation_matrix.png') plt.savefig(plot_path, dpi=300, bbox_inches='tight') plt.close() plot_files['correlation_matrix'] = plot_path # Plot 7: Volatility Forecasts forecast_results = self.volatility_forecast() for ticker in self.tickers: plt.figure(figsize=(12, 8)) ticker_data = forecast_results['individual_forecasts'][ticker] horizons = ticker_data['horizons'] point_forecasts = ticker_data['point_forecasts'] # Plot point forecasts plt.plot(horizons, point_forecasts, 'o-', linewidth=2, markersize=8, label='Point Forecast', color='blue') # Plot confidence intervals colors = ['red', 'orange', 'green'] for i, conf_level in enumerate(ticker_data['confidence_levels']): lower = ticker_data['confidence_intervals'][conf_level]['lower'] upper = ticker_data['confidence_intervals'][conf_level]['upper'] plt.fill_between(horizons, lower, upper, alpha=0.3, color=colors[i], label=f'{(1-conf_level)*100:.0f}% Confidence Interval') plt.xlabel('Forecast Horizon (Days)', fontsize=12) plt.ylabel('Volatility', fontsize=12) plt.title(f'{ticker} - Volatility Forecast with Confidence Intervals', fontsize=16, fontweight='bold') plt.legend(fontsize=10) plt.grid(True, alpha=0.3) plt.tight_layout() plot_path = os.path.join(save_directory, f'07_volatility_forecast_{ticker}.png') plt.savefig(plot_path, dpi=300, bbox_inches='tight') plt.close() plot_files[f'volatility_forecast_{ticker}'] = plot_path # Plot 8: Portfolio VaR Forecasts plt.figure(figsize=(12, 8)) portfolio_data = forecast_results['portfolio_forecasts'] horizons = portfolio_data['horizons'] colors = ['red', 'orange', 'green'] for i, conf_level in enumerate([0.05, 0.10, 0.25]): if conf_level in portfolio_data['portfolio_var']: var_forecasts = portfolio_data['portfolio_var'][conf_level] plt.plot(horizons, var_forecasts, 'o-', linewidth=2, markersize=8, color=colors[i], label=f'VaR {conf_level*100:.0f}%') plt.xlabel('Forecast Horizon (Days)', fontsize=12) plt.ylabel('Value at Risk', fontsize=12) plt.title('Portfolio VaR Forecasts', fontsize=16, fontweight='bold') plt.legend(fontsize=12) plt.grid(True, alpha=0.3) plt.tight_layout() plot_path = os.path.join(save_directory, '08_portfolio_var_forecasts.png') plt.savefig(plot_path, dpi=300, bbox_inches='tight') plt.close() plot_files['portfolio_var_forecasts'] = plot_path # Plot 9: Model Performance Summary for ticker in self.tickers: plt.figure(figsize=(12, 8)) ensemble_info = self.ensemble_models[ticker] exc_rates = [] reg_losses = [] uc_pvalues = [] model_names = [] for model_name, model_data in ensemble_info['acceptable_models'].items(): exc_rates.append(model_data['exceedance_rate']) reg_losses.append(model_data['avg_loss_regulatory']) uc_pvalues.append(model_data['christoffersen']['p_value_uc']) model_names.append(model_name.split('_')[0]) # Create scatter plot scatter = plt.scatter(exc_rates, reg_losses, c=uc_pvalues, s=100, alpha=0.7, cmap='RdYlGn', vmin=0, vmax=0.1) # Add reference lines plt.axvline(x=0.05, color='red', linestyle='--', alpha=0.5, label='Expected 5% VaR') plt.xlabel('Exceedance Rate', fontsize=12) plt.ylabel('Regulatory Loss', fontsize=12) plt.title(f'{ticker} - Model Performance (Color = UC p-value)', fontsize=16, fontweight='bold') plt.colorbar(scatter, label='UC p-value') plt.legend(fontsize=10) plt.grid(True, alpha=0.3) plt.tight_layout() plot_path = os.path.join(save_directory, f'09_model_performance_{ticker}.png') plt.savefig(plot_path, dpi=300, bbox_inches='tight') plt.close() plot_files[f'model_performance_{ticker}'] = plot_path # Plot 10: Ensemble Validation Summary plt.figure(figsize=(14, 8)) # Collect summary data total_tested = [] successful = [] acceptable = [] final_fitted = [] for ticker in self.tickers: summary = self.ensemble_models[ticker]['validation_summary'] total_tested.append(summary['total_models_tested']) successful.append(summary['successful_models']) acceptable.append(summary['acceptable_models']) final_fitted.append(summary['final_fitted_models']) x = np.arange(len(self.tickers)) width = 0.2 plt.bar(x - 1.5*width, total_tested, width, label='Total Tested', alpha=0.8) plt.bar(x - 0.5*width, successful, width, label='Successful', alpha=0.8) plt.bar(x + 0.5*width, acceptable, width, label='Acceptable (Coverage)', alpha=0.8) plt.bar(x + 1.5*width, final_fitted, width, label='Final Fitted', alpha=0.8) plt.xticks(x, self.tickers, fontsize=12) plt.ylabel('Number of Models', fontsize=12) plt.title('Ensemble Model Validation Summary by Asset', fontsize=16, fontweight='bold') plt.legend(fontsize=10) plt.grid(True, alpha=0.3, axis='y') plt.tight_layout() plot_path = os.path.join(save_directory, '10_validation_summary.png') plt.savefig(plot_path, dpi=300, bbox_inches='tight') plt.close() plot_files['validation_summary'] = plot_path logger.info(f"Created {len(plot_files)} plots in {save_directory}") return plot_files def run_portfolio_analysis(portfolio_data, start_date='2014-01-01', end_date=None, confidence_level=0.05, train_ratio=0.8, save_results=True) -> Dict: """ Standalone function to run complete portfolio analysis Args: portfolio_data: List of (ticker, amount) tuples start_date: Start date for historical data end_date: End date for historical data confidence_level: VaR confidence level train_ratio: Training/testing split ratio save_results: Whether to save results to file Returns: Dictionary with modeling results """ # Create model instance model = PortfolioCCCEnsembleModel( portfolio_data=portfolio_data, start_date=start_date, end_date=end_date, confidence_level=confidence_level, train_ratio=train_ratio ) # Run complete analysis results = model.run_complete_analysis() # Save results if requested if save_results: model.save_results_to_file() return results class PortfolioDCCEnsembleModel: """ This class implements Dynamic Conditional Correlation (DCC) GARCH modeling for portfolios. Unlike CCC, DCC allows correlations to change over time, providing more flexible modeling of time-varying dependencies between assets. """ def __init__(self, portfolio_data, start_date='2014-01-01', end_date=None, confidence_level=0.05, train_ratio=0.8) -> None: """ Initialize the DCC portfolio modeling system Args: portfolio_data: List of (ticker, amount) tuples start_date: Start date for historical data end_date: End date for historical data (None = today) confidence_level: VaR confidence level train_ratio: Training/testing split ratio """ # Initialize the portfolio data self.portfolio_data = portfolio_data self.start_date = start_date self.end_date = end_date or datetime.now().strftime('%Y-%m-%d') self.confidence_level = confidence_level self.train_ratio = train_ratio # Portfolio composition self.tickers = [ticker for ticker, amount in portfolio_data] self.amounts = [float(amount) for ticker, amount in portfolio_data] self.total_value = sum(self.amounts) self.weights = [amount / self.total_value for amount in self.amounts] # Data storage self.returns = None self.prices = None self.ensemble_models = {} self.conditional_volatilities = {} self.standardized_residuals = None # DCC-specific parameters self.dcc_params = None # alpha, beta parameters self.dynamic_correlations = None # Time-varying correlation matrices self.unconditional_correlation = None # Long-run correlation matrix self.q_matrices = None # DCC Q matrices # Results storage self.modeling_results = {} self.optimization_results = {} logger.info(f"Initialized DCC portfolio model with {len(self.tickers)} stocks") logger.info(f"Portfolio value: ${self.total_value:,.2f}") logger.info(f"Data period: {self.start_date} to {self.end_date}") # Inherit basic methods from CCC model def download_portfolio_data(self) -> bool: """ Download historical data for all portfolio stocks Returns: bool: True if data was downloaded successfully, False otherwise """ logger.info("Downloading portfolio data...") try: # Download data for all tickers data = yf.download(self.tickers, start=self.start_date, end=self.end_date) if len(self.tickers) == 1: # Handle single stock case - data['Close'] is already a Series if isinstance(data['Close'], pd.Series): prices = data['Close'].to_frame() prices.columns = self.tickers else: # If it's already a DataFrame prices = data['Close'] if prices.shape[1] == 1: prices.columns = self.tickers else: # Multiple stocks - data['Close'] should be a DataFrame prices = data['Close'] if hasattr(prices, 'columns'): prices.columns = self.tickers # Calculate log returns returns = np.log(prices / prices.shift(1)).dropna() # Store data self.prices = prices self.returns = returns logger.info(f"Downloaded {len(returns)} observations") logger.info(f"Date range: {returns.index[0]} to {returns.index[-1]}") return True except Exception as e: logger.error(f"Error downloading data: {str(e)}") return False def create_garch_model_specifications(self) -> List[Dict]: """ Create GARCH model specifications Returns: List[Dict]: List of GARCH model specifications """ models = [] model_configs = [ ('GARCH', [(1,1), (1,2), (2,1), (2,2)]), ('EGARCH', [(1,1), (1,2), (2,1)]), ('TGARCH', [(1,1), (1,2), (2,1)]) ] distributions = ['normal', 't', 'ged'] for model_type, param_combos in model_configs: for p, q in param_combos: for dist in distributions: model_spec = { 'model_type': model_type, 'p': p, 'q': q, 'distribution': dist, 'name': f"{model_type}({p},{q})_{dist}" } if model_type == 'EGARCH': model_spec['vol_type'] = 'EGARCH' elif model_type == 'TGARCH': model_spec['vol_type'] = 'GARCH' model_spec['power'] = 2.0 else: model_spec['vol_type'] = 'Garch' models.append(model_spec) return models def estimate_dcc_parameters(self, max_iter=1000, tolerance=1e-6) -> Dict: """ Estimate DCC parameters using maximum likelihood estimation This implements the two-step DCC estimation procedure: 1. Estimate univariate GARCH models (already done in ensemble fitting) 2. Estimate DCC parameters (alpha, beta) using standardized residuals Args: max_iter (int): Maximum number of iterations tolerance (float): Tolerance for convergence Returns: Dict: Dictionary containing the DCC parameters """ logger.info("Estimating DCC parameters...") if self.standardized_residuals is None: raise ValueError("Must fit ensemble models first") # Calculate unconditional correlation matrix self.unconditional_correlation = self.standardized_residuals.corr().values logger.info("Unconditional correlation matrix estimated") # Initial parameter values initial_params = [0.01, 0.95] # alpha, beta # Bounds for parameters bounds = [(1e-6, 1-1e-6), (1e-6, 1-1e-6)] # alpha, beta must be positive and sum < 1 # Constraint: alpha + beta < 1 def constraint_func(params): return 1 - (params[0] + params[1]) constraint = {'type': 'ineq', 'fun': constraint_func} # Objective function: negative log-likelihood def negative_log_likelihood(params): try: alpha, beta = params return -self.dcc_log_likelihood(alpha, beta) except: return 1e10 # Return large value if calculation fails # Optimize result = minimize(negative_log_likelihood, initial_params, method='SLSQP', bounds=bounds, constraints=constraint, options={'maxiter': max_iter, 'ftol': tolerance}) if result.success: self.dcc_params = {'alpha': result.x[0], 'beta': result.x[1]} logger.info(f"DCC parameters estimated successfully:") logger.info(f"Alpha: {self.dcc_params['alpha']:.6f}") logger.info(f"Beta: {self.dcc_params['beta']:.6f}") logger.info(f"Alpha + Beta: {sum(result.x):.6f}") # Generate dynamic correlation matrices self.generate_dynamic_correlations() else: logger.error("DCC parameter estimation failed!") logger.error(f"Reason: {result.message}") # Use simple fallback self.dcc_params = {'alpha': 0.01, 'beta': 0.95} self.generate_dynamic_correlations() return self.dcc_params def dcc_log_likelihood(self, alpha, beta) -> float: """ Calculate DCC log-likelihood for given parameters Args: alpha (float): Alpha parameter beta (float): Beta parameter Returns: float: The DCC log-likelihood """ n_obs = len(self.standardized_residuals) # Initialize Q matrix with unconditional correlation Q_bar = self.unconditional_correlation.copy() Q_t = Q_bar.copy() log_likelihood = 0.0 for t in range(n_obs): # Get standardized residuals for time t z_t = self.standardized_residuals.iloc[t].values.reshape(-1, 1) # Update Q matrix if t > 0: z_prev = self.standardized_residuals.iloc[t-1].values.reshape(-1, 1) Q_t = (1 - alpha - beta) * Q_bar + alpha * (z_prev @ z_prev.T) + beta * Q_t # Calculate correlation matrix D_t = np.sqrt(np.diag(Q_t)) D_t_inv = np.diag(1.0 / D_t) R_t = D_t_inv @ Q_t @ D_t_inv # Ensure positive definiteness try: R_t_inv = np.linalg.inv(R_t) det_R_t = np.linalg.det(R_t) if det_R_t <= 0: return -1e10 # Add to log-likelihood log_likelihood += -0.5 * (np.log(det_R_t) + z_t.T @ R_t_inv @ z_t) except np.linalg.LinAlgError: return -1e10 return float(log_likelihood) def generate_dynamic_correlations(self) -> None: """ Generate time-varying correlation matrices using estimated DCC parameters """ logger.info("Generating dynamic correlation matrices...") alpha = self.dcc_params['alpha'] beta = self.dcc_params['beta'] n_assets = len(self.tickers) n_obs = len(self.standardized_residuals) # Initialize storage self.q_matrices = np.zeros((n_obs, n_assets, n_assets)) self.dynamic_correlations = np.zeros((n_obs, n_assets, n_assets)) # Initialize Q matrix Q_bar = self.unconditional_correlation.copy() Q_t = Q_bar.copy() for t in range(n_obs): # Get standardized residuals for time t z_t = self.standardized_residuals.iloc[t].values.reshape(-1, 1) # Update Q matrix if t > 0: z_prev = self.standardized_residuals.iloc[t-1].values.reshape(-1, 1) Q_t = (1 - alpha - beta) * Q_bar + alpha * (z_prev @ z_prev.T) + beta * Q_t # Store Q matrix self.q_matrices[t] = Q_t.copy() # Calculate correlation matrix D_t = np.sqrt(np.diag(Q_t)) D_t_inv = np.diag(1.0 / np.maximum(D_t, 1e-8)) # Avoid division by zero R_t = D_t_inv @ Q_t @ D_t_inv # Ensure correlation matrix properties np.fill_diagonal(R_t, 1.0) # Diagonal elements = 1 # Store correlation matrix self.dynamic_correlations[t] = R_t.copy() logger.info(f"Generated {n_obs} dynamic correlation matrices") # Log some statistics avg_corr = np.mean([self.dynamic_correlations[t][0, 1] for t in range(n_obs) if n_assets > 1]) if n_assets > 1: logger.info(f"Average correlation (first pair): {avg_corr:.4f}") def forecast_dcc_correlation(self, horizon=1) -> np.ndarray: """ Forecast correlation matrix using DCC model Args: horizon (int): The horizon for the forecast Returns: np.ndarray: The forecasted correlation matrix """ if self.dynamic_correlations is None: raise ValueError("Must estimate DCC parameters first") alpha = self.dcc_params['alpha'] beta = self.dcc_params['beta'] # Use last Q matrix as starting point Q_t = self.q_matrices[-1].copy() # Get last standardized residuals z_last = self.standardized_residuals.iloc[-1].values.reshape(-1, 1) # Forecast Q matrix Q_forecast = (1 - alpha - beta) * self.unconditional_correlation + \ alpha * (z_last @ z_last.T) + beta * Q_t # For multi-step ahead, use persistence if horizon > 1: Q_forecast = (1 - beta**horizon) * self.unconditional_correlation + beta**horizon * Q_forecast # Convert to correlation matrix D_forecast = np.sqrt(np.diag(Q_forecast)) D_forecast_inv = np.diag(1.0 / np.maximum(D_forecast, 1e-8)) R_forecast = D_forecast_inv @ Q_forecast @ D_forecast_inv # Ensure correlation matrix properties np.fill_diagonal(R_forecast, 1.0) return R_forecast # Inherit and adapt methods from CCC for DCC-specific functionality def fit_garch_model(self, returns, model_spec, verbose=False) -> Optional[arch_model]: """ Fit GARCH model with given specification (inherited from CCC) Args: returns (pd.Series): The returns to fit the model to model_spec (Dict): The model specification verbose (bool): Whether to print verbose output Returns: Optional[arch.arch_model.results.results.Results]: The fitted model """ try: if model_spec['model_type'] == 'GARCH': model = arch_model(returns, vol='Garch', p=model_spec['p'], q=model_spec['q'], dist=model_spec['distribution']) elif model_spec['model_type'] == 'EGARCH': model = arch_model(returns, vol='EGARCH', p=model_spec['p'], q=model_spec['q'], dist=model_spec['distribution']) elif model_spec['model_type'] == 'TGARCH': model = arch_model(returns, vol='GARCH', p=model_spec['p'], q=model_spec['q'], power=model_spec.get('power', 2.0), dist=model_spec['distribution']) result = model.fit(disp='off', show_warning=False) if verbose: logger.info(f"Successfully fitted {model_spec['name']}") return result except Exception as e: if verbose: logger.warning(f"Failed to fit {model_spec['name']}: {str(e)}") return None def calculate_var(self, returns_forecast, volatility_forecast, confidence_level=0.05, distribution='normal', dist_params=None) -> float: """ Calculate Value at Risk using GARCH forecasts (inherited from CCC) Args: returns_forecast (float): The forecasted return volatility_forecast (float): The forecasted volatility confidence_level (float): The confidence level distribution (str): The distribution to use dist_params (Dict): The distribution parameters Returns: float: The Value at Risk """ if distribution == 'normal': z_score = stats.norm.ppf(confidence_level) elif distribution == 't': if dist_params is not None and 'nu' in dist_params: nu = dist_params['nu'] z_score = stats.t.ppf(confidence_level, df=nu) else: z_score = stats.norm.ppf(confidence_level) elif distribution == 'ged': if dist_params is not None and 'nu' in dist_params: nu = dist_params['nu'] z_score = stats.gennorm.ppf(confidence_level, beta=nu) else: z_score = stats.norm.ppf(confidence_level) else: z_score = stats.norm.ppf(confidence_level) var = -(returns_forecast + z_score * volatility_forecast) return var def asymmetric_quantile_loss(self, actual_return, var_forecast, confidence_level=0.05) -> float: """ Asymmetric quantile loss function (inherited from CCC) Args: actual_return (float): The actual return var_forecast (float): The forecasted Value at Risk confidence_level (float): The confidence level Returns: float: The asymmetric quantile loss """ if hasattr(actual_return, 'item'): actual_return = actual_return.item() if hasattr(var_forecast, 'item'): var_forecast = var_forecast.item() actual_loss = -actual_return var_loss = var_forecast exceedance = actual_loss > var_loss if exceedance: loss = confidence_level * (actual_loss - var_loss) else: loss = (1 - confidence_level) * (var_loss - actual_loss) return loss def regulatory_loss_function(self, actual_return, var_forecast, confidence_level=0.05) -> float: """ Regulatory-style loss function (inherited from CCC) Args: actual_return (float): The actual return var_forecast (float): The forecasted Value at Risk confidence_level (float): The confidence level Returns: float: The regulatory loss """ if hasattr(actual_return, 'item'): actual_return = actual_return.item() if hasattr(var_forecast, 'item'): var_forecast = var_forecast.item() actual_loss = -actual_return var_loss = var_forecast exceedance = actual_loss > var_loss base_loss = self.asymmetric_quantile_loss(actual_return, var_forecast, confidence_level) if exceedance: excess_magnitude = (actual_loss - var_loss) / var_loss if var_loss != 0 else 0 magnitude_penalty = excess_magnitude ** 2 base_loss += magnitude_penalty return base_loss def christoffersen_test(self, exceedances, confidence_level=0.05) -> Dict: """ Christoffersen conditional coverage test (inherited from CCC) Args: exceedances (List[int]): The exceedances confidence_level (float): The confidence level Returns: Dict: The Christoffersen test results """ exceedances = np.array(exceedances, dtype=int) n = len(exceedances) n_exceedances = np.sum(exceedances) if n_exceedances == 0 or n_exceedances == n: return { 'lr_uc': np.nan, 'lr_cc': np.nan, 'lr_ind': np.nan, 'p_value_uc': np.nan, 'p_value_cc': np.nan, 'p_value_ind': np.nan, 'status': 'Cannot compute (extreme exceedance rate)' } p_obs = n_exceedances / n p_exp = confidence_level # Unconditional coverage test lr_uc = -2 * (n_exceedances * np.log(p_exp) + (n - n_exceedances) * np.log(1 - p_exp) - n_exceedances * np.log(p_obs) - (n - n_exceedances) * np.log(1 - p_obs)) # Independence test n00 = n01 = n10 = n11 = 0 for i in range(n - 1): if exceedances[i] == 0 and exceedances[i + 1] == 0: n00 += 1 elif exceedances[i] == 0 and exceedances[i + 1] == 1: n01 += 1 elif exceedances[i] == 1 and exceedances[i + 1] == 0: n10 += 1 elif exceedances[i] == 1 and exceedances[i + 1] == 1: n11 += 1 n_0 = n00 + n01 n_1 = n10 + n11 if n_0 > 0 and n_1 > 0 and n01 > 0 and n10 > 0: p01 = n01 / n_0 p11 = n11 / n_1 lr_ind = -2 * (n01 * np.log(p_obs) + n11 * np.log(p_obs) + n00 * np.log(1 - p_obs) + n10 * np.log(1 - p_obs) - n01 * np.log(p01) - n11 * np.log(p11) - n00 * np.log(1 - p01) - n10 * np.log(1 - p11)) else: lr_ind = np.nan lr_cc = lr_uc + lr_ind if not np.isnan(lr_ind) else np.nan p_value_uc = 1 - stats.chi2.cdf(lr_uc, df=1) if not np.isnan(lr_uc) else np.nan p_value_ind = 1 - stats.chi2.cdf(lr_ind, df=1) if not np.isnan(lr_ind) else np.nan p_value_cc = 1 - stats.chi2.cdf(lr_cc, df=2) if not np.isnan(lr_cc) else np.nan return { 'lr_uc': lr_uc, 'lr_cc': lr_cc, 'lr_ind': lr_ind, 'p_value_uc': p_value_uc, 'p_value_cc': p_value_cc, 'p_value_ind': p_value_ind, 'exceedance_rate': p_obs, 'status': 'computed' } def rolling_var_backtest_single_model(self, returns, model_spec, min_train_window=500, rebalance_freq=REBALANCE_FREQUENCY) -> pd.DataFrame: """ Rolling window VaR backtesting for a single model (inherited from CCC) Args: returns (pd.Series): The returns to backtest model_spec (Dict): The model specification """ n_obs = len(returns) initial_train_size = max(int(n_obs * self.train_ratio), min_train_window) results = [] for i in range(initial_train_size, n_obs, rebalance_freq): train_start = 0 train_end = i train_returns = returns.iloc[train_start:train_end] try: garch_result = self.fit_garch_model(train_returns, model_spec, verbose=False) if garch_result is None: continue forecast = garch_result.forecast(horizon=1) vol_forecast = float(np.sqrt(forecast.variance.iloc[-1, 0])) expected_return = float(train_returns.mean()) dist_params = {} if hasattr(garch_result, 'params'): if 'nu' in garch_result.params.index: dist_params['nu'] = garch_result.params['nu'] var_forecast = self.calculate_var(expected_return, vol_forecast, self.confidence_level, model_spec['distribution'], dist_params) for j in range(min(rebalance_freq, n_obs - i)): if i + j < n_obs: actual_return = float(returns.iloc[i + j]) actual_loss = -actual_return exceedance = actual_loss > var_forecast loss_asym = self.asymmetric_quantile_loss(actual_return, var_forecast, self.confidence_level) loss_reg = self.regulatory_loss_function(actual_return, var_forecast, self.confidence_level) results.append({ 'date': returns.index[i + j], 'var_forecast': var_forecast, 'actual_return': actual_return, 'exceedance': bool(exceedance), 'loss_asymmetric': loss_asym, 'loss_regulatory': loss_reg }) except Exception as e: logger.error(f"Error fitting model {model_spec['name']}: {str(e)}") continue return pd.DataFrame(results) def select_best_models_by_quality(self, all_model_results, max_models=5) -> Dict: """ Select the best models based on quality scores (inherited from CCC) Args: all_model_results (Dict): The model results max_models (int): The maximum number of models to select Returns: Dict: The best models """ if len(all_model_results) == 0: return {} # If we have fewer models than max_models, return all if len(all_model_results) <= max_models: logger.info(f"Only {len(all_model_results)} models available - using all") return all_model_results logger.info(f"Selecting best {max_models} models from {len(all_model_results)} available models") # Calculate quality scores using the same approach as ensemble weighting model_names = list(all_model_results.keys()) alpha, beta, gamma, delta = 0.4, 0.3, 0.2, 0.1 # Same weights as ensemble method uc_pvalues = [] cc_pvalues = [] reg_losses = [] exc_rates = [] asym_losses = [] for model_name in model_names: model_data = all_model_results[model_name] uc_pvalues.append(model_data['christoffersen']['p_value_uc']) cc_pvalues.append(model_data['christoffersen']['p_value_cc']) reg_losses.append(model_data['avg_loss_regulatory']) exc_rates.append(abs(model_data['exceedance_rate'] - self.confidence_level)) asym_losses.append(model_data['avg_loss_asymmetric']) uc_pvalues = np.array(uc_pvalues) cc_pvalues = np.array(cc_pvalues) reg_losses = np.array(reg_losses) exc_rates = np.array(exc_rates) asym_losses = np.array(asym_losses) # Handle NaN values uc_pvalues = np.where(np.isnan(uc_pvalues), 1e-6, uc_pvalues) cc_pvalues = np.where(np.isnan(cc_pvalues), 1e-6, cc_pvalues) # Calculate component scores stat_scores = np.sqrt(np.maximum(uc_pvalues, 1e-6) * np.maximum(cc_pvalues, 1e-6)) stat_weights = stat_scores / np.sum(stat_scores) reg_scores = 1.0 / (reg_losses + 1e-10) reg_weights = reg_scores / np.sum(reg_scores) coverage_scores = 1.0 / (exc_rates + 1e-10) coverage_weights = coverage_scores / np.sum(coverage_scores) asym_scores = 1.0 / (asym_losses + 1e-10) asym_weights = asym_scores / np.sum(asym_scores) # Calculate final quality scores final_scores = (alpha * stat_weights + beta * reg_weights + gamma * coverage_weights + delta * asym_weights) # Get indices of top models top_indices = np.argsort(final_scores)[::-1][:max_models] # Sort descending, take top N # Select best models best_models = {} selected_names = [] for idx in top_indices: model_name = model_names[idx] best_models[model_name] = all_model_results[model_name] selected_names.append(model_name) logger.info(f"Selected models: {selected_names}") logger.info(f"Quality scores: {[f'{final_scores[idx]:.4f}' for idx in top_indices]}") return best_models def calculate_ensemble_weights_multicriteria(self, acceptable_models, alpha=0.4, beta=0.3, gamma=0.2, delta=0.1) -> Tuple[Dict, Dict]: """ Multi-criteria ensemble weighting (inherited from CCC) Args: acceptable_models (Dict): The acceptable models alpha (float): The weight for the statistical component beta (float): The weight for the regulatory component gamma (float): The weight for the coverage component delta (float): The weight for the asymmetric component Returns: Tuple[Dict, Dict]: The ensemble weights and the weight breakdown """ model_names = list(acceptable_models.keys()) uc_pvalues = [] cc_pvalues = [] reg_losses = [] exc_rates = [] asym_losses = [] for model_name in model_names: model_data = acceptable_models[model_name] uc_pvalues.append(model_data['christoffersen']['p_value_uc']) cc_pvalues.append(model_data['christoffersen']['p_value_cc']) reg_losses.append(model_data['avg_loss_regulatory']) exc_rates.append(abs(model_data['exceedance_rate'] - self.confidence_level)) asym_losses.append(model_data['avg_loss_asymmetric']) uc_pvalues = np.array(uc_pvalues) cc_pvalues = np.array(cc_pvalues) reg_losses = np.array(reg_losses) exc_rates = np.array(exc_rates) asym_losses = np.array(asym_losses) uc_pvalues = np.where(np.isnan(uc_pvalues), 1e-6, uc_pvalues) cc_pvalues = np.where(np.isnan(cc_pvalues), 1e-6, cc_pvalues) # Calculate component weights stat_scores = np.sqrt(np.maximum(uc_pvalues, 1e-6) * np.maximum(cc_pvalues, 1e-6)) stat_weights = stat_scores / np.sum(stat_scores) reg_scores = 1.0 / (reg_losses + 1e-10) reg_weights = reg_scores / np.sum(reg_scores) coverage_scores = 1.0 / (exc_rates + 1e-10) coverage_weights = coverage_scores / np.sum(coverage_scores) asym_scores = 1.0 / (asym_losses + 1e-10) asym_weights = asym_scores / np.sum(asym_scores) # Combine weights final_weights = (alpha * stat_weights + beta * reg_weights + gamma * coverage_weights + delta * asym_weights) final_weights = final_weights / np.sum(final_weights) weight_breakdown = { 'statistical': dict(zip(model_names, stat_weights)), 'regulatory': dict(zip(model_names, reg_weights)), 'coverage': dict(zip(model_names, coverage_weights)), 'asymmetric': dict(zip(model_names, asym_weights)), 'final': dict(zip(model_names, final_weights)) } return dict(zip(model_names, final_weights)), weight_breakdown def fit_ensemble_univariate_model(self, returns, ticker) -> Dict: """ Fit ensemble of univariate GARCH models for a single asset using the same methodology as CCC Args: returns (pd.Series): The returns to fit the model to ticker (str): The ticker of the asset Returns: Dict: The ensemble model information """ logger.info(f"Fitting ensemble models for {ticker}...") model_specs = self.create_garch_model_specifications() all_model_results = {} # Test each model specification with rolling window VaR backtesting (same as CCC) for i, model_spec in enumerate(model_specs): try: model_results = self.rolling_var_backtest_single_model( returns, model_spec, min_train_window=500, rebalance_freq=REBALANCE_FREQUENCY ) if len(model_results) > 0: exceedances = model_results['exceedance'].values total_obs = len(model_results) total_exceedances = np.sum(exceedances) exceedance_rate = total_exceedances / total_obs avg_loss_asym = model_results['loss_asymmetric'].mean() avg_loss_reg = model_results['loss_regulatory'].mean() christoffersen_results = self.christoffersen_test(exceedances, self.confidence_level) all_model_results[model_spec['name']] = { 'model_spec': model_spec, 'results_df': model_results, 'total_obs': total_obs, 'total_exceedances': total_exceedances, 'exceedance_rate': exceedance_rate, 'avg_loss_asymmetric': avg_loss_asym, 'avg_loss_regulatory': avg_loss_reg, 'christoffersen': christoffersen_results } except Exception as e: logger.error(f"Error fitting model {model_spec['name']}: {str(e)}") continue logger.info(f"Successfully tested {len(all_model_results)} models for {ticker}") # Filter acceptable models using same criteria as CCC acceptable_models = {} for model_name, results in all_model_results.items(): uc_pvalue = results['christoffersen']['p_value_uc'] cc_pvalue = results['christoffersen']['p_value_cc'] if (not np.isnan(uc_pvalue) and not np.isnan(cc_pvalue) and uc_pvalue > 0.05 and cc_pvalue > 0.05): acceptable_models[model_name] = results if len(acceptable_models) == 0: logger.warning(f"No models pass coverage tests for {ticker} - selecting best 5 models based on quality scores") # Use same quality-based selection as CCC best_models = self.select_best_models_by_quality(all_model_results, max_models=5) acceptable_models = best_models logger.info(f"Using {len(acceptable_models)} models for {ticker} ensemble") # Calculate ensemble weights using same multi-criteria approach as CCC weight_dict, weight_breakdown = self.calculate_ensemble_weights_multicriteria(acceptable_models) # Fit final models on full data (same as CCC) final_fitted_models = {} for model_name in acceptable_models.keys(): model_spec = acceptable_models[model_name]['model_spec'] try: fitted_model = self.fit_garch_model(returns, model_spec, verbose=False) if fitted_model is not None: final_fitted_models[model_name] = fitted_model except Exception as e: continue ensemble_info = { 'acceptable_models': acceptable_models, 'fitted_models': final_fitted_models, 'weights': weight_dict, 'weight_breakdown': weight_breakdown, 'ticker': ticker, 'validation_summary': { 'total_models_tested': len(model_specs), 'successful_models': len(all_model_results), 'acceptable_models': len(acceptable_models), 'final_fitted_models': len(final_fitted_models) } } return ensemble_info def fit_all_ensemble_models(self) -> None: """Fit ensemble models for all assets""" logger.info("Fitting ensemble models for all assets...") for ticker in self.tickers: self.ensemble_models[ticker] = self.fit_ensemble_univariate_model( self.returns[ticker], ticker ) # Generate ensemble volatilities and residuals for ticker in self.tickers: self.conditional_volatilities[ticker] = self.get_ensemble_conditional_volatility( self.ensemble_models[ticker], self.returns[ticker] ) # Generate standardized residuals for DCC estimation residuals_dict = {} for ticker in self.tickers: residuals_dict[ticker] = self.get_ensemble_standardized_residuals( self.ensemble_models[ticker], self.returns[ticker] ) self.standardized_residuals = pd.DataFrame(residuals_dict) logger.info("All ensemble models fitted successfully") def get_ensemble_conditional_volatility(self, ensemble_info, returns) -> pd.Series: """ Get ensemble conditional volatility (inherited from CCC) Args: ensemble_info (Dict): The ensemble information returns (pd.Series): The returns to get the conditional volatility for Returns: pd.Series: The ensemble conditional volatility """ weighted_volatility = pd.Series(0.0, index=returns.index) total_weight = 0.0 for model_name, fitted_model in ensemble_info['fitted_models'].items(): try: weight = ensemble_info['weights'][model_name] conditional_vol = fitted_model.conditional_volatility aligned_vol = conditional_vol.reindex(returns.index, fill_value=conditional_vol.mean()) weighted_volatility += weight * aligned_vol total_weight += weight except: continue if total_weight > 0: weighted_volatility = weighted_volatility / total_weight return weighted_volatility def get_ensemble_standardized_residuals(self, ensemble_info, returns) -> pd.Series: """ Get ensemble standardized residuals (inherited from CCC) Args: ensemble_info (Dict): The ensemble information returns (pd.Series): The returns to get the standardized residuals for Returns: pd.Series: The ensemble standardized residuals """ weighted_residuals = pd.Series(0.0, index=returns.index) total_weight = 0.0 for model_name, fitted_model in ensemble_info['fitted_models'].items(): try: weight = ensemble_info['weights'][model_name] residuals = fitted_model.resid conditional_vol = fitted_model.conditional_volatility std_residuals = residuals / conditional_vol aligned_residuals = std_residuals.reindex(returns.index, fill_value=0) weighted_residuals += weight * aligned_residuals total_weight += weight except: continue if total_weight > 0: weighted_residuals = weighted_residuals / total_weight return weighted_residuals def forecast_ensemble_volatility(self, ensemble_info, horizon=1) -> float: """ Generate volatility forecast from ensemble (inherited from CCC) Args: ensemble_info (Dict): The ensemble information """ weighted_forecast = 0.0 total_weight = 0.0 for model_name, fitted_model in ensemble_info['fitted_models'].items(): try: weight = ensemble_info['weights'][model_name] forecast = fitted_model.forecast(horizon=horizon) vol_forecast = np.sqrt(forecast.variance.iloc[-1, 0]) weighted_forecast += weight * vol_forecast total_weight += weight except: continue if total_weight > 0: weighted_forecast = weighted_forecast / total_weight return weighted_forecast def forecast_covariance_matrix(self, horizon=1) -> Tuple[np.ndarray, List[float], np.ndarray]: """ Forecast covariance matrix using DCC model (Difference from CCC: uses dynamic correlations) Args: horizon (int): The horizon for the forecast Returns: Tuple[np.ndarray, List[float], np.ndarray]: The forecasted covariance matrix, volatilities, and correlation matrix """ logger.info(f"Forecasting DCC covariance matrix for {horizon} period(s) ahead...") # Get volatility forecasts vol_forecasts = [] for ticker in self.tickers: vol_forecast = self.forecast_ensemble_volatility( self.ensemble_models[ticker], horizon ) vol_forecasts.append(vol_forecast) logger.info(f"{ticker} forecasted volatility: {vol_forecast:.4f}") # Get correlation forecast using DCC correlation_forecast = self.forecast_dcc_correlation(horizon) # Construct covariance matrix vol_matrix = np.diag(vol_forecasts) forecast_cov = vol_matrix @ correlation_forecast @ vol_matrix logger.info("DCC covariance matrix forecasted") return forecast_cov, vol_forecasts, correlation_forecast def optimize_portfolio(self) -> Dict: """ Portfolio optimization using DCC ensemble model Returns: Dict: The optimization results """ logger.info("Optimizing portfolio with DCC model...") # Get forecasted covariance matrix forecast_cov, vol_forecasts, corr_forecast = self.forecast_covariance_matrix() logger.info("Forecasted volatilities:") for i, ticker in enumerate(self.tickers): logger.info(f" {ticker}: {vol_forecasts[i]:.4f}") # Minimum variance portfolio n_assets = len(self.tickers) def portfolio_variance(weights, cov_matrix): return weights.T @ cov_matrix @ weights # Constraints and bounds constraints = {'type': 'eq', 'fun': lambda w: np.sum(w) - 1} bounds = [(0, 1) for _ in range(n_assets)] initial_weights = np.array([1/n_assets] * n_assets) # Optimize result = minimize(portfolio_variance, initial_weights, args=(forecast_cov,), method='SLSQP', bounds=bounds, constraints=constraints, options={'maxiter': 1000, 'ftol': 1e-10}) if result.success: optimal_weights = result.x min_variance = result.fun logger.info("DCC portfolio optimization successful!") logger.info("Optimal weights:") for i, ticker in enumerate(self.tickers): logger.info(f" {ticker}: {optimal_weights[i]:.4f}") logger.info(f"Portfolio Variance: {min_variance:.6f}") logger.info(f"Portfolio Volatility: {np.sqrt(min_variance):.4f}") # Store optimization results self.optimization_results = { 'optimal_weights': optimal_weights, 'portfolio_variance': min_variance, 'portfolio_volatility': np.sqrt(min_variance), 'forecasted_volatilities': vol_forecasts, 'forecasted_covariance': forecast_cov, 'forecasted_correlation': corr_forecast, 'current_weights': self.weights, 'optimization_success': True, 'model_type': 'DCC' } else: logger.error("DCC portfolio optimization failed!") logger.error(f"Reason: {result.message}") # Use current weights as fallback optimal_weights = np.array(self.weights) min_variance = portfolio_variance(optimal_weights, forecast_cov) self.optimization_results = { 'optimal_weights': optimal_weights, 'portfolio_variance': min_variance, 'portfolio_volatility': np.sqrt(min_variance), 'forecasted_volatilities': vol_forecasts, 'forecasted_covariance': forecast_cov, 'forecasted_correlation': corr_forecast, 'current_weights': self.weights, 'optimization_success': False, 'model_type': 'DCC' } return self.optimization_results def run_complete_analysis(self) -> Dict: """ Run complete DCC analysis pipeline Returns: Dict: The analysis results """ logger.info("STARTING DCC PORTFOLIO ANALYSIS") logger.info("\n") try: # Step 1: Download data if not self.download_portfolio_data(): raise Exception("Failed to download portfolio data") # Step 2: Fit ensemble models self.fit_all_ensemble_models() # Step 3: Estimate DCC parameters self.estimate_dcc_parameters() # Step 4: Optimize portfolio self.optimize_portfolio() # Step 5: Generate results results = self.generate_results() logger.info("DCC PORTFOLIO ANALYSIS COMPLETED SUCCESSFULLY") logger.info("\n") return results except Exception as e: logger.error(f"DCC Analysis failed: {str(e)}") raise def generate_results(self) -> Dict: """ Generate DCC modeling results Returns: Dict: The results """ logger.info("Generating DCC results...") # Generate forecasts forecast_results = self.volatility_forecast() # Create plots plot_files = self.create_plots() # Portfolio summary portfolio_summary = { 'total_value': self.total_value, 'num_stocks': len(self.tickers), 'tickers': self.tickers, 'amounts': self.amounts, 'weights': self.weights, 'model_type': 'DCC', 'data_period': { 'start_date': self.start_date, 'end_date': self.end_date, 'total_observations': len(self.returns) } } # Ensemble model summary ensemble_summary = {} for ticker in self.tickers: ensemble_info = self.ensemble_models[ticker] ensemble_summary[ticker] = { 'total_models_tested': ensemble_info['validation_summary']['total_models_tested'], 'successful_models': ensemble_info['validation_summary']['successful_models'], 'acceptable_models': ensemble_info['validation_summary']['acceptable_models'], 'final_fitted_models': ensemble_info['validation_summary']['final_fitted_models'], 'top_3_models': dict(sorted(ensemble_info['weights'].items(), key=lambda x: x[1], reverse=True)[:3]), 'avg_conditional_volatility': self.conditional_volatilities[ticker].mean(), 'model_performance_stats': { 'avg_exceedance_rate': np.mean([ model_data['exceedance_rate'] for model_data in ensemble_info['acceptable_models'].values() ]), 'avg_regulatory_loss': np.mean([ model_data['avg_loss_regulatory'] for model_data in ensemble_info['acceptable_models'].values() ]), 'avg_uc_pvalue': np.mean([ model_data['christoffersen']['p_value_uc'] for model_data in ensemble_info['acceptable_models'].values() if not np.isnan(model_data['christoffersen']['p_value_uc']) ]) } } # DCC-specific results dcc_results = { 'dcc_parameters': self.dcc_params, 'unconditional_correlation': self.unconditional_correlation.tolist(), 'final_correlation': self.dynamic_correlations[-1].tolist() if self.dynamic_correlations is not None else None, 'average_correlations': np.mean(self.dynamic_correlations, axis=0).tolist() if self.dynamic_correlations is not None else None, 'correlation_time_series': { 'dates': self.returns.index.strftime('%Y-%m-%d').tolist(), 'correlations': self.dynamic_correlations.tolist() if self.dynamic_correlations is not None else None }, 'correlation_analysis': { 'min_correlation': float(np.min(self.dynamic_correlations)) if self.dynamic_correlations is not None else None, 'max_correlation': float(np.max(self.dynamic_correlations)) if self.dynamic_correlations is not None else None, 'correlation_std': float(np.std(self.dynamic_correlations)) if self.dynamic_correlations is not None else None, 'time_varying_correlations': True } } # Risk metrics (same structure as CCC model for compatibility) risk_metrics = { 'correlation_matrix': self.unconditional_correlation.tolist() if self.unconditional_correlation is not None else None, 'average_correlation': np.mean(self.dynamic_correlations) if self.dynamic_correlations is not None else None, 'portfolio_volatility_forecast': self.optimization_results['portfolio_volatility'], 'individual_volatility_forecasts': dict(zip(self.tickers, self.optimization_results['forecasted_volatilities'])), 'comprehensive_forecasts': forecast_results, 'correlation_analysis': { 'min_correlation': float(np.min(self.dynamic_correlations)) if self.dynamic_correlations is not None else None, 'max_correlation': float(np.max(self.dynamic_correlations)) if self.dynamic_correlations is not None else None, 'correlation_std': float(np.std(self.dynamic_correlations)) if self.dynamic_correlations is not None else None }, # DCC-specific additions 'dynamic_correlation_matrices': self.dynamic_correlations.tolist() if self.dynamic_correlations is not None else None, 'dcc_correlation_dynamics': { 'correlation_persistence': self.dcc_params['alpha'] + self.dcc_params['beta'] if self.dcc_params else None, 'shock_sensitivity': self.dcc_params['alpha'] if self.dcc_params else None, 'mean_reversion_speed': 1 - (self.dcc_params['alpha'] + self.dcc_params['beta']) if self.dcc_params else None } } # Optimization results optimization_summary = { 'optimization_successful': self.optimization_results['optimization_success'], 'optimal_weights': dict(zip(self.tickers, self.optimization_results['optimal_weights'])), 'current_weights': dict(zip(self.tickers, self.optimization_results['current_weights'])), 'portfolio_variance': self.optimization_results['portfolio_variance'], 'portfolio_volatility': self.optimization_results['portfolio_volatility'], 'weight_changes': { ticker: self.optimization_results['optimal_weights'][i] - self.optimization_results['current_weights'][i] for i, ticker in enumerate(self.tickers) }, 'rebalancing_analysis': { 'significant_changes': [ (ticker, change) for ticker, change in [(ticker, self.optimization_results['optimal_weights'][i] - self.optimization_results['current_weights'][i]) for i, ticker in enumerate(self.tickers)] if abs(change) > 0.05 ], 'total_turnover': sum([ abs(self.optimization_results['optimal_weights'][i] - self.optimization_results['current_weights'][i]) for i in range(len(self.tickers)) ]) / 2 } } # Compile all results self.modeling_results = { 'portfolio_summary': portfolio_summary, 'ensemble_summary': ensemble_summary, 'risk_metrics': risk_metrics, 'optimization_results': optimization_summary, 'dcc_results': dcc_results, 'plot_files': plot_files, 'model_metadata': { 'confidence_level': self.confidence_level, 'train_ratio': self.train_ratio, 'modeling_approach': 'DCC-GARCH with Ensemble Univariate Models', 'rebalance_frequency': REBALANCE_FREQUENCY, 'generation_timestamp': datetime.now().isoformat(), 'technical_details': { 'garch_specifications': len(self.create_garch_model_specifications()), 'distributions_used': ['normal', 't', 'ged'], 'model_types': ['GARCH', 'EGARCH', 'TGARCH'], 'validation_method': 'Rolling window VaR backtesting with Christoffersen tests', 'ensemble_weighting': 'Multi-criteria: statistical validity (40%), regulatory loss (30%), coverage accuracy (20%), asymmetric loss (10%)', 'correlation_estimation': 'Dynamic Conditional Correlation with time-varying correlation matrices', 'optimization_method': 'Minimum variance with full investment constraint using DCC forecasted correlations', 'dcc_specification': 'Q_t = (1-α-β)Q̄ + α(z_{t-1}z\'_{t-1}) + βQ_{t-1}, R_t = D_t^{-1}Q_tD_t^{-1}' } } } logger.info("DCC results generated successfully") return self.modeling_results def volatility_forecast(self, horizons=[1, 5, 10, 22], confidence_levels=[0.05, 0.10, 0.25]) -> Dict: """ Generate volatility forecasts with multiple horizons and confidence intervals Args: horizons: List of forecast horizons (in days) confidence_levels: List of confidence levels for intervals Returns: Dictionary with forecast results """ logger.info("Generating DCC volatility forecasts...") forecast_results = {} for ticker in self.tickers: logger.info(f"Forecasting volatility for {ticker}...") ticker_forecasts = { 'horizons': horizons, 'confidence_levels': confidence_levels, 'point_forecasts': [], 'confidence_intervals': {}, 'var_forecasts': {} } ensemble_info = self.ensemble_models[ticker] # Initialize confidence intervals storage for conf_level in confidence_levels: ticker_forecasts['confidence_intervals'][conf_level] = { 'lower': [], 'upper': [] } ticker_forecasts['var_forecasts'][conf_level] = [] # Generate forecasts for each horizon for horizon in horizons: # Collect forecasts from all models in ensemble model_forecasts = [] var_forecasts_by_conf = {conf: [] for conf in confidence_levels} for model_name, fitted_model in ensemble_info['fitted_models'].items(): try: weight = ensemble_info['weights'][model_name] # Get model forecast forecast = fitted_model.forecast(horizon=horizon) vol_forecast = np.sqrt(forecast.variance.iloc[-1, 0]) # Add to collection (weighted) for _ in range(int(weight * 1000)): # Weight by adding multiple copies model_forecasts.append(vol_forecast) # Calculate VaR for each confidence level model_spec = ensemble_info['acceptable_models'][model_name]['model_spec'] expected_return = self.returns[ticker].mean() for conf_level in confidence_levels: var_forecast = self.calculate_var( expected_return, vol_forecast, conf_level, model_spec['distribution'] ) for _ in range(int(weight * 1000)): var_forecasts_by_conf[conf_level].append(var_forecast) except Exception as e: continue # Calculate ensemble statistics if model_forecasts: point_forecast = np.mean(model_forecasts) ticker_forecasts['point_forecasts'].append(point_forecast) # Calculate confidence intervals for volatility for conf_level in confidence_levels: lower_bound = np.percentile(model_forecasts, (conf_level/2) * 100) upper_bound = np.percentile(model_forecasts, (1 - conf_level/2) * 100) ticker_forecasts['confidence_intervals'][conf_level]['lower'].append(lower_bound) ticker_forecasts['confidence_intervals'][conf_level]['upper'].append(upper_bound) # VaR forecast var_point = np.mean(var_forecasts_by_conf[conf_level]) ticker_forecasts['var_forecasts'][conf_level].append(var_point) else: # Fallback to single forecast vol_forecast = self.forecast_ensemble_volatility(ensemble_info, horizon) ticker_forecasts['point_forecasts'].append(vol_forecast) for conf_level in confidence_levels: # Simple confidence intervals (±20% of point forecast) margin = vol_forecast * 0.2 ticker_forecasts['confidence_intervals'][conf_level]['lower'].append(vol_forecast - margin) ticker_forecasts['confidence_intervals'][conf_level]['upper'].append(vol_forecast + margin) # VaR calculation expected_return = self.returns[ticker].mean() var_forecast = self.calculate_var(expected_return, vol_forecast, conf_level) ticker_forecasts['var_forecasts'][conf_level].append(var_forecast) forecast_results[ticker] = ticker_forecasts # Portfolio-level forecasts with DCC dynamics portfolio_forecasts = { 'horizons': horizons, 'portfolio_volatility': [], 'portfolio_var': {}, 'correlation_forecasts': [] } for conf_level in confidence_levels: portfolio_forecasts['portfolio_var'][conf_level] = [] for horizon in horizons: # Get individual volatility forecasts for this horizon individual_vols = [forecast_results[ticker]['point_forecasts'][horizons.index(horizon)] for ticker in self.tickers] # Get DCC correlation forecast for this horizon correlation_forecast = self.forecast_dcc_correlation(horizon) portfolio_forecasts['correlation_forecasts'].append(correlation_forecast.tolist()) # Calculate portfolio volatility using DCC correlations vol_matrix = np.diag(individual_vols) portfolio_cov = vol_matrix @ correlation_forecast @ vol_matrix portfolio_vol = np.sqrt(np.array(self.weights) @ portfolio_cov @ np.array(self.weights)) portfolio_forecasts['portfolio_volatility'].append(portfolio_vol) # Portfolio VaR portfolio_return = np.sum([w * self.returns[ticker].mean() for w, ticker in zip(self.weights, self.tickers)]) for conf_level in confidence_levels: portfolio_var = self.calculate_var(portfolio_return, portfolio_vol, conf_level) portfolio_forecasts['portfolio_var'][conf_level].append(portfolio_var) results = { 'individual_forecasts': forecast_results, 'portfolio_forecasts': portfolio_forecasts, 'dcc_correlation_forecasts': { 'horizons': horizons, 'correlation_matrices': portfolio_forecasts['correlation_forecasts'], 'correlation_persistence': self.dcc_params['alpha'] + self.dcc_params['beta'] if self.dcc_params else None }, 'metadata': { 'horizons': horizons, 'confidence_levels': confidence_levels, 'forecast_date': datetime.now().isoformat(), 'base_currency': 'USD', 'model_type': 'DCC-GARCH Ensemble' } } logger.info("DCC volatility forecasting completed") return results def create_plots(self, save_directory='dcc_plots') -> Dict: """ Create visualization plots for DCC model Args: save_directory: Directory to save plots Returns: Dictionary with plot file paths """ logger.info("Creating DCC plots...") # Create plots directory if not os.path.exists(save_directory): os.makedirs(save_directory) plot_files = {} # Plot 1: Cumulative Returns Evolution plt.figure(figsize=(12, 8)) cumulative_returns = self.returns.cumsum().apply(np.exp) for ticker in self.tickers: plt.plot(cumulative_returns.index, cumulative_returns[ticker], label=ticker, linewidth=2) plt.title('Portfolio Cumulative Returns Evolution (DCC Analysis)', fontsize=16, fontweight='bold') plt.xlabel('Date', fontsize=12) plt.ylabel('Cumulative Return', fontsize=12) plt.legend(fontsize=12) plt.grid(True, alpha=0.3) plt.tight_layout() plot_path = os.path.join(save_directory, '01_dcc_cumulative_returns.png') plt.savefig(plot_path, dpi=300, bbox_inches='tight') plt.close() plot_files['cumulative_returns'] = plot_path # Plot 2: Daily Returns plt.figure(figsize=(14, 8)) for ticker in self.tickers: plt.plot(self.returns.index, self.returns[ticker], label=ticker, alpha=0.7, linewidth=0.8) plt.title('Daily Returns Time Series (DCC Analysis)', fontsize=16, fontweight='bold') plt.xlabel('Date', fontsize=12) plt.ylabel('Daily Return', fontsize=12) plt.legend(fontsize=12) plt.grid(True, alpha=0.3) plt.tight_layout() plot_path = os.path.join(save_directory, '02_dcc_daily_returns.png') plt.savefig(plot_path, dpi=300, bbox_inches='tight') plt.close() plot_files['daily_returns'] = plot_path # Plot 3: Ensemble Conditional Volatilities plt.figure(figsize=(14, 8)) for ticker in self.tickers: vol_series = self.conditional_volatilities[ticker] plt.plot(vol_series.index, vol_series, label=f'{ticker} DCC Ensemble Volatility', linewidth=2) plt.title('DCC Ensemble Conditional Volatilities', fontsize=16, fontweight='bold') plt.xlabel('Date', fontsize=12) plt.ylabel('Conditional Volatility', fontsize=12) plt.legend(fontsize=12) plt.grid(True, alpha=0.3) plt.tight_layout() plot_path = os.path.join(save_directory, '03_dcc_ensemble_volatilities.png') plt.savefig(plot_path, dpi=300, bbox_inches='tight') plt.close() plot_files['ensemble_volatilities'] = plot_path # Plot 4: Dynamic Correlations Evolution (DCC-specific) if self.dynamic_correlations is not None and len(self.tickers) > 1: plt.figure(figsize=(14, 10)) # Plot all pairwise correlations n_assets = len(self.tickers) pair_idx = 0 colors = plt.cm.tab10(np.linspace(0, 1, n_assets*(n_assets-1)//2)) for i in range(n_assets): for j in range(i+1, n_assets): corr_series = [self.dynamic_correlations[t][i, j] for t in range(len(self.dynamic_correlations))] plt.plot(self.returns.index, corr_series, label=f'{self.tickers[i]}-{self.tickers[j]}', color=colors[pair_idx], linewidth=2) pair_idx += 1 plt.title('DCC Dynamic Correlations Evolution', fontsize=16, fontweight='bold') plt.xlabel('Date', fontsize=12) plt.ylabel('Correlation Coefficient', fontsize=12) plt.legend(fontsize=10) plt.grid(True, alpha=0.3) plt.ylim(-1, 1) plt.tight_layout() plot_path = os.path.join(save_directory, '04_dcc_dynamic_correlations.png') plt.savefig(plot_path, dpi=300, bbox_inches='tight') plt.close() plot_files['dynamic_correlations'] = plot_path # Plot 5: Correlation Heatmaps Comparison (Unconditional vs Final) if len(self.tickers) > 1: fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(18, 7)) # Unconditional correlation im1 = ax1.imshow(self.unconditional_correlation, cmap='RdBu_r', vmin=-1, vmax=1) ax1.set_title('Unconditional Correlation Matrix', fontsize=14, fontweight='bold') ax1.set_xticks(range(len(self.tickers))) ax1.set_yticks(range(len(self.tickers))) ax1.set_xticklabels(self.tickers) ax1.set_yticklabels(self.tickers) # Add text annotations for i in range(len(self.tickers)): for j in range(len(self.tickers)): ax1.text(j, i, f'{self.unconditional_correlation[i, j]:.3f}', ha='center', va='center', fontsize=10, fontweight='bold') # Final period correlation final_corr = self.dynamic_correlations[-1] if self.dynamic_correlations is not None else self.unconditional_correlation im2 = ax2.imshow(final_corr, cmap='RdBu_r', vmin=-1, vmax=1) ax2.set_title('Final Period DCC Correlation Matrix', fontsize=14, fontweight='bold') ax2.set_xticks(range(len(self.tickers))) ax2.set_yticks(range(len(self.tickers))) ax2.set_xticklabels(self.tickers) ax2.set_yticklabels(self.tickers) # Add text annotations for i in range(len(self.tickers)): for j in range(len(self.tickers)): ax2.text(j, i, f'{final_corr[i, j]:.3f}', ha='center', va='center', fontsize=10, fontweight='bold') # Adjust subplot spacing to make room for colorbar plt.subplots_adjust(left=0.05, right=0.85, top=0.95, bottom=0.05, wspace=0.3) # Add colorbar in a dedicated space cbar_ax = fig.add_axes([0.87, 0.15, 0.02, 0.7]) # [left, bottom, width, height] plt.colorbar(im2, cax=cbar_ax, label='Correlation Coefficient') plot_path = os.path.join(save_directory, '05_dcc_correlation_comparison.png') plt.savefig(plot_path, dpi=300, bbox_inches='tight') plt.close() plot_files['correlation_comparison'] = plot_path # Plot 6: DCC Parameters and Persistence if self.dcc_params: fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6)) # DCC Parameters params = ['Alpha (α)', 'Beta (β)', 'Persistence (α+β)'] values = [self.dcc_params['alpha'], self.dcc_params['beta'], self.dcc_params['alpha'] + self.dcc_params['beta']] colors = ['red', 'blue', 'green'] bars = ax1.bar(params, values, color=colors, alpha=0.7) ax1.set_title('DCC Parameters', fontsize=14, fontweight='bold') ax1.set_ylabel('Parameter Value', fontsize=12) ax1.grid(True, alpha=0.3, axis='y') # Add value labels on bars for bar, value in zip(bars, values): ax1.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.01, f'{value:.4f}', ha='center', va='bottom', fontweight='bold') # Interpretation guide ax2.text(0.1, 0.9, 'DCC Parameter Interpretation:', fontsize=14, fontweight='bold', transform=ax2.transAxes) ax2.text(0.1, 0.8, f'α = {self.dcc_params["alpha"]:.4f} (Shock Sensitivity)', fontsize=12, transform=ax2.transAxes) ax2.text(0.1, 0.7, f'β = {self.dcc_params["beta"]:.4f} (Correlation Persistence)', fontsize=12, transform=ax2.transAxes) ax2.text(0.1, 0.6, f'α + β = {values[2]:.4f} (Overall Persistence)', fontsize=12, transform=ax2.transAxes) if self.dcc_params['alpha'] > 0.05: shock_interp = 'High sensitivity to recent shocks' elif self.dcc_params['alpha'] < 0.01: shock_interp = 'Low sensitivity to recent shocks' else: shock_interp = 'Moderate sensitivity to recent shocks' if values[2] > 0.99: persist_interp = 'Very high correlation persistence' elif values[2] > 0.95: persist_interp = 'High correlation persistence' else: persist_interp = 'Moderate correlation persistence' ax2.text(0.1, 0.4, shock_interp, fontsize=11, style='italic', transform=ax2.transAxes) ax2.text(0.1, 0.3, persist_interp, fontsize=11, style='italic', transform=ax2.transAxes) ax2.set_xlim(0, 1) ax2.set_ylim(0, 1) ax2.axis('off') plt.tight_layout() plot_path = os.path.join(save_directory, '06_dcc_parameters.png') plt.savefig(plot_path, dpi=300, bbox_inches='tight') plt.close() plot_files['dcc_parameters'] = plot_path # Plot 7: Model Weights for each asset for ticker in self.tickers: plt.figure(figsize=(12, 8)) weights = self.ensemble_models[ticker]['weights'] # Get all models and their weights, sort by weight sorted_weights = sorted(weights.items(), key=lambda x: x[1], reverse=True) models, weight_values = zip(*sorted_weights) # Clean model names for display clean_models = [m.replace('_', '\n').replace('(', '\n(') for m in models] colors = plt.cm.viridis(np.linspace(0, 1, len(models))) bars = plt.barh(range(len(models)), weight_values, color=colors) plt.yticks(range(len(models)), clean_models, fontsize=8) plt.xlabel('Weight', fontsize=12) plt.title(f'{ticker} - DCC Ensemble Model Weights', fontsize=16, fontweight='bold') plt.grid(True, alpha=0.3, axis='x') # Add weight values on bars for i, (bar, weight) in enumerate(zip(bars, weight_values)): plt.text(weight + 0.001, i, f'{weight:.3f}', va='center', fontsize=8) plt.tight_layout() plot_path = os.path.join(save_directory, f'07_dcc_model_weights_{ticker}.png') plt.savefig(plot_path, dpi=300, bbox_inches='tight') plt.close() plot_files[f'model_weights_{ticker}'] = plot_path # Plot 8: Volatility Forecasts with DCC forecast_results = self.volatility_forecast() for ticker in self.tickers: plt.figure(figsize=(12, 8)) ticker_data = forecast_results['individual_forecasts'][ticker] horizons = ticker_data['horizons'] point_forecasts = ticker_data['point_forecasts'] # Plot point forecasts plt.plot(horizons, point_forecasts, 'o-', linewidth=2, markersize=8, label='DCC Point Forecast', color='blue') # Collect all y-values to determine appropriate y-axis limits all_y_values = list(point_forecasts) # Plot confidence intervals with improved visibility colors = ['red', 'orange', 'green'] alphas = [0.2, 0.3, 0.4] # Different transparencies for better visibility for i, conf_level in enumerate(ticker_data['confidence_levels']): lower = ticker_data['confidence_intervals'][conf_level]['lower'] upper = ticker_data['confidence_intervals'][conf_level]['upper'] # Add confidence interval bounds to y-values for scaling all_y_values.extend(lower) all_y_values.extend(upper) # Check if confidence intervals are meaningful (not too narrow) interval_width = np.mean([u - l for u, l in zip(upper, lower)]) point_level = np.mean(point_forecasts) # Only plot if interval is at least 1% of the point forecast if interval_width > point_level * 0.01: plt.fill_between(horizons, lower, upper, alpha=alphas[i], color=colors[i], label=f'{(1-conf_level)*100:.0f}% Confidence Interval') # Also plot the bounds as lines for better visibility plt.plot(horizons, lower, '--', color=colors[i], alpha=0.7, linewidth=1) plt.plot(horizons, upper, '--', color=colors[i], alpha=0.7, linewidth=1) else: # If intervals are too narrow, just indicate uncertainty with error bars # Calculate error bar values ensuring they are always positive point_arr = np.array(point_forecasts) lower_arr = np.array(lower) upper_arr = np.array(upper) # Validate confidence intervals if np.any(lower_arr > point_arr): logger.warning(f"Lower confidence bound exceeds point forecast for {ticker}") # Correct invalid bounds lower_arr = np.minimum(lower_arr, point_arr) if np.any(upper_arr < point_arr): logger.warning(f"Upper confidence bound below point forecast for {ticker}") # Correct invalid bounds upper_arr = np.maximum(upper_arr, point_arr) lower_err = np.maximum(0, point_arr - lower_arr) upper_err = np.maximum(0, upper_arr - point_arr) # Only plot error bars if there's meaningful uncertainty if np.any(lower_err > 0) or np.any(upper_err > 0): plt.errorbar(horizons, point_forecasts, yerr=[lower_err, upper_err], fmt='none', color=colors[i], alpha=0.6, capsize=3, label=f'{(1-conf_level)*100:.0f}% Confidence Range') # Set appropriate y-axis limits with padding if all_y_values: y_min, y_max = min(all_y_values), max(all_y_values) y_range = y_max - y_min # Add 10% padding to ensure intervals are visible padding = max(y_range * 0.1, y_max * 0.05) # At least 5% of max value plt.ylim(max(0, y_min - padding), y_max + padding) plt.xlabel('Forecast Horizon (Days)', fontsize=12) plt.ylabel('Volatility', fontsize=12) plt.title(f'{ticker} - DCC Volatility Forecast with Confidence Intervals', fontsize=16, fontweight='bold') plt.legend(fontsize=10, loc='best') plt.grid(True, alpha=0.3) plt.tight_layout() plot_path = os.path.join(save_directory, f'08_dcc_volatility_forecast_{ticker}.png') plt.savefig(plot_path, dpi=300, bbox_inches='tight') plt.close() plot_files[f'volatility_forecast_{ticker}'] = plot_path # Plot 9: Portfolio VaR Forecasts with DCC plt.figure(figsize=(12, 8)) portfolio_data = forecast_results['portfolio_forecasts'] horizons = portfolio_data['horizons'] colors = ['red', 'orange', 'green'] for i, conf_level in enumerate([0.05, 0.10, 0.25]): if conf_level in portfolio_data['portfolio_var']: var_forecasts = portfolio_data['portfolio_var'][conf_level] plt.plot(horizons, var_forecasts, 'o-', linewidth=2, markersize=8, color=colors[i], label=f'DCC VaR {conf_level*100:.0f}%') plt.xlabel('Forecast Horizon (Days)', fontsize=12) plt.ylabel('Value at Risk', fontsize=12) plt.title('Portfolio VaR Forecasts (DCC Model)', fontsize=16, fontweight='bold') plt.legend(fontsize=12) plt.grid(True, alpha=0.3) plt.tight_layout() plot_path = os.path.join(save_directory, '09_dcc_portfolio_var_forecasts.png') plt.savefig(plot_path, dpi=300, bbox_inches='tight') plt.close() plot_files['portfolio_var_forecasts'] = plot_path # Plot 10: Ensemble Validation Summary plt.figure(figsize=(14, 8)) # Collect summary data total_tested = [] successful = [] acceptable = [] final_fitted = [] for ticker in self.tickers: summary = self.ensemble_models[ticker]['validation_summary'] total_tested.append(summary['total_models_tested']) successful.append(summary['successful_models']) acceptable.append(summary['acceptable_models']) final_fitted.append(summary['final_fitted_models']) x = np.arange(len(self.tickers)) width = 0.2 plt.bar(x - 1.5*width, total_tested, width, label='Total Tested', alpha=0.8) plt.bar(x - 0.5*width, successful, width, label='Successful', alpha=0.8) plt.bar(x + 0.5*width, acceptable, width, label='Acceptable (Coverage)', alpha=0.8) plt.bar(x + 1.5*width, final_fitted, width, label='Final Fitted', alpha=0.8) plt.xticks(x, self.tickers, fontsize=12) plt.ylabel('Number of Models', fontsize=12) plt.title('DCC Ensemble Model Validation Summary by Asset', fontsize=16, fontweight='bold') plt.legend(fontsize=10) plt.grid(True, alpha=0.3, axis='y') plt.tight_layout() plot_path = os.path.join(save_directory, '10_dcc_validation_summary.png') plt.savefig(plot_path, dpi=300, bbox_inches='tight') plt.close() plot_files['validation_summary'] = plot_path logger.info(f"Created {len(plot_files)} DCC plots in {save_directory}") return plot_files def run_dcc_portfolio_analysis(portfolio_data, start_date='2014-01-01', end_date=None, confidence_level=0.05, train_ratio=0.8, save_results=True) -> Dict: """ Standalone function to run complete DCC portfolio analysis Args: portfolio_data: List of (ticker, amount) tuples start_date: Start date for historical data end_date: End date for historical data confidence_level: VaR confidence level train_ratio: Training/testing split ratio save_results: Whether to save results to file Returns: Dictionary with DCC modeling results """ # Create DCC model instance model = PortfolioDCCEnsembleModel( portfolio_data=portfolio_data, start_date=start_date, end_date=end_date, confidence_level=confidence_level, train_ratio=train_ratio ) # Run complete analysis results = model.run_complete_analysis() # Save results if requested if save_results: filename = f"dcc_portfolio_analysis_{datetime.now().strftime('%Y%m%d_%H%M%S')}.json" with open(filename, 'w') as f: json.dump(results, f, indent=2, default=str) logger.info(f"DCC results saved to {filename}") return results # Example usage and testing if __name__ == "__main__": # Example portfolio (similar to what users might submit) example_portfolio = [ ('AAPL', 5000), ('GOOGL', 3000), ('MSFT', 4000), ] print("PORTFOLIO MODELING - CCC vs DCC COMPARISON") try: # Run CCC analysis print("\n") print("RUNNING CCC ANALYSIS") ccc_results = run_portfolio_analysis( portfolio_data=example_portfolio, start_date='2010-01-01', end_date='2025-01-01', confidence_level=0.05, train_ratio=0.8, save_results=False ) print("\nCCC Analysis completed successfully!") print(f"Portfolio value: ${ccc_results['portfolio_summary']['total_value']:,.2f}") print(f"Number of stocks: {ccc_results['portfolio_summary']['num_stocks']}") print(f"Portfolio volatility forecast: {ccc_results['optimization_results']['portfolio_volatility']:.6f}") print("\nCCC Optimal weights:") # Handle both numpy array and dictionary formats (caused some troubles...) ccc_weights = ccc_results['optimization_results']['optimal_weights'] if isinstance(ccc_weights, dict): for ticker in ccc_results['portfolio_summary']['tickers']: print(f" {ticker}: {ccc_weights[ticker]:.4f}") else: # numpy array format for i, ticker in enumerate(ccc_results['portfolio_summary']['tickers']): print(f" {ticker}: {ccc_weights[i]:.4f}") # Run DCC analysis print("\n") print("RUNNING DCC ANALYSIS") dcc_results = run_dcc_portfolio_analysis( portfolio_data=example_portfolio, start_date='2010-01-01', end_date='2025-01-01', confidence_level=0.05, train_ratio=0.8, save_results=False ) print("\nDCC Analysis completed successfully!") print(f"Portfolio value: ${dcc_results['portfolio_summary']['total_value']:,.2f}") print(f"Number of stocks: {dcc_results['portfolio_summary']['num_stocks']}") print(f"Portfolio volatility forecast: {dcc_results['optimization_results']['portfolio_volatility']:.6f}") print("\nDCC Optimal weights:") # Handle both numpy array and dictionary formats dcc_weights = dcc_results['optimization_results']['optimal_weights'] if isinstance(dcc_weights, dict): for ticker in dcc_results['portfolio_summary']['tickers']: print(f" {ticker}: {dcc_weights[ticker]:.4f}") else: # numpy array format for i, ticker in enumerate(dcc_results['portfolio_summary']['tickers']): print(f" {ticker}: {dcc_weights[i]:.4f}") print("\nDCC Parameters:") print(f"Alpha: {dcc_results['dcc_results']['dcc_parameters']['alpha']:.6f}") print(f"Beta: {dcc_results['dcc_results']['dcc_parameters']['beta']:.6f}") print(f"Alpha + Beta: {dcc_results['dcc_results']['dcc_parameters']['alpha'] + dcc_results['dcc_results']['dcc_parameters']['beta']:.6f}") # Model Comparison print("\n") print("MODEL COMPARISON") # 1. Portfolio Risk Comparison ccc_vol = ccc_results['optimization_results']['portfolio_volatility'] dcc_vol = dcc_results['optimization_results']['portfolio_volatility'] vol_diff = dcc_vol - ccc_vol vol_improvement = (vol_diff / ccc_vol) * 100 print(f"\n1. PORTFOLIO RISK METRICS:") print(f"CCC Portfolio Volatility: {ccc_vol:.6f}") print(f"DCC Portfolio Volatility: {dcc_vol:.6f}") print(f"Difference (DCC - CCC): {vol_diff:.6f}") print(f"Relative Change: {vol_improvement:+.2f}%") if abs(vol_diff) > 0.001: # Significant difference threshold better_model = "DCC" if vol_diff < 0 else "CCC" print(f"*** {better_model} provides {'lower' if vol_diff < 0 else 'higher'} portfolio risk ***") else: print(" *** Models provide similar risk estimates ***") # 2. Weight Allocation Comparison print(f"\n2. PORTFOLIO ALLOCATION COMPARISON:") print(f"{'Asset':<8} {'CCC Weight':<12} {'DCC Weight':<12} {'Difference':<12} {'Change %':<10}") print(f"{'-'*8} {'-'*12} {'-'*12} {'-'*12} {'-'*10}") total_turnover = 0 significant_changes = [] # Get weights in consistent format ccc_weights = ccc_results['optimization_results']['optimal_weights'] dcc_weights = dcc_results['optimization_results']['optimal_weights'] for i, ticker in enumerate(ccc_results['portfolio_summary']['tickers']): # Handle both dictionary and array formats if isinstance(ccc_weights, dict): ccc_weight = ccc_weights[ticker] else: ccc_weight = ccc_weights[i] if isinstance(dcc_weights, dict): dcc_weight = dcc_weights[ticker] else: dcc_weight = dcc_weights[i] weight_diff = dcc_weight - ccc_weight weight_change_pct = (weight_diff / ccc_weight) * 100 if ccc_weight > 0 else 0 total_turnover += abs(weight_diff) if abs(weight_diff) > 0.05: # 5% threshold for significant changes significant_changes.append((ticker, weight_diff, weight_change_pct)) print(f"{ticker:<8} {ccc_weight:<12.4f} {dcc_weight:<12.4f} {weight_diff:<+12.4f} {weight_change_pct:<+10.1f}%") print(f"\nTotal Portfolio Turnover: {total_turnover/2:.4f} ({(total_turnover/2)*100:.1f}%)") if significant_changes: print(f"Significant Changes (>5%):") for ticker, diff, pct in significant_changes: print(f"{ticker}: {diff:+.4f} ({pct:+.1f}%)") else: print("No significant allocation changes between models") # 3. Correlation Analysis print(f"\n3. CORRELATION STRUCTURE COMPARISON:") # CCC uses constant correlation ccc_corr = np.array(ccc_results['risk_metrics']['correlation_matrix']) avg_ccc_corr = np.mean(ccc_corr[np.triu_indices_from(ccc_corr, k=1)]) # DCC has time-varying correlations dcc_avg_corr = np.array(dcc_results['dcc_results']['average_correlations']) avg_dcc_corr = np.mean(dcc_avg_corr[np.triu_indices_from(dcc_avg_corr, k=1)]) dcc_final_corr = np.array(dcc_results['dcc_results']['final_correlation']) avg_final_corr = np.mean(dcc_final_corr[np.triu_indices_from(dcc_final_corr, k=1)]) print(f"CCC Average Correlation: {avg_ccc_corr:.4f}") print(f"DCC Time-Average Correlation: {avg_dcc_corr:.4f}") print(f"DCC Final Period Correlation: {avg_final_corr:.4f}") corr_variability = np.std([avg_dcc_corr, avg_final_corr]) print(f"DCC Correlation Variability: {corr_variability:.4f}") if corr_variability > 0.02: # 2% threshold for significant variability print("*** DCC shows significant time variation in correlations ***") else: print("*** Correlations are relatively stable over time ***") # 4. Model Complexity and Efficiency print(f"\n4. MODEL COMPLEXITY COMPARISON:") print(f"CCC Parameters: 0 (uses sample correlation)") print(f"DCC Parameters: 2 (alpha, beta)") print(f"DCC Alpha: {dcc_results['dcc_results']['dcc_parameters']['alpha']:.6f}") print(f"DCC Beta: {dcc_results['dcc_results']['dcc_parameters']['beta']:.6f}") print(f"DCC Persistence: {dcc_results['dcc_results']['dcc_parameters']['alpha'] + dcc_results['dcc_results']['dcc_parameters']['beta']:.6f}") # Interpretation of DCC parameters alpha = dcc_results['dcc_results']['dcc_parameters']['alpha'] beta = dcc_results['dcc_results']['dcc_parameters']['beta'] persistence = alpha + beta if alpha > 0.05: print("*** High sensitivity to recent shocks (alpha > 0.05) ***") elif alpha < 0.01: print("*** Low sensitivity to recent shocks (alpha < 0.01) ***") else: print("*** Moderate sensitivity to recent shocks ***") if persistence > 0.99: print("*** Very high correlation persistence (near unit root) ***") elif persistence > 0.95: print("*** High correlation persistence ***") else: print("*** Moderate correlation persistence ***") # 5. Statistical Model Selection Tests print(f"\n5. MODEL SELECTION CRITERIA:") # Simple AIC-like comparison (pseudo-AIC based on portfolio performance) # Lower volatility = better fit, but penalize DCC for extra parameters ccc_pseudo_aic = 2 * np.log(ccc_vol) # 0 parameters dcc_pseudo_aic = 2 * np.log(dcc_vol) + 2 * 2 # 2 parameters, penalty factor of 2 print(f"CCC Pseudo-AIC: {ccc_pseudo_aic:.6f}") print(f"DCC Pseudo-AIC: {dcc_pseudo_aic:.6f}") if dcc_pseudo_aic < ccc_pseudo_aic: aic_winner = "DCC" aic_improvement = ccc_pseudo_aic - dcc_pseudo_aic else: aic_winner = "CCC" aic_improvement = dcc_pseudo_aic - ccc_pseudo_aic print(f"AIC Winner: {aic_winner} (improvement: {aic_improvement:.6f})") # Likelihood ratio test approximation # If DCC provides significantly better fit, it justifies the extra complexity lr_statistic = 2 * (np.log(ccc_vol) - np.log(dcc_vol)) * len(ccc_results['portfolio_summary']['tickers']) * 1000 # Scaled for visibility lr_critical = 5.991 # Chi-square critical value for 2 df at 5% significance print(f"Likelihood Ratio Test:") print(f"LR Statistic: {lr_statistic:.3f}") print(f"Critical Value: {lr_critical:.3f} (5% significance)") if lr_statistic > lr_critical: print("*** DCC significantly outperforms CCC (reject H0: CCC adequate) ***") else: print("*** No significant evidence that DCC outperforms CCC ***") # 6. Practical Recommendations print(f"\n6. PRACTICAL RECOMMENDATIONS:") # Decision matrix dcc_advantages = 0 ccc_advantages = 0 # Risk reduction if vol_diff < -0.001: # DCC reduces risk by more than 0.1% dcc_advantages += 1 print("DCC: Lower portfolio risk") elif vol_diff > 0.001: # CCC has lower risk ccc_advantages += 1 print("CCC: Lower portfolio risk") else: print("Similar portfolio risk") # Correlation dynamics if corr_variability > 0.02: dcc_advantages += 1 print("DCC: Captures correlation dynamics") else: ccc_advantages += 1 print("CCC: Stable correlations sufficient") # Model complexity if lr_statistic > lr_critical: dcc_advantages += 1 print("DCC: Statistical significance justifies complexity") else: ccc_advantages += 1 print("CCC: Simpler model adequate") # Portfolio turnover if total_turnover/2 < 0.1: # Less than 10% turnover print("Similar portfolio allocations") elif len(significant_changes) > 0: print("Significant allocation changes - consider transaction costs") # Final recommendation print(f"\nFINAL RECOMMENDATION:") if dcc_advantages > ccc_advantages: print("*** USE DCC MODEL ***") print("Reasons: Better risk-return profile and captures market dynamics") elif ccc_advantages > dcc_advantages: print("*** USE CCC MODEL ***") print("Reasons: Simpler, more robust, and adequate for current data") else: print("*** EITHER MODEL ACCEPTABLE ***") print("Consider: CCC for simplicity, DCC for sophistication") # 7. Save Comparison Results comparison_results = { 'comparison_timestamp': datetime.now().isoformat(), 'portfolio_summary': ccc_results['portfolio_summary'], 'risk_comparison': { 'ccc_volatility': float(ccc_vol), 'dcc_volatility': float(dcc_vol), 'volatility_difference': float(vol_diff), 'relative_improvement_pct': float(vol_improvement) }, 'allocation_comparison': { 'total_turnover': float(total_turnover/2), 'significant_changes': significant_changes, 'weight_differences': [ { 'ticker': ticker, 'ccc_weight': float(ccc_weights[ticker] if isinstance(ccc_weights, dict) else ccc_weights[i]), 'dcc_weight': float(dcc_weights[ticker] if isinstance(dcc_weights, dict) else dcc_weights[i]), 'difference': float((dcc_weights[ticker] if isinstance(dcc_weights, dict) else dcc_weights[i]) - (ccc_weights[ticker] if isinstance(ccc_weights, dict) else ccc_weights[i])) } for i, ticker in enumerate(ccc_results['portfolio_summary']['tickers']) ] }, 'correlation_comparison': { 'ccc_avg_correlation': float(avg_ccc_corr), 'dcc_avg_correlation': float(avg_dcc_corr), 'dcc_final_correlation': float(avg_final_corr), 'correlation_variability': float(corr_variability) }, 'model_selection': { 'ccc_pseudo_aic': float(ccc_pseudo_aic), 'dcc_pseudo_aic': float(dcc_pseudo_aic), 'aic_winner': aic_winner, 'lr_statistic': float(lr_statistic), 'lr_critical_value': float(lr_critical), 'lr_test_significant': lr_statistic > lr_critical }, 'dcc_parameters': dcc_results['dcc_results']['dcc_parameters'], 'recommendation': { 'dcc_advantages': dcc_advantages, 'ccc_advantages': ccc_advantages, 'recommended_model': 'DCC' if dcc_advantages > ccc_advantages else 'CCC' if ccc_advantages > dcc_advantages else 'Either' } } # Save comparison results #comparison_filename = f"model_comparison_{datetime.now().strftime('%Y%m%d_%H%M%S')}.json" #with open(comparison_filename, 'w') as f: # json.dump(comparison_results, f, indent=2, default=str) #print(f"\n Detailed comparison results saved to: {comparison_filename}") print("\nANALYSIS COMPLETED SUCCESSFULLY") except Exception as e: print(f"Error in analysis: {str(e)}")