Value at Risk (VaR): From Theory to Production

1import numpy as np
2import pandas as pd
3from typing import Tuple, Optional
4from dataclasses import dataclass
5
6@dataclass
7class VaRResult:
8    """Container for VaR calculation results"""
9    var: float
10    confidence_level: float
11    horizon_days: int
12    methodology: str
13    calculation_time: float
14    num_scenarios: int
15
16class HistoricalSimulationVaR:
17    """
18    Production-grade Historical Simulation VaR calculator
19    
20    Features:
21    - Efficient numpy-based calculations
22    - Handles missing data
23    - Supports multiple confidence levels
24    - Provides detailed diagnostics
25    """
26    
27    def __init__(self, 
28                 returns_window: int = 250,
29                 min_observations: int = 100):
30        """
31        Args:
32            returns_window: Number of historical days to use
33            min_observations: Minimum required observations
34        """
35        self.returns_window = returns_window
36        self.min_observations = min_observations
37        
38    def calculate(self,
39                  positions: pd.DataFrame,
40                  historical_prices: pd.DataFrame,
41                  confidence_level: float = 0.99,
42                  horizon_days: int = 1) -> VaRResult:
43        """
44        Calculate VaR using historical simulation
45        
46        Args:
47            positions: DataFrame with columns ['asset', 'quantity', 'current_price']
48            historical_prices: DataFrame with assets as columns, dates as index
49            confidence_level: VaR confidence level (e.g., 0.99 for 99%)
50            horizon_days: VaR horizon in days
51            
52        Returns:
53            VaRResult object with VaR and diagnostics
54        """
55        import time
56        start_time = time.time()
57        
58        # Validate inputs
59        self._validate_inputs(positions, historical_prices)
60        
61        # Calculate historical returns
62        returns = historical_prices.pct_change().dropna()
63        
64        # Use only the most recent returns_window observations
65        returns = returns.tail(self.returns_window)
66        
67        if len(returns) < self.min_observations:
68            raise ValueError(f"Insufficient data: {len(returns)} < {self.min_observations}")
69        
70        # Scale returns to horizon
71        if horizon_days > 1:
72            returns = returns * np.sqrt(horizon_days)
73        
74        # Calculate current portfolio value
75        positions_dict = positions.set_index('asset')['quantity'].to_dict()
76        current_value = (positions['quantity'] * positions['current_price']).sum()
77        
78        # Calculate hypothetical P&L for each historical scenario
79        pnl_scenarios = np.zeros(len(returns))
80        
81        for asset in positions['asset']:
82            if asset not in returns.columns:
83                raise ValueError(f"Asset {asset} not found in historical prices")
84            
85            quantity = positions_dict[asset]
86            current_price = positions[positions['asset'] == asset]['current_price'].iloc[0]
87            
88            # Hypothetical P&L = quantity * current_price * return
89            pnl_scenarios += quantity * current_price * returns[asset].values
90        
91        # VaR is the negative of the alpha-quantile
92        var = -np.percentile(pnl_scenarios, (1 - confidence_level) * 100)
93        
94        calculation_time = time.time() - start_time
95        
96        return VaRResult(
97            var=var,
98            confidence_level=confidence_level,
99            horizon_days=horizon_days,
100            methodology="Historical Simulation",
101            calculation_time=calculation_time,
102            num_scenarios=len(returns)
103        )
104    
105    def _validate_inputs(self, 
106                        positions: pd.DataFrame, 
107                        historical_prices: pd.DataFrame):
108        """Validate input data"""
109        required_cols = ['asset', 'quantity', 'current_price']
110        if not all(col in positions.columns for col in required_cols):
111            raise ValueError(f"positions must have columns: {required_cols}")
112        
113        if historical_prices.empty:
114            raise ValueError("historical_prices cannot be empty")
115        
116        if not isinstance(historical_prices.index, pd.DatetimeIndex):
117            raise ValueError("historical_prices must have DatetimeIndex")
118
119
120# Example usage
121if __name__ == "__main__":
122    # Sample portfolio: Long 100 shares of AAPL, 200 shares of MSFT
123    positions = pd.DataFrame({
124        'asset': ['AAPL', 'MSFT'],
125        'quantity': [100, 200],
126        'current_price': [150.0, 300.0]
127    })
128    
129    # Generate sample historical prices (in production, load from database)
130    np.random.seed(42)
131    dates = pd.date_range(end='2025-11-25', periods=300, freq='D')
132    historical_prices = pd.DataFrame({
133        'AAPL': 150 * np.exp(np.random.randn(300).cumsum() * 0.02),
134        'MSFT': 300 * np.exp(np.random.randn(300).cumsum() * 0.02)
135    }, index=dates)
136    
137    # Calculate VaR
138    var_calculator = HistoricalSimulationVaR(returns_window=250)
139    result = var_calculator.calculate(
140        positions=positions,
141        historical_prices=historical_prices,
142        confidence_level=0.99,
143        horizon_days=1
144    )
145    
146    print(f"1-Day 99% VaR: ${result.var:,.2f}")
147    print(f"Calculation time: {result.calculation_time*1000:.2f}ms")
148    print(f"Number of scenarios: {result.num_scenarios}")
149

