rustkernel_risk/
market.rs

1//! Market risk kernels.
2//!
3//! This module provides market risk analytics:
4//! - Monte Carlo Value at Risk (VaR)
5//! - Portfolio risk aggregation
6//! - Expected Shortfall (CVaR)
7
8use crate::messages::{
9    MonteCarloVaRInput, MonteCarloVaROutput, PortfolioRiskAggregationInput,
10    PortfolioRiskAggregationOutput,
11};
12use crate::ring_messages::{
13    K2KMarketUpdate, K2KMarketUpdateAck, K2KVaRAggregation, K2KVaRAggregationResponse,
14    QueryVaRResponse, QueryVaRRing, RecalculateVaRResponse, RecalculateVaRRing,
15    UpdatePositionResponse, UpdatePositionRing, from_currency_fp, from_fixed_point, to_currency_fp,
16};
17use crate::types::{Portfolio, PortfolioRiskResult, VaRParams, VaRResult};
18use async_trait::async_trait;
19use ringkernel_core::RingContext;
20use rustkernel_core::error::Result;
21use rustkernel_core::traits::{BatchKernel, RingKernelHandler};
22use rustkernel_core::{domain::Domain, kernel::KernelMetadata, traits::GpuKernel};
23use std::time::Instant;
24
25// ============================================================================
26// Monte Carlo VaR Kernel
27// ============================================================================
28
29/// Monte Carlo VaR state for Ring mode operations.
30#[derive(Debug, Clone, Default)]
31pub struct MonteCarloVaRState {
32    /// Current portfolio.
33    pub portfolio: Option<Portfolio>,
34    /// Cached VaR result.
35    pub var: f64,
36    /// Cached Expected Shortfall.
37    pub es: f64,
38    /// Confidence level used for cached result.
39    pub confidence_level: f64,
40    /// Holding period used for cached result.
41    pub holding_period: u32,
42    /// Whether cached result is stale.
43    pub is_stale: bool,
44    /// Number of simulations used.
45    pub n_simulations: u32,
46}
47
48/// Monte Carlo Value at Risk kernel.
49///
50/// Simulates portfolio P&L using correlated random variates.
51#[derive(Debug)]
52pub struct MonteCarloVaR {
53    metadata: KernelMetadata,
54    /// Internal state for Ring mode operations.
55    state: std::sync::RwLock<MonteCarloVaRState>,
56}
57
58impl Clone for MonteCarloVaR {
59    fn clone(&self) -> Self {
60        Self {
61            metadata: self.metadata.clone(),
62            state: std::sync::RwLock::new(self.state.read().unwrap().clone()),
63        }
64    }
65}
66
67impl Default for MonteCarloVaR {
68    fn default() -> Self {
69        Self::new()
70    }
71}
72
73impl MonteCarloVaR {
74    /// Create a new Monte Carlo VaR kernel.
75    #[must_use]
76    pub fn new() -> Self {
77        Self {
78            metadata: KernelMetadata::ring("risk/monte-carlo-var", Domain::RiskAnalytics)
79                .with_description("Monte Carlo Value at Risk simulation")
80                .with_throughput(100_000)
81                .with_latency_us(1000.0),
82            state: std::sync::RwLock::new(MonteCarloVaRState::default()),
83        }
84    }
85
86    /// Initialize the kernel with a portfolio for Ring mode operations.
87    pub fn initialize(&self, portfolio: Portfolio) {
88        let mut state = self.state.write().unwrap();
89        state.portfolio = Some(portfolio);
90        state.is_stale = true;
91    }
92
93    /// Update a position value in the portfolio.
94    pub fn update_position(&self, asset_id: u64, new_value: f64) -> bool {
95        let mut state = self.state.write().unwrap();
96        if let Some(ref mut portfolio) = state.portfolio {
97            // Find and update the asset
98            for (i, &id) in portfolio.asset_ids.iter().enumerate() {
99                if id == asset_id {
100                    portfolio.values[i] = new_value;
101                    state.is_stale = true;
102                    return true;
103                }
104            }
105        }
106        false
107    }
108
109    /// Get cached VaR result.
110    pub fn cached_var(&self) -> (f64, f64, bool) {
111        let state = self.state.read().unwrap();
112        (state.var, state.es, !state.is_stale)
113    }
114
115    /// Recalculate VaR using given parameters.
116    pub fn recalculate(
117        &self,
118        confidence_level: f64,
119        holding_period: u32,
120        n_simulations: u32,
121    ) -> (f64, f64) {
122        let mut state = self.state.write().unwrap();
123        let Some(ref portfolio) = state.portfolio else {
124            return (0.0, 0.0);
125        };
126
127        let params = VaRParams::new(confidence_level, holding_period, n_simulations);
128        let result = Self::compute(portfolio, params);
129
130        state.var = result.var;
131        state.es = result.expected_shortfall;
132        state.confidence_level = confidence_level;
133        state.holding_period = holding_period;
134        state.n_simulations = n_simulations;
135        state.is_stale = false;
136
137        (result.var, result.expected_shortfall)
138    }
139
140    /// Calculate VaR using Monte Carlo simulation.
141    ///
142    /// # Arguments
143    /// * `portfolio` - Portfolio to analyze
144    /// * `params` - VaR calculation parameters
145    pub fn compute(portfolio: &Portfolio, params: VaRParams) -> VaRResult {
146        if portfolio.n_assets() == 0 {
147            return VaRResult {
148                var: 0.0,
149                expected_shortfall: 0.0,
150                confidence_level: params.confidence_level,
151                holding_period: params.holding_period,
152                component_var: Vec::new(),
153                marginal_var: Vec::new(),
154                percentiles: Vec::new(),
155            };
156        }
157
158        let n_sims = params.n_simulations as usize;
159        let holding_factor = (params.holding_period as f64).sqrt();
160
161        // Generate correlated random returns using Cholesky decomposition
162        let cholesky = Self::cholesky_decomposition(portfolio);
163
164        // Simulate P&L scenarios
165        let mut pnl_scenarios = Vec::with_capacity(n_sims);
166        let mut rng = SimpleRng::new(42); // Deterministic for reproducibility
167
168        for _ in 0..n_sims {
169            // Generate independent standard normal variates
170            let z: Vec<f64> = (0..portfolio.n_assets()).map(|_| rng.normal()).collect();
171
172            // Correlate using Cholesky
173            let correlated = Self::apply_cholesky(&cholesky, &z, portfolio.n_assets());
174
175            // Calculate portfolio P&L
176            let mut scenario_pnl = 0.0;
177            for (i, (&z_corr, (&vol, &value))) in correlated
178                .iter()
179                .zip(portfolio.volatilities.iter().zip(portfolio.values.iter()))
180                .enumerate()
181            {
182                let ret = portfolio.expected_returns[i] * params.holding_period as f64 / 252.0
183                    + vol * holding_factor / (252.0_f64).sqrt() * z_corr;
184                scenario_pnl += value * ret;
185            }
186
187            pnl_scenarios.push(scenario_pnl);
188        }
189
190        // Sort scenarios for percentile calculation
191        pnl_scenarios.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
192
193        // Calculate VaR (loss at confidence level)
194        let var_idx = ((1.0 - params.confidence_level) * n_sims as f64) as usize;
195        let var = -pnl_scenarios[var_idx.min(n_sims - 1)]; // Positive VaR = loss
196
197        // Calculate Expected Shortfall (average of tail losses)
198        let tail_start = var_idx.min(n_sims - 1);
199        let expected_shortfall = if tail_start > 0 {
200            -pnl_scenarios[..tail_start].iter().sum::<f64>() / tail_start as f64
201        } else {
202            var
203        };
204
205        // Calculate component VaR (marginal contribution)
206        let component_var = Self::calculate_component_var(portfolio, var, &cholesky);
207
208        // Calculate marginal VaR
209        let marginal_var = Self::calculate_marginal_var(portfolio, params, &cholesky);
210
211        // Calculate percentiles
212        let percentiles: Vec<(f64, f64)> = [0.01, 0.05, 0.10, 0.25, 0.50, 0.75, 0.90, 0.95, 0.99]
213            .iter()
214            .map(|&p| {
215                let idx = ((1.0 - p) * n_sims as f64) as usize;
216                (p, pnl_scenarios[idx.min(n_sims - 1)])
217            })
218            .collect();
219
220        VaRResult {
221            var,
222            expected_shortfall,
223            confidence_level: params.confidence_level,
224            holding_period: params.holding_period,
225            component_var,
226            marginal_var,
227            percentiles,
228        }
229    }
230
231    /// Perform Cholesky decomposition of correlation matrix.
232    fn cholesky_decomposition(portfolio: &Portfolio) -> Vec<f64> {
233        let n = portfolio.n_assets();
234        let mut l = vec![0.0; n * n];
235
236        for i in 0..n {
237            for j in 0..=i {
238                let mut sum = 0.0;
239                for k in 0..j {
240                    sum += l[i * n + k] * l[j * n + k];
241                }
242
243                if i == j {
244                    let diag = portfolio.correlation(i, i) - sum;
245                    l[i * n + j] = if diag > 0.0 { diag.sqrt() } else { 0.0 };
246                } else {
247                    let l_jj = l[j * n + j];
248                    l[i * n + j] = if l_jj.abs() > 1e-10 {
249                        (portfolio.correlation(i, j) - sum) / l_jj
250                    } else {
251                        0.0
252                    };
253                }
254            }
255        }
256
257        l
258    }
259
260    /// Apply Cholesky factor to independent normals.
261    fn apply_cholesky(l: &[f64], z: &[f64], n: usize) -> Vec<f64> {
262        let mut result = vec![0.0; n];
263        for i in 0..n {
264            for j in 0..=i {
265                result[i] += l[i * n + j] * z[j];
266            }
267        }
268        result
269    }
270
271    /// Calculate component VaR.
272    #[allow(clippy::needless_range_loop)]
273    fn calculate_component_var(
274        portfolio: &Portfolio,
275        total_var: f64,
276        _cholesky: &[f64],
277    ) -> Vec<f64> {
278        // Component VaR = weight_i * sigma_i * rho_i,p * VaR / sigma_p
279        let weights = portfolio.weights();
280        let n = portfolio.n_assets();
281
282        // Calculate portfolio volatility
283        let port_var = Self::portfolio_variance(portfolio);
284        let port_vol = port_var.sqrt();
285
286        if port_vol < 1e-10 {
287            return vec![0.0; n];
288        }
289
290        // Calculate covariance of each asset with portfolio
291        let mut component_vars = Vec::with_capacity(n);
292
293        for i in 0..n {
294            let mut cov_i_p = 0.0;
295            for j in 0..n {
296                cov_i_p += weights[j]
297                    * portfolio.volatilities[i]
298                    * portfolio.volatilities[j]
299                    * portfolio.correlation(i, j);
300            }
301
302            let beta_i = cov_i_p / port_var;
303            let component_var_i = weights[i] * beta_i * total_var;
304            component_vars.push(component_var_i);
305        }
306
307        component_vars
308    }
309
310    /// Calculate marginal VaR.
311    #[allow(clippy::needless_range_loop)]
312    fn calculate_marginal_var(
313        portfolio: &Portfolio,
314        params: VaRParams,
315        _cholesky: &[f64],
316    ) -> Vec<f64> {
317        let n = portfolio.n_assets();
318        let weights = portfolio.weights();
319        let holding_factor = (params.holding_period as f64).sqrt();
320
321        // Z-score for confidence level
322        let z = Self::norm_inv(params.confidence_level);
323
324        // Portfolio volatility
325        let port_vol = Self::portfolio_variance(portfolio).sqrt();
326
327        if port_vol < 1e-10 {
328            return vec![0.0; n];
329        }
330
331        // Marginal VaR = dVaR/dw_i = z * sigma_p * d_sigma_p/dw_i
332        let mut marginal_vars = Vec::with_capacity(n);
333
334        for i in 0..n {
335            let mut cov_i_p = 0.0;
336            for j in 0..n {
337                cov_i_p += weights[j]
338                    * portfolio.volatilities[i]
339                    * portfolio.volatilities[j]
340                    * portfolio.correlation(i, j);
341            }
342
343            let d_sigma_dw = cov_i_p / port_vol;
344            let marginal_var_i = z * d_sigma_dw * holding_factor / (252.0_f64).sqrt();
345            marginal_vars.push(marginal_var_i * portfolio.total_value());
346        }
347
348        marginal_vars
349    }
350
351    /// Calculate portfolio variance.
352    fn portfolio_variance(portfolio: &Portfolio) -> f64 {
353        let weights = portfolio.weights();
354        let n = portfolio.n_assets();
355        let mut var = 0.0;
356
357        for i in 0..n {
358            for j in 0..n {
359                var += weights[i]
360                    * weights[j]
361                    * portfolio.volatilities[i]
362                    * portfolio.volatilities[j]
363                    * portfolio.correlation(i, j);
364            }
365        }
366
367        var
368    }
369
370    /// Standard normal inverse CDF.
371    fn norm_inv(p: f64) -> f64 {
372        // Rational approximation
373        let p_clamped = p.clamp(1e-10, 1.0 - 1e-10);
374        let t = if p_clamped < 0.5 {
375            (-2.0 * p_clamped.ln()).sqrt()
376        } else {
377            (-2.0 * (1.0 - p_clamped).ln()).sqrt()
378        };
379
380        let c0 = 2.515517;
381        let c1 = 0.802853;
382        let c2 = 0.010328;
383        let d1 = 1.432788;
384        let d2 = 0.189269;
385        let d3 = 0.001308;
386
387        let result = t - (c0 + c1 * t + c2 * t * t) / (1.0 + d1 * t + d2 * t * t + d3 * t * t * t);
388
389        if p_clamped < 0.5 { -result } else { result }
390    }
391}
392
393impl GpuKernel for MonteCarloVaR {
394    fn metadata(&self) -> &KernelMetadata {
395        &self.metadata
396    }
397}
398
399// ============================================================================
400// Monte Carlo VaR RingKernelHandler Implementations
401// ============================================================================
402
403/// RingKernelHandler for position updates.
404///
405/// Enables streaming position updates for real-time VaR.
406#[async_trait]
407impl RingKernelHandler<UpdatePositionRing, UpdatePositionResponse> for MonteCarloVaR {
408    async fn handle(
409        &self,
410        _ctx: &mut RingContext,
411        msg: UpdatePositionRing,
412    ) -> Result<UpdatePositionResponse> {
413        // Update position in internal state
414        let new_value = from_currency_fp(msg.value_fp);
415        let updated = self.update_position(msg.asset_id, new_value);
416
417        // Position changed, VaR needs recalculation
418        Ok(UpdatePositionResponse {
419            request_id: msg.id.0,
420            asset_id: msg.asset_id,
421            var_stale: updated,
422        })
423    }
424}
425
426/// RingKernelHandler for VaR queries.
427///
428/// Returns current cached VaR or triggers recalculation.
429#[async_trait]
430impl RingKernelHandler<QueryVaRRing, QueryVaRResponse> for MonteCarloVaR {
431    async fn handle(&self, _ctx: &mut RingContext, msg: QueryVaRRing) -> Result<QueryVaRResponse> {
432        // Query cached VaR from internal state
433        let (var, es, is_fresh) = self.cached_var();
434
435        Ok(QueryVaRResponse {
436            request_id: msg.id.0,
437            var_fp: to_currency_fp(var),
438            es_fp: to_currency_fp(es),
439            confidence_fp: msg.confidence_fp,
440            holding_period: msg.holding_period,
441            is_fresh,
442        })
443    }
444}
445
446/// RingKernelHandler for VaR recalculation.
447///
448/// Triggers full Monte Carlo simulation with given parameters.
449#[async_trait]
450impl RingKernelHandler<RecalculateVaRRing, RecalculateVaRResponse> for MonteCarloVaR {
451    async fn handle(
452        &self,
453        _ctx: &mut RingContext,
454        msg: RecalculateVaRRing,
455    ) -> Result<RecalculateVaRResponse> {
456        let start = Instant::now();
457        let confidence = from_fixed_point(msg.confidence_fp);
458
459        // Run full Monte Carlo simulation on internal state
460        let (var, es) = self.recalculate(confidence, msg.holding_period, msg.n_simulations);
461        let compute_time_us = start.elapsed().as_micros() as u64;
462
463        Ok(RecalculateVaRResponse {
464            request_id: msg.id.0,
465            var_fp: to_currency_fp(var),
466            es_fp: to_currency_fp(es),
467            compute_time_us,
468            n_simulations: msg.n_simulations,
469        })
470    }
471}
472
473/// RingKernelHandler for K2K market updates.
474///
475/// Processes streaming market data for real-time VaR adjustments.
476#[async_trait]
477impl RingKernelHandler<K2KMarketUpdate, K2KMarketUpdateAck> for MonteCarloVaR {
478    async fn handle(
479        &self,
480        _ctx: &mut RingContext,
481        msg: K2KMarketUpdate,
482    ) -> Result<K2KMarketUpdateAck> {
483        // Get current portfolio value for impact estimation
484        let portfolio_value = {
485            let state = self.state.read().unwrap();
486            state
487                .portfolio
488                .as_ref()
489                .map(|p| p.total_value())
490                .unwrap_or(0.0)
491        };
492
493        // Estimate VaR impact from volatility change
494        // VaR impact ≈ portfolio_value * z * Δvolatility * sqrt(holding_period)
495        let vol_delta = from_fixed_point(msg.vol_delta_fp);
496        let z_95 = 1.645; // 95% confidence z-score
497        let var_impact = portfolio_value * z_95 * vol_delta.abs();
498
499        // Mark state as stale since market conditions changed
500        {
501            let mut state = self.state.write().unwrap();
502            state.is_stale = true;
503        }
504
505        Ok(K2KMarketUpdateAck {
506            request_id: msg.id.0,
507            worker_id: 0, // Single instance worker
508            var_impact_fp: to_currency_fp(var_impact),
509        })
510    }
511}
512
513/// RingKernelHandler for K2K partial VaR aggregation.
514///
515/// Used in distributed VaR calculation to aggregate partial results.
516#[async_trait]
517impl RingKernelHandler<K2KVaRAggregation, K2KVaRAggregationResponse> for MonteCarloVaR {
518    async fn handle(
519        &self,
520        _ctx: &mut RingContext,
521        msg: K2KVaRAggregation,
522    ) -> Result<K2KVaRAggregationResponse> {
523        let complete = msg.workers_reported >= msg.expected_workers;
524        let aggregated_var = from_currency_fp(msg.aggregated_var_fp);
525
526        // Calculate diversification benefit based on portfolio correlation structure
527        // With n workers, diversification benefit scales with 1 - sqrt(1/n) for uncorrelated
528        // For partially correlated, use average correlation from portfolio if available
529        let diversification_factor = {
530            let state = self.state.read().unwrap();
531            if let Some(ref portfolio) = state.portfolio {
532                let n = portfolio.n_assets();
533                if n > 1 {
534                    // Calculate average off-diagonal correlation
535                    let mut avg_corr = 0.0;
536                    let mut count = 0;
537                    for i in 0..n {
538                        for j in (i + 1)..n {
539                            avg_corr += portfolio.correlation(i, j);
540                            count += 1;
541                        }
542                    }
543                    if count > 0 {
544                        avg_corr /= count as f64;
545                    }
546                    // Diversification benefit = 1 - sqrt(average correlation adjustment)
547                    // Higher correlation = less benefit
548                    (1.0 - avg_corr.max(0.0).sqrt()) * 0.3 // Cap at 30%
549                } else {
550                    0.0
551                }
552            } else {
553                0.15 // Default 15% when no portfolio available
554            }
555        };
556
557        let diversification_benefit = aggregated_var * diversification_factor;
558        let final_var = aggregated_var - diversification_benefit;
559
560        // ES ratio depends on confidence level (1.25x for 95%, 1.15x for 99%)
561        let es_ratio = 1.25;
562        let final_es = final_var * es_ratio;
563
564        Ok(K2KVaRAggregationResponse {
565            correlation_id: msg.correlation_id,
566            complete,
567            final_var_fp: to_currency_fp(final_var),
568            final_es_fp: to_currency_fp(final_es),
569            diversification_benefit_fp: to_currency_fp(diversification_benefit),
570        })
571    }
572}
573
574#[async_trait]
575impl BatchKernel<MonteCarloVaRInput, MonteCarloVaROutput> for MonteCarloVaR {
576    async fn execute(&self, input: MonteCarloVaRInput) -> Result<MonteCarloVaROutput> {
577        let start = Instant::now();
578        let result = Self::compute(&input.portfolio, input.params);
579        Ok(MonteCarloVaROutput {
580            result,
581            compute_time_us: start.elapsed().as_micros() as u64,
582        })
583    }
584}
585
586// ============================================================================
587// Portfolio Risk Aggregation Kernel
588// ============================================================================
589
590/// Portfolio risk aggregation kernel.
591///
592/// Aggregates risk measures across portfolio with correlation adjustment.
593#[derive(Debug, Clone)]
594pub struct PortfolioRiskAggregation {
595    metadata: KernelMetadata,
596}
597
598impl Default for PortfolioRiskAggregation {
599    fn default() -> Self {
600        Self::new()
601    }
602}
603
604impl PortfolioRiskAggregation {
605    /// Create a new portfolio risk aggregation kernel.
606    #[must_use]
607    pub fn new() -> Self {
608        Self {
609            metadata: KernelMetadata::ring("risk/portfolio-aggregation", Domain::RiskAnalytics)
610                .with_description("Correlation-adjusted portfolio risk")
611                .with_throughput(10_000)
612                .with_latency_us(500.0),
613        }
614    }
615
616    /// Aggregate portfolio risk.
617    ///
618    /// # Arguments
619    /// * `portfolio` - Portfolio to analyze
620    /// * `confidence_level` - Confidence level for VaR
621    /// * `holding_period` - Holding period in days
622    pub fn compute(
623        portfolio: &Portfolio,
624        confidence_level: f64,
625        holding_period: u32,
626    ) -> PortfolioRiskResult {
627        let n = portfolio.n_assets();
628        if n == 0 {
629            return PortfolioRiskResult {
630                portfolio_var: 0.0,
631                portfolio_es: 0.0,
632                undiversified_var: 0.0,
633                diversification_benefit: 0.0,
634                asset_vars: Vec::new(),
635                risk_contributions: Vec::new(),
636                covariance_matrix: Vec::new(),
637            };
638        }
639
640        let z = Self::norm_inv(confidence_level);
641        let holding_factor = (holding_period as f64).sqrt() / (252.0_f64).sqrt();
642
643        // Calculate covariance matrix
644        let mut cov_matrix = vec![0.0; n * n];
645        for i in 0..n {
646            for j in 0..n {
647                cov_matrix[i * n + j] = portfolio.volatilities[i]
648                    * portfolio.volatilities[j]
649                    * portfolio.correlation(i, j);
650            }
651        }
652
653        // Calculate individual asset VaRs
654        let asset_vars: Vec<f64> = portfolio
655            .values
656            .iter()
657            .zip(portfolio.volatilities.iter())
658            .map(|(&value, &vol)| value * z * vol * holding_factor)
659            .collect();
660
661        // Undiversified VaR (sum of individual VaRs)
662        let undiversified_var: f64 = asset_vars.iter().sum();
663
664        // Portfolio variance
665        let weights = portfolio.weights();
666        let mut portfolio_variance = 0.0;
667        for i in 0..n {
668            for j in 0..n {
669                portfolio_variance += weights[i] * weights[j] * cov_matrix[i * n + j];
670            }
671        }
672
673        // Portfolio VaR
674        let portfolio_var =
675            portfolio.total_value() * z * portfolio_variance.sqrt() * holding_factor;
676
677        // Diversification benefit
678        let diversification_benefit = undiversified_var - portfolio_var;
679
680        // Risk contributions (Euler allocation)
681        let port_vol = portfolio_variance.sqrt();
682        let risk_contributions: Vec<f64> = if port_vol > 1e-10 {
683            (0..n)
684                .map(|i| {
685                    let mut cov_i_p = 0.0;
686                    for j in 0..n {
687                        cov_i_p += weights[j] * cov_matrix[i * n + j];
688                    }
689                    weights[i] * cov_i_p / port_vol * z * holding_factor * portfolio.total_value()
690                })
691                .collect()
692        } else {
693            vec![0.0; n]
694        };
695
696        // Expected Shortfall (using normal approximation)
697        let pdf_at_z = (-z * z / 2.0).exp() / (2.0 * std::f64::consts::PI).sqrt();
698        let portfolio_es = portfolio.total_value() * port_vol * holding_factor * pdf_at_z
699            / (1.0 - confidence_level);
700
701        PortfolioRiskResult {
702            portfolio_var,
703            portfolio_es,
704            undiversified_var,
705            diversification_benefit,
706            asset_vars,
707            risk_contributions,
708            covariance_matrix: cov_matrix,
709        }
710    }
711
712    /// Standard normal inverse CDF.
713    fn norm_inv(p: f64) -> f64 {
714        MonteCarloVaR::norm_inv(p)
715    }
716}
717
718impl GpuKernel for PortfolioRiskAggregation {
719    fn metadata(&self) -> &KernelMetadata {
720        &self.metadata
721    }
722}
723
724#[async_trait]
725impl BatchKernel<PortfolioRiskAggregationInput, PortfolioRiskAggregationOutput>
726    for PortfolioRiskAggregation
727{
728    async fn execute(
729        &self,
730        input: PortfolioRiskAggregationInput,
731    ) -> Result<PortfolioRiskAggregationOutput> {
732        let start = Instant::now();
733        let result = Self::compute(
734            &input.portfolio,
735            input.confidence_level,
736            input.holding_period,
737        );
738        Ok(PortfolioRiskAggregationOutput {
739            result,
740            compute_time_us: start.elapsed().as_micros() as u64,
741        })
742    }
743}
744
745// ============================================================================
746// Simple RNG for Monte Carlo
747// ============================================================================
748
749/// Simple pseudo-random number generator (xorshift64).
750struct SimpleRng {
751    state: u64,
752}
753
754impl SimpleRng {
755    fn new(seed: u64) -> Self {
756        Self { state: seed.max(1) }
757    }
758
759    fn next_u64(&mut self) -> u64 {
760        self.state ^= self.state << 13;
761        self.state ^= self.state >> 7;
762        self.state ^= self.state << 17;
763        self.state
764    }
765
766    fn next_f64(&mut self) -> f64 {
767        (self.next_u64() as f64) / (u64::MAX as f64)
768    }
769
770    /// Generate standard normal using Box-Muller transform.
771    fn normal(&mut self) -> f64 {
772        let u1 = self.next_f64().max(1e-10);
773        let u2 = self.next_f64();
774        (-2.0 * u1.ln()).sqrt() * (2.0 * std::f64::consts::PI * u2).cos()
775    }
776}
777
778#[cfg(test)]
779mod tests {
780    use super::*;
781
782    fn create_simple_portfolio() -> Portfolio {
783        // Two-asset portfolio with known properties
784        Portfolio::new(
785            vec![1, 2],
786            vec![100_000.0, 100_000.0],
787            vec![0.08, 0.10],         // 8%, 10% expected returns
788            vec![0.15, 0.20],         // 15%, 20% volatility
789            vec![1.0, 0.5, 0.5, 1.0], // Correlation matrix
790        )
791    }
792
793    fn create_diversified_portfolio() -> Portfolio {
794        // Four-asset portfolio with low correlations
795        Portfolio::new(
796            vec![1, 2, 3, 4],
797            vec![50_000.0, 30_000.0, 15_000.0, 5_000.0],
798            vec![0.06, 0.08, 0.10, 0.12],
799            vec![0.10, 0.15, 0.20, 0.25],
800            vec![
801                1.0, 0.2, 0.1, 0.0, 0.2, 1.0, 0.3, 0.1, 0.1, 0.3, 1.0, 0.2, 0.0, 0.1, 0.2, 1.0,
802            ],
803        )
804    }
805
806    #[test]
807    fn test_monte_carlo_var_metadata() {
808        let kernel = MonteCarloVaR::new();
809        assert_eq!(kernel.metadata().id, "risk/monte-carlo-var");
810        assert_eq!(kernel.metadata().domain, Domain::RiskAnalytics);
811    }
812
813    #[test]
814    fn test_monte_carlo_var_calculation() {
815        let portfolio = create_simple_portfolio();
816        let params = VaRParams::new(0.95, 10, 10_000);
817
818        let result = MonteCarloVaR::compute(&portfolio, params);
819
820        assert!(result.var > 0.0, "VaR should be positive");
821        assert!(result.expected_shortfall >= result.var, "ES >= VaR");
822        assert_eq!(result.confidence_level, 0.95);
823        assert_eq!(result.holding_period, 10);
824
825        // VaR should be reasonable (not more than 50% of portfolio)
826        assert!(
827            result.var < 100_000.0,
828            "VaR seems too large: {}",
829            result.var
830        );
831    }
832
833    #[test]
834    fn test_var_increases_with_holding_period() {
835        let portfolio = create_simple_portfolio();
836
837        let var_1d = MonteCarloVaR::compute(&portfolio, VaRParams::new(0.95, 1, 10_000));
838        let var_10d = MonteCarloVaR::compute(&portfolio, VaRParams::new(0.95, 10, 10_000));
839
840        // 10-day VaR should be roughly sqrt(10) times 1-day VaR
841        let ratio = var_10d.var / var_1d.var;
842        assert!(
843            ratio > 2.5 && ratio < 4.0,
844            "VaR scaling ratio should be ~sqrt(10): {}",
845            ratio
846        );
847    }
848
849    #[test]
850    fn test_var_increases_with_confidence() {
851        let portfolio = create_simple_portfolio();
852
853        let var_95 = MonteCarloVaR::compute(&portfolio, VaRParams::new(0.95, 10, 10_000));
854        let var_99 = MonteCarloVaR::compute(&portfolio, VaRParams::new(0.99, 10, 10_000));
855
856        assert!(
857            var_99.var > var_95.var,
858            "99% VaR should exceed 95% VaR: {} vs {}",
859            var_99.var,
860            var_95.var
861        );
862    }
863
864    #[test]
865    fn test_component_var_sums_to_total() {
866        let portfolio = create_diversified_portfolio();
867        let params = VaRParams::new(0.95, 10, 10_000);
868
869        let result = MonteCarloVaR::compute(&portfolio, params);
870
871        // Component VaRs should roughly sum to total VaR
872        let component_sum: f64 = result.component_var.iter().sum();
873        let diff = (component_sum - result.var).abs() / result.var;
874
875        assert!(
876            diff < 0.20, // Allow 20% tolerance due to Monte Carlo noise
877            "Component VaR sum should be close to total: {} vs {}",
878            component_sum,
879            result.var
880        );
881    }
882
883    #[test]
884    fn test_portfolio_aggregation_metadata() {
885        let kernel = PortfolioRiskAggregation::new();
886        assert_eq!(kernel.metadata().id, "risk/portfolio-aggregation");
887    }
888
889    #[test]
890    fn test_portfolio_aggregation() {
891        let portfolio = create_diversified_portfolio();
892        let result = PortfolioRiskAggregation::compute(&portfolio, 0.95, 10);
893
894        assert!(result.portfolio_var > 0.0);
895        assert!(result.undiversified_var > result.portfolio_var);
896        assert!(result.diversification_benefit > 0.0);
897    }
898
899    #[test]
900    fn test_diversification_benefit() {
901        let portfolio = create_diversified_portfolio();
902        let result = PortfolioRiskAggregation::compute(&portfolio, 0.95, 10);
903
904        // Diversification benefit = undiversified - portfolio
905        assert!(
906            (result.diversification_benefit - (result.undiversified_var - result.portfolio_var))
907                .abs()
908                < 1.0
909        );
910
911        // With low correlations, should have significant benefit
912        let benefit_pct = result.diversification_benefit / result.undiversified_var;
913        assert!(
914            benefit_pct > 0.10,
915            "Should have >10% diversification benefit: {}%",
916            benefit_pct * 100.0
917        );
918    }
919
920    #[test]
921    fn test_risk_contributions_sum() {
922        let portfolio = create_diversified_portfolio();
923        let result = PortfolioRiskAggregation::compute(&portfolio, 0.95, 10);
924
925        let contrib_sum: f64 = result.risk_contributions.iter().sum();
926
927        // Risk contributions should sum to portfolio VaR
928        let diff = (contrib_sum - result.portfolio_var).abs() / result.portfolio_var;
929        assert!(
930            diff < 0.01,
931            "Risk contributions should sum to portfolio VaR: {} vs {}",
932            contrib_sum,
933            result.portfolio_var
934        );
935    }
936
937    #[test]
938    fn test_empty_portfolio() {
939        let empty = Portfolio::new(Vec::new(), Vec::new(), Vec::new(), Vec::new(), Vec::new());
940
941        let var_result = MonteCarloVaR::compute(&empty, VaRParams::default());
942        assert_eq!(var_result.var, 0.0);
943
944        let agg_result = PortfolioRiskAggregation::compute(&empty, 0.95, 10);
945        assert_eq!(agg_result.portfolio_var, 0.0);
946    }
947
948    #[test]
949    fn test_covariance_matrix() {
950        let portfolio = create_simple_portfolio();
951        let result = PortfolioRiskAggregation::compute(&portfolio, 0.95, 10);
952
953        // Covariance matrix should be symmetric
954        assert_eq!(result.covariance_matrix.len(), 4);
955        assert!((result.covariance_matrix[1] - result.covariance_matrix[2]).abs() < 1e-10);
956
957        // Diagonal should be variance
958        let var1 = 0.15 * 0.15;
959        let var2 = 0.20 * 0.20;
960        assert!((result.covariance_matrix[0] - var1).abs() < 1e-10);
961        assert!((result.covariance_matrix[3] - var2).abs() < 1e-10);
962    }
963}