1import multiprocessing as mp
2from functools import partial
3
4def calculate_scenario_pnl(scenario_idx: int,
5                          returns: np.ndarray,
6                          positions: np.ndarray,
7                          prices: np.ndarray) -> float:
8    """Calculate P&L for a single scenario (for parallel processing)"""
9    return np.sum(positions * prices * returns[scenario_idx])
10
11class ParallelHistoricalVaR:
12    """Multi-core historical VaR calculation"""
13    
14    def __init__(self, num_processes: Optional[int] = None):
15        self.num_processes = num_processes or mp.cpu_count()
16    
17    def calculate(self, positions, historical_prices, confidence_level=0.99):
18        returns = historical_prices.pct_change().dropna().values
19        positions_array = positions['quantity'].values
20        prices_array = positions['current_price'].values
21        
22        # Parallel scenario calculation
23        with mp.Pool(self.num_processes) as pool:
24            calc_func = partial(
25                calculate_scenario_pnl,
26                returns=returns,
27                positions=positions_array,
28                prices=prices_array
29            )
30            pnl_scenarios = pool.map(calc_func, range(len(returns)))
31        
32        var = -np.percentile(pnl_scenarios, (1 - confidence_level) * 100)
33        return var
34

1from scipy import stats
2from scipy.linalg import cholesky
3
4class ParametricVaR:
5    """
6    Variance-Covariance (Parametric) VaR calculator
7    
8    Assumes multivariate normal returns
9    """
10    
11    def __init__(self, 
12                 returns_window: int = 250,
13                 use_exponential_weighting: bool = False,
14                 lambda_decay: float = 0.94):
15        """
16        Args:
17            returns_window: Historical window for covariance estimation
18            use_exponential_weighting: Use EWMA for covariance
19            lambda_decay: Decay factor for EWMA (RiskMetrics uses 0.94)
20        """
21        self.returns_window = returns_window
22        self.use_exponential_weighting = use_exponential_weighting
23        self.lambda_decay = lambda_decay
24    
25    def calculate(self,
26                  positions: pd.DataFrame,
27                  historical_prices: pd.DataFrame,
28                  confidence_level: float = 0.99,
29                  horizon_days: int = 1) -> VaRResult:
30        """Calculate parametric VaR"""
31        import time
32        start_time = time.time()
33        
34        # Calculate returns
35        returns = historical_prices.pct_change().dropna()
36        returns = returns.tail(self.returns_window)
37        
38        # Estimate covariance matrix
39        if self.use_exponential_weighting:
40            cov_matrix = self._ewma_covariance(returns)
41        else:
42            cov_matrix = returns.cov()
43        
44        # Calculate portfolio weights (in dollar terms)
45        positions_dict = positions.set_index('asset')
46        portfolio_value = (positions['quantity'] * positions['current_price']).sum()
47        
48        weights = {}
49        for asset in returns.columns:
50            if asset in positions_dict.index:
51                qty = positions_dict.loc[asset, 'quantity']
52                price = positions_dict.loc[asset, 'current_price']
53                weights[asset] = (qty * price) / portfolio_value
54            else:
55                weights[asset] = 0.0
56        
57        weights_array = np.array([weights[asset] for asset in returns.columns])
58        
59        # Portfolio variance: w^T * Sigma * w
60        portfolio_variance = weights_array @ cov_matrix.values @ weights_array
61        portfolio_std = np.sqrt(portfolio_variance)
62        
63        # Scale to horizon
64        if horizon_days > 1:
65            portfolio_std *= np.sqrt(horizon_days)
66        
67        # VaR = z_alpha * sigma * portfolio_value
68        z_alpha = stats.norm.ppf(confidence_level)
69        var = z_alpha * portfolio_std * portfolio_value
70        
71        calculation_time = time.time() - start_time
72        
73        return VaRResult(
74            var=var,
75            confidence_level=confidence_level,
76            horizon_days=horizon_days,
77            methodology="Parametric (Variance-Covariance)",
78            calculation_time=calculation_time,
79            num_scenarios=len(returns)
80        )
81    
82    def _ewma_covariance(self, returns: pd.DataFrame) -> pd.DataFrame:
83        """
84        Calculate exponentially weighted moving average covariance
85        (RiskMetrics methodology)
86        """
87        n_assets = len(returns.columns)
88        n_obs = len(returns)
89        
90        # Initialize covariance matrix with first observation
91        cov = np.outer(returns.iloc[0].values, returns.iloc[0].values)
92        
93        # EWMA recursion: Cov_t = lambda * Cov_{t-1} + (1-lambda) * r_t * r_t^T
94        for t in range(1, n_obs):
95            r_t = returns.iloc[t].values
96            cov = self.lambda_decay * cov + (1 - self.lambda_decay) * np.outer(r_t, r_t)
97        
98        return pd.DataFrame(cov, index=returns.columns, columns=returns.columns)
99    
100    def calculate_component_var(self,
101                               positions: pd.DataFrame,
102                               historical_prices: pd.DataFrame,
103                               confidence_level: float = 0.99) -> pd.DataFrame:
104        """
105        Calculate component VaR (contribution of each asset to total VaR)
106        
107        Component VaR_i = w_i * (dVaR/dw_i) = w_i * z_alpha * (Sigma * w)_i / sigma_p
108        """
109        returns = historical_prices.pct_change().dropna().tail(self.returns_window)
110        cov_matrix = returns.cov()
111        
112        portfolio_value = (positions['quantity'] * positions['current_price']).sum()
113        positions_dict = positions.set_index('asset')
114        
115        weights = np.array([
116            (positions_dict.loc[asset, 'quantity'] * 
117             positions_dict.loc[asset, 'current_price'] / portfolio_value)
118            if asset in positions_dict.index else 0.0
119            for asset in returns.columns
120        ])
121        
122        portfolio_variance = weights @ cov_matrix.values @ weights
123        portfolio_std = np.sqrt(portfolio_variance)
124        
125        z_alpha = stats.norm.ppf(confidence_level)
126        
127        # Marginal VaR: dVaR/dw_i
128        marginal_var = z_alpha * (cov_matrix.values @ weights) / portfolio_std
129        
130        # Component VaR: w_i * marginal_var_i
131        component_var = weights * marginal_var * portfolio_value
132        
133        results = pd.DataFrame({
134            'asset': returns.columns,
135            'weight': weights,
136            'marginal_var': marginal_var,
137            'component_var': component_var,
138            'pct_contribution': component_var / component_var.sum() * 100
139        })
140        
141        return results[results['weight'] > 0]
142

1from scipy.stats import multivariate_normal
2
3class MonteCarloVaR:
4    """
5    Monte Carlo VaR with variance reduction techniques
6    """
7    
8    def __init__(self,
9                 num_simulations: int = 10000,
10                 use_antithetic: bool = True,
11                 use_control_variate: bool = False,
12                 random_seed: Optional[int] = None):
13        """
14        Args:
15            num_simulations: Number of Monte Carlo paths
16            use_antithetic: Use antithetic variates for variance reduction
17            use_control_variate: Use control variates (requires analytical benchmark)
18            random_seed: Random seed for reproducibility
19        """
20        self.num_simulations = num_simulations
21        self.use_antithetic = use_antithetic
22        self.use_control_variate = use_control_variate
23        self.random_seed = random_seed
24        
25        if random_seed is not None:
26            np.random.seed(random_seed)
27    
28    def calculate(self,
29                  positions: pd.DataFrame,
30                  historical_prices: pd.DataFrame,
31                  confidence_level: float = 0.99,
32                  horizon_days: int = 1,
33                  model: str = 'gbm') -> VaRResult:
34        """
35        Calculate VaR using Monte Carlo simulation
36        
37        Args:
38            model: 'gbm' (Geometric Brownian Motion) or 'garch'
39        """
40        import time
41        start_time = time.time()
42        
43        returns = historical_prices.pct_change().dropna()
44        
45        # Estimate parameters
46        if model == 'gbm':
47            mu = returns.mean().values
48            cov = returns.cov().values
49        elif model == 'garch':
50            # Simplified: use GARCH(1,1) forecasted volatility
51            mu = returns.mean().values
52            cov = self._garch_forecast_cov(returns)
53        else:
54            raise ValueError(f"Unknown model: {model}")
55        
56        # Generate scenarios
57        scenarios = self._generate_scenarios(mu, cov, horizon_days)
58        
59        # Calculate P&L for each scenario
60        pnl_scenarios = self._calculate_pnl(scenarios, positions, historical_prices)
61        
62        # VaR is the negative of the alpha-quantile
63        var = -np.percentile(pnl_scenarios, (1 - confidence_level) * 100)
64        
65        calculation_time = time.time() - start_time
66        
67        return VaRResult(
68            var=var,
69            confidence_level=confidence_level,
70            horizon_days=horizon_days,
71            methodology=f"Monte Carlo ({model.upper()})",
72            calculation_time=calculation_time,
73            num_scenarios=self.num_simulations
74        )
75    
76    def _generate_scenarios(self, 
77                           mu: np.ndarray, 
78                           cov: np.ndarray,
79                           horizon_days: int) -> np.ndarray:
80        """Generate correlated random scenarios"""
81        n_assets = len(mu)
82        
83        if self.use_antithetic:
84            # Generate half the scenarios, create antithetic pairs
85            n_half = self.num_simulations // 2
86            Z = np.random.multivariate_normal(
87                mean=np.zeros(n_assets),
88                cov=cov,
89                size=n_half
90            )
91            # Antithetic variates: if Z ~ N(0, Sigma), then -Z ~ N(0, Sigma)
92            Z_anti = -Z
93            Z_combined = np.vstack([Z, Z_anti])
94            
95            # Add drift
96            scenarios = mu + Z_combined * np.sqrt(horizon_days)
97        else:
98            scenarios = np.random.multivariate_normal(
99                mean=mu * horizon_days,
100                cov=cov * horizon_days,
101                size=self.num_simulations
102            )
103        
104        return scenarios
105    
106    def _calculate_pnl(self,
107                      scenarios: np.ndarray,
108                      positions: pd.DataFrame,
109                      historical_prices: pd.DataFrame) -> np.ndarray:
110        """Calculate P&L for each scenario"""
111        positions_dict = positions.set_index('asset')
112        asset_list = list(historical_prices.columns)
113        
114        pnl_scenarios = np.zeros(len(scenarios))
115        
116        for i, asset in enumerate(asset_list):
117            if asset in positions_dict.index:
118                qty = positions_dict.loc[asset, 'quantity']
119                price = positions_dict.loc[asset, 'current_price']
120                
121                # P&L = quantity * price * return
122                pnl_scenarios += qty * price * scenarios[:, i]
123        
124        return pnl_scenarios
125    
126    def _garch_forecast_cov(self, returns: pd.DataFrame) -> np.ndarray:
127        """
128        Forecast covariance using GARCH(1,1)
129        Simplified implementation - in production, use arch library
130        """
131        # For simplicity, use exponentially weighted covariance
132        # In production, fit proper GARCH(1,1) model
133        lambda_decay = 0.94
134        cov = returns.ewm(alpha=1-lambda_decay).cov().iloc[-len(returns.columns):]
135        return cov.values
136

1from arch import arch_model
2
3class GARCHVaR:
4    """VaR calculation using GARCH(1,1) volatility forecasting"""
5    
6    def calculate(self,
7                  returns: pd.Series,
8                  confidence_level: float = 0.99,
9                  horizon_days: int = 1) -> float:
10        """
11        Calculate VaR using GARCH(1,1) volatility forecast
12        
13        Args:
14            returns: Historical returns (single asset or portfolio)
15            confidence_level: VaR confidence level
16            horizon_days: VaR horizon
17        """
18        # Fit GARCH(1,1) model
19        model = arch_model(returns * 100, vol='Garch', p=1, q=1)
20        fitted = model.fit(disp='off')
21        
22        # Forecast volatility for next horizon_days
23        forecast = fitted.forecast(horizon=horizon_days)
24        forecasted_variance = forecast.variance.values[-1, :]
25        
26        # Multi-period volatility (assuming i.i.d. innovations)
27        total_variance = forecasted_variance.sum()
28        forecasted_vol = np.sqrt(total_variance) / 100  # Convert back from percentage
29        
30        # VaR = z_alpha * sigma
31        z_alpha = stats.norm.ppf(confidence_level)
32        var = z_alpha * forecasted_vol
33        
34        return var
35

1from scipy.stats import chi2
2
3class VaRBacktester:
4    """
5    Comprehensive VaR backtesting framework
6    
7    Implements:
8    - Kupiec POF test
9    - Christoffersen conditional coverage test
10    - Traffic light approach (Basel)
11    """
12    
13    def __init__(self, var_series: pd.Series, returns_series: pd.Series):
14        """
15        Args:
16            var_series: Time series of VaR estimates
17            returns_series: Time series of actual returns (negative = loss)
18        """
19        if len(var_series) != len(returns_series):
20            raise ValueError("VaR and returns series must have same length")
21        
22        self.var_series = var_series
23        self.returns_series = returns_series
24        self.breaches = (returns_series < -var_series).astype(int)
25        self.num_breaches = self.breaches.sum()
26        self.num_observations = len(var_series)
27    
28    def kupiec_pof_test(self, confidence_level: float = 0.99) -> dict:
29        """
30        Kupiec Proportion of Failures test
31        
32        Returns:
33            dict with test results
34        """
35        alpha = 1 - confidence_level
36        x = self.num_breaches
37        N = self.num_observations
38        p_hat = x / N
39        
40        # Likelihood ratio test statistic
41        if x == 0:
42            lr_stat = -2 * N * np.log(1 - alpha)
43        elif x == N:
44            lr_stat = -2 * N * np.log(alpha)
45        else:
46            lr_stat = -2 * (
47                (N - x) * np.log((1 - alpha) / (1 - p_hat)) +
48                x * np.log(alpha / p_hat)
49            )
50        
51        # P-value from chi-squared distribution with 1 degree of freedom
52        p_value = 1 - chi2.cdf(lr_stat, df=1)
53        
54        # Decision: reject if p-value < 0.05
55        reject = p_value < 0.05
56        
57        return {
58            'test': 'Kupiec POF',
59            'num_breaches': x,
60            'num_observations': N,
61            'failure_rate': p_hat,
62            'expected_failure_rate': alpha,
63            'lr_statistic': lr_stat,
64            'p_value': p_value,
65            'reject_null': reject,
66            'conclusion': 'Model is REJECTED' if reject else 'Model is ACCEPTED'
67        }
68    
69    def christoffersen_test(self, confidence_level: float = 0.99) -> dict:
70        """
71        Christoffersen conditional coverage test
72        
73        Tests both:
74        1. Unconditional coverage (like Kupiec)
75        2. Independence of breaches (no clustering)
76        """
77        alpha = 1 - confidence_level
78        
79        # Unconditional coverage (same as Kupiec)
80        uc_result = self.kupiec_pof_test(confidence_level)
81        lr_uc = uc_result['lr_statistic']
82        
83        # Independence test
84        # Count transitions: 00, 01, 10, 11
85        breaches_array = self.breaches.values
86        n_00 = np.sum((breaches_array[:-1] == 0) & (breaches_array[1:] == 0))
87        n_01 = np.sum((breaches_array[:-1] == 0) & (breaches_array[1:] == 1))
88        n_10 = np.sum((breaches_array[:-1] == 1) & (breaches_array[1:] == 0))
89        n_11 = np.sum((breaches_array[:-1] == 1) & (breaches_array[1:] == 1))
90        
91        # Transition probabilities
92        pi_0 = n_01 / (n_00 + n_01) if (n_00 + n_01) > 0 else 0
93        pi_1 = n_11 / (n_10 + n_11) if (n_10 + n_11) > 0 else 0
94        pi = (n_01 + n_11) / (n_00 + n_01 + n_10 + n_11)
95        
96        # Independence LR statistic
97        if pi_0 == 0 or pi_1 == 0 or pi == 0:
98            lr_ind = 0
99        else:
100            lr_ind = -2 * (
101                n_00 * np.log(1 - pi) + n_01 * np.log(pi) +
102                n_10 * np.log(1 - pi) + n_11 * np.log(pi) -
103                n_00 * np.log(1 - pi_0) - n_01 * np.log(pi_0) -
104                n_10 * np.log(1 - pi_1) - n_11 * np.log(pi_1)
105            )
106        
107        # Conditional coverage = UC + IND
108        lr_cc = lr_uc + lr_ind
109        p_value_cc = 1 - chi2.cdf(lr_cc, df=2)
110        reject_cc = p_value_cc < 0.05
111        
112        return {
113            'test': 'Christoffersen Conditional Coverage',
114            'lr_uc': lr_uc,
115            'lr_ind': lr_ind,
116            'lr_cc': lr_cc,
117            'p_value': p_value_cc,
118            'reject_null': reject_cc,
119            'transition_probs': {'pi_0': pi_0, 'pi_1': pi_1},
120            'conclusion': 'Model is REJECTED' if reject_cc else 'Model is ACCEPTED'
121        }
122    
123    def basel_traffic_light(self, confidence_level: float = 0.99) -> dict:
124        """
125        Basel traffic light approach for VaR backtesting
126        
127        Green zone: 0-4 breaches (acceptable)
128        Yellow zone: 5-9 breaches (concerns)
129        Red zone: 10+ breaches (unacceptable)
130        
131        Assumes 250 trading days (1 year) and 99% VaR
132        """
133        if self.num_observations < 250:
134            raise ValueError("Traffic light test requires at least 250 observations")
135        
136        x = self.num_breaches
137        
138        if x <= 4:
139            zone = 'GREEN'
140            action = 'No action required'
141        elif x <= 9:
142            zone = 'YELLOW'
143            action = 'Increase capital multiplier'
144        else:
145            zone = 'RED'
146            action = 'Model must be revised'
147        
148        return {
149            'test': 'Basel Traffic Light',
150            'num_breaches': x,
151            'num_observations': self.num_observations,
152            'zone': zone,
153            'action': action
154        }
155    
156    def generate_report(self, confidence_level: float = 0.99) -> str:
157        """Generate comprehensive backtesting report"""
158        kupiec = self.kupiec_pof_test(confidence_level)
159        christoffersen = self.christoffersen_test(confidence_level)
160        
161        report = f"""
162VaR Backtesting Report
163{'=' * 60}
164
165Data Summary:
166- Number of observations: {self.num_observations}
167- Number of VaR breaches: {self.num_breaches}
168- Observed failure rate: {kupiec['failure_rate']:.4f}
169- Expected failure rate: {kupiec['expected_failure_rate']:.4f}
170
171Kupiec POF Test:
172- LR Statistic: {kupiec['lr_statistic']:.4f}
173- P-value: {kupiec['p_value']:.4f}
174- Conclusion: {kupiec['conclusion']}
175
176Christoffersen Conditional Coverage Test:
177- LR Statistic (CC): {christoffersen['lr_cc']:.4f}
178- P-value: {christoffersen['p_value']:.4f}
179- Conclusion: {christoffersen['conclusion']}
180- Independence: {'PASSED' if not christoffersen['reject_null'] else 'FAILED'}
181
182"""
183        
184        if self.num_observations >= 250:
185            traffic = self.basel_traffic_light(confidence_level)
186            report += f"""Basel Traffic Light:
187- Zone: {traffic['zone']}
188- Action: {traffic['action']}
189"""
190        
191        return report
192
193
194# Example: Backtest a VaR model
195if __name__ == "__main__":
196    # Generate sample data
197    np.random.seed(42)
198    dates = pd.date_range(end='2025-11-25', periods=300, freq='D')
199    
200    # Simulate returns with occasional large losses
201    returns = np.random.standard_t(df=5, size=300) * 0.02  # Fat tails
202    
203    # Calculate historical VaR
204    var_estimates = []
205    for i in range(250, 300):
206        historical_returns = returns[i-250:i]
207        var_99 = -np.percentile(historical_returns, 1)
208        var_estimates.append(var_99)
209    
210    var_series = pd.Series(var_estimates, index=dates[-50:])
211    returns_series = pd.Series(returns[-50:], index=dates[-50:])
212    
213    # Backtest
214    backtester = VaRBacktester(var_series, returns_series)
215    print(backtester.generate_report(confidence_level=0.99))
216

1#include <vector>
2#include <algorithm>
3#include <cmath>
4#include <numeric>
5#include <Eigen/Dense>
6
7using namespace Eigen;
8
9class HistoricalVaRCalculator {
10public:
11    struct VaRResult {
12        double var;
13        double confidence_level;
14        int horizon_days;
15        int num_scenarios;
16        double calculation_time_ms;
17    };
18    
19    HistoricalVaRCalculator(int returns_window = 250)
20        : returns_window_(returns_window) {}
21    
22    VaRResult calculate(
23        const VectorXd& positions,
24        const VectorXd& current_prices,
25        const MatrixXd& historical_prices,
26        double confidence_level = 0.99,
27        int horizon_days = 1
28    ) {
29        auto start = std::chrono::high_resolution_clock::now();
30        
31        // Calculate returns matrix
32        MatrixXd returns = calculateReturns(historical_prices);
33        
34        // Use only last returns_window observations
35        int start_idx = std::max(0, (int)returns.rows() - returns_window_);
36        MatrixXd recent_returns = returns.bottomRows(returns.rows() - start_idx);
37        
38        // Scale to horizon
39        if (horizon_days > 1) {
40            recent_returns *= std::sqrt(horizon_days);
41        }
42        
43        // Calculate P&L scenarios
44        VectorXd pnl_scenarios = calculatePnLScenarios(
45            positions, current_prices, recent_returns
46        );
47        
48        // Calculate VaR as negative of alpha-quantile
49        double var = -percentile(pnl_scenarios, (1.0 - confidence_level) * 100.0);
50        
51        auto end = std::chrono::high_resolution_clock::now();
52        double calc_time = std::chrono::duration<double, std::milli>(end - start).count();
53        
54        return VaRResult{
55            var,
56            confidence_level,
57            horizon_days,
58            (int)recent_returns.rows(),
59            calc_time
60        };
61    }
62    
63private:
64    int returns_window_;
65    
66    MatrixXd calculateReturns(const MatrixXd& prices) {
67        MatrixXd returns(prices.rows() - 1, prices.cols());
68        
69        for (int i = 0; i < prices.rows() - 1; ++i) {
70            for (int j = 0; j < prices.cols(); ++j) {
71                returns(i, j) = (prices(i + 1, j) - prices(i, j)) / prices(i, j);
72            }
73        }
74        
75        return returns;
76    }
77    
78    VectorXd calculatePnLScenarios(
79        const VectorXd& positions,
80        const VectorXd& current_prices,
81        const MatrixXd& returns
82    ) {
83        // P&L for each scenario = sum over assets of (position * price * return)
84        VectorXd position_values = positions.array() * current_prices.array();
85        VectorXd pnl_scenarios = returns * position_values;
86        
87        return pnl_scenarios;
88    }
89    
90    double percentile(const VectorXd& data, double pct) {
91        std::vector<double> sorted(data.data(), data.data() + data.size());
92        std::sort(sorted.begin(), sorted.end());
93        
94        double idx = pct / 100.0 * (sorted.size() - 1);
95        int lower = (int)std::floor(idx);
96        int upper = (int)std::ceil(idx);
97        double weight = idx - lower;
98        
99        return sorted[lower] * (1.0 - weight) + sorted[upper] * weight;
100    }
101};
102
103// Example usage
104int main() {
105    // Portfolio: 100 shares AAPL @ $150, 200 shares MSFT @ $300
106    VectorXd positions(2);
107    positions << 100, 200;
108    
109    VectorXd current_prices(2);
110    current_prices << 150.0, 300.0;
111    
112    // Historical prices (300 days x 2 assets)
113    MatrixXd historical_prices = MatrixXd::Random(300, 2) * 50 + 200;
114    
115    HistoricalVaRCalculator calculator(250);
116    auto result = calculator.calculate(
117        positions, current_prices, historical_prices, 0.99, 1
118    );
119    
120    std::cout << "1-Day 99% VaR: $" << result.var << std::endl;
121    std::cout << "Calculation time: " << result.calculation_time_ms << "ms" << std::endl;
122    
123    return 0;
124}
125

1def scale_var_to_10_day(var_1_day: float, method: str = 'sqrt') -> float:
2    """Scale 1-day VaR to 10-day VaR"""
3    if method == 'sqrt':
4        return var_1_day * np.sqrt(10)
5    elif method == 'conservative':
6        # More conservative: assume no diversification benefit
7        return var_1_day * 10
8    else:
9        raise ValueError(f"Unknown method: {method}")
10

1def calculate_incremental_var(
2    base_var: float,
3    positions: pd.DataFrame,
4    historical_prices: pd.DataFrame,
5    new_position: dict,  # {'asset': 'AAPL', 'quantity': 100}
6    var_calculator: HistoricalSimulationVaR
7) -> float:
8    """Calculate incremental VaR of adding a new position"""
9    
10    # Add new position to portfolio
11    new_positions = positions.copy()
12    new_row = pd.DataFrame([new_position])
13    new_positions = pd.concat([new_positions, new_row], ignore_index=True)
14    
15    # Calculate VaR with new position
16    new_var_result = var_calculator.calculate(
17        positions=new_positions,
18        historical_prices=historical_prices
19    )
20    
21    incremental_var = new_var_result.var - base_var
22    
23    return incremental_var
24

1def calculate_cvar(pnl_scenarios: np.ndarray, confidence_level: float = 0.99) -> float:
2    """Calculate Conditional VaR (Expected Shortfall)"""
3    var = -np.percentile(pnl_scenarios, (1 - confidence_level) * 100)
4    
5    # CVaR = average of losses beyond VaR
6    tail_losses = pnl_scenarios[pnl_scenarios <= -var]
7    cvar = -tail_losses.mean()
8    
9    return cvar
10