Add mutable accessors and fix floating asian options (#51)

* Add mutable accessors and fix floating asian options

* Add simulation methods for arithmetic and geometric averages with mutable access

* Update tests

* Update tests

* Update tests

Author: carlobortolan

README.md

@@ -54,7 +54,7 @@ Quantrs supports options pricing with various models for both vanilla and exotic
 
 > ¹ _"Exotic" options with standard exercise style; only differ in their payoff value_\
 > ² _Non-vanilla path-dependent "exotic" options_\
-> ³ _MC simulates underlying price paths based on geometric Brownian motion for Black-Scholes models and geometric average price paths for Asian and Lookback options_\
+> ³ _MC simulates underlying price paths based on geometric Brownian motion for Black-Scholes models and both arithmetic or geometric average price paths for Asian and Lookback options_\
 > ✅ = Supported, ⏳ = Planned / In progress, ❌ = Not supported / Not applicable
 
 <!--Bachelier and Modified Bachelier-->

examples/options_pricing.rs

@@ -278,15 +278,15 @@ fn example_strategy() {
     println!(
         "[Covered Call: {:?}], given stock: {}, call: {}",
         model.covered_call(&instrument, &call)(50.0),
-        instrument.spot,
+        instrument.spot(),
         model.price(&call)
     );
 
     let put = EuropeanOption::new(instrument.clone(), 40.0, 1.0, Put);
     println!(
         "[Protective Put: {:?}], given stock: {}, put: {}",
         model.protective_put(&instrument, &put)(50.0),
-        instrument.spot,
+        instrument.spot(),
         model.price(&put)
     );
 
@@ -295,7 +295,7 @@ fn example_strategy() {
     println!(
         "[Collar: {:?}], given stock: {}, put: {}, call: {}",
         model.collar(&instrument, &otm_put, &otm_call)(50.0),
-        instrument.spot,
+        instrument.spot(),
         model.price(&otm_put),
         model.price(&otm_call)
     );
@@ -306,7 +306,7 @@ fn example_strategy() {
     println!(
         "[Fence: {:?}], given stock: {}, put: {}, put: {}, call: {}",
         model.fence(&instrument, &atm_put, &otm_put, &otm_call)(50.0),
-        instrument.spot,
+        instrument.spot(),
         model.price(&atm_put),
         model.price(&otm_call),
         model.price(&otm_put)

src/options/instrument.rs

@@ -29,7 +29,7 @@ use rand_distr::{Distribution, Normal};
 #[derive(Debug, Default, Clone)]
 pub struct Instrument {
     /// Current price of the underlying asset or future price at time 0.
-    pub spot: f64,
+    pub spot: Vec<f64>,
     /// Maximum spot price of the underlying asset.
     pub max_spot: f64,
     /// Minimum spot price of the underlying asset.
@@ -54,7 +54,7 @@ impl Instrument {
     /// A new `Instrument`.
     pub fn new() -> Self {
         Self {
-            spot: 0.0,
+            spot: vec![0.0],
             max_spot: 0.0,
             min_spot: 0.0,
             continuous_dividend_yield: 0.0,
@@ -75,6 +75,20 @@ impl Instrument {
     ///
     /// The instrument with the spot price set.
     pub fn with_spot(mut self, spot: f64) -> Self {
+        self.spot = vec![spot];
+        self
+    }
+
+    /// Set the spot price of the instrument.
+    ///
+    /// # Arguments
+    ///
+    /// * `spot` - The spot price of the instrument (i.e., the current price of the underlying asset or future price at time 0).
+    ///
+    /// # Returns
+    ///
+    /// The instrument with the spot price set.
+    pub fn with_spots(mut self, spot: Vec<f64>) -> Self {
         self.spot = spot;
         self
     }
@@ -173,7 +187,9 @@ impl Instrument {
 
         let weight = 1.0 / assets.len() as f64;
         self.assets = assets.iter().map(|asset| (asset.clone(), weight)).collect();
-        self.spot = self.assets.iter().map(|(a, w)| a.spot * w).sum::<f64>();
+        let new_spot = self.assets.iter().map(|(a, w)| a.spot() * w).sum::<f64>();
+        self.spot = vec![new_spot];
+
         self.sort_assets_by_performance();
         self
     }
@@ -201,7 +217,7 @@ impl Instrument {
     pub fn sort_assets_by_performance(&mut self) {
         self.assets
             .sort_by(|a, b| b.0.spot.partial_cmp(&a.0.spot).unwrap());
-        self.spot = self.assets.iter().map(|(a, w)| a.spot * w).sum::<f64>();
+        self.spot = vec![(self.assets.iter().map(|(a, w)| a.spot() * w).sum::<f64>())];
         self.sorted = true;
     }
 
@@ -236,6 +252,36 @@ impl Instrument {
         &self.assets.last().unwrap().0
     }
 
+    /// Get the best-performing asset (mutable).
+    ///
+    /// # Returns
+    ///
+    /// The best-performing asset.
+    pub fn best_performer_mut(&mut self) -> &mut Instrument {
+        if self.assets.is_empty() {
+            return self;
+        }
+        if !self.sorted {
+            panic!("Assets are not sorted");
+        }
+        &mut self.assets.first_mut().unwrap().0
+    }
+
+    /// Get the worst-performing asset (mutable).
+    ///
+    /// # Returns
+    ///
+    /// The worst-performing asset.
+    pub fn worst_performer_mut(&mut self) -> &mut Instrument {
+        if self.assets.is_empty() {
+            return self;
+        }
+        if !self.sorted {
+            panic!("Assets are not sorted");
+        }
+        &mut self.assets.last_mut().unwrap().0
+    }
+
     /// Calculate the adjusted spot price.
     ///
     /// # Arguments
@@ -248,7 +294,25 @@ impl Instrument {
     /// The adjusted spot price.
     pub fn calculate_adjusted_spot(&self, ttm: f64) -> f64 {
         let n_dividends = self.dividend_times.iter().filter(|&&t| t <= ttm).count() as f64;
-        self.spot * (1.0 - self.discrete_dividend_yield).powf(n_dividends)
+        self.spot() * (1.0 - self.discrete_dividend_yield).powf(n_dividends)
+    }
+
+    /// Return current spot price.
+    ///
+    /// # Returns
+    ///
+    /// The current spot price.
+    pub fn spot(&self) -> f64 {
+        *self.spot.first().unwrap()
+    }
+
+    /// Return terminal spot price.
+    ///
+    /// # Returns
+    ///
+    /// The terminal spot price.
+    pub fn terminal_spot(&self) -> f64 {
+        *self.spot.last().unwrap()
     }
 
     /// Simulate random asset prices (Euler method)
@@ -267,11 +331,12 @@ impl Instrument {
         rng: &mut ThreadRng,
         risk_free_rate: f64,
         volatility: f64,
+        steps: usize,
     ) -> Vec<f64> {
         let normal = Normal::new(0.0, 1.0).unwrap();
-        let dt: f64 = 1.0 / 252.0; // Daily time step
-        let mut prices = vec![self.spot; 252];
-        for i in 1..252 {
+        let dt: f64 = 1.0 / steps as f64; // Daily time step
+        let mut prices = vec![self.spot(); steps];
+        for i in 1..steps {
             let z = normal.sample(rng);
             prices[i] = prices[i - 1]
                 * (1.0
@@ -304,14 +369,14 @@ impl Instrument {
     ) -> Vec<f64> {
         let dt = time_to_maturity / steps as f64; // Time step
         let normal: Normal<f64> = Normal::new(0.0, dt.sqrt()).unwrap(); // Adjusted standard deviation
-        let mut logs = vec![self.spot.ln(); steps];
+        let mut logs = vec![self.spot().ln(); steps];
         for i in 1..steps {
             let z = normal.sample(rng);
             logs[i] = logs[i - 1]
                 + (risk_free_rate - self.continuous_dividend_yield - 0.5 * volatility.powi(2)) * dt
                 + volatility * z;
         }
-        logs
+        logs.iter().map(|log| log.exp()).collect()
     }
 
     /// Average asset prices
@@ -337,19 +402,15 @@ impl Instrument {
         risk_free_rate: f64,
         steps: usize,
     ) -> f64 {
-        let prices: Vec<f64> = match method {
+        let prices = match method {
             SimMethod::Milstein => unimplemented!("Milstein method not implemented"),
-            SimMethod::Euler => self.euler_simulation(rng, risk_free_rate, volatility),
+            SimMethod::Euler => self.euler_simulation(rng, risk_free_rate, volatility, steps),
             SimMethod::Log => {
                 self.log_simulation(rng, volatility, time_to_maturity, risk_free_rate, steps)
             }
         };
 
-        let res = prices.iter().sum::<f64>() / (prices.len()) as f64;
-        match method {
-            SimMethod::Log => res.exp(),
-            _ => res,
-        }
+        prices.iter().sum::<f64>() / (prices.len()) as f64
     }
 
     /// Geometric average asset prices
@@ -375,18 +436,83 @@ impl Instrument {
         risk_free_rate: f64,
         steps: usize,
     ) -> f64 {
-        let prices: Vec<f64> = match method {
+        let prices = match method {
             SimMethod::Milstein => unimplemented!("Milstein method not implemented"),
-            SimMethod::Euler => self.euler_simulation(rng, risk_free_rate, volatility),
+            SimMethod::Euler => self.euler_simulation(rng, risk_free_rate, volatility, steps),
             SimMethod::Log => {
                 self.log_simulation(rng, volatility, time_to_maturity, risk_free_rate, steps)
             }
         };
 
-        match method {
-            SimMethod::Log => (prices.iter().sum::<f64>() / prices.len() as f64).exp(),
-            _ => (prices.iter().map(|price| price.ln()).sum::<f64>() / prices.len() as f64).exp(),
-        }
+        (self.spot.iter().map(|price| price.ln()).sum::<f64>() / self.spot.len() as f64).exp()
+    }
+
+    /// Average asset prices (mutates the underlying instrument)
+    ///
+    /// # Arguments
+    ///
+    /// * `rng` - Random number generator.
+    /// * `method` - Simulation method.
+    /// * `volatility` - Volatility.
+    /// * `time_to_maturity` - Time to maturity.
+    /// * `risk_free_rate` - Risk-free rate.
+    /// * `steps` - Number of steps.
+    ///
+    /// # Returns
+    ///
+    /// The average asset price.
+    pub fn simulate_arithmetic_average_mut(
+        &mut self,
+        rng: &mut ThreadRng,
+        method: SimMethod,
+        volatility: f64,
+        time_to_maturity: f64,
+        risk_free_rate: f64,
+        steps: usize,
+    ) -> f64 {
+        self.spot = match method {
+            SimMethod::Milstein => unimplemented!("Milstein method not implemented"),
+            SimMethod::Euler => self.euler_simulation(rng, risk_free_rate, volatility, steps),
+            SimMethod::Log => {
+                self.log_simulation(rng, volatility, time_to_maturity, risk_free_rate, steps)
+            }
+        };
+
+        self.spot.iter().sum::<f64>() / (self.spot.len()) as f64
+    }
+
+    /// Geometric asset prices (mutates the underlying instrument)
+    ///
+    /// # Arguments
+    ///
+    /// * `rng` - Random number generator.
+    /// * `method` - Simulation method.
+    /// * `volatility` - Volatility.
+    /// * `time_to_maturity` - Time to maturity.
+    /// * `risk_free_rate` - Risk-free rate.
+    /// * `steps` - Number of steps.
+    ///
+    /// # Returns
+    ///
+    /// The geometric average asset price.
+    pub fn simulate_geometric_average_mut(
+        &mut self,
+        rng: &mut ThreadRng,
+        method: SimMethod,
+        volatility: f64,
+        time_to_maturity: f64,
+        risk_free_rate: f64,
+        steps: usize,
+    ) -> f64 {
+        self.spot = match method {
+            SimMethod::Milstein => unimplemented!("Milstein method not implemented"),
+            SimMethod::Euler => self.euler_simulation(rng, risk_free_rate, volatility, steps),
+            SimMethod::Log => {
+                self.log_simulation(rng, volatility, time_to_maturity, risk_free_rate, steps)
+            }
+        };
+
+        (self.spot.iter().map(|price| price.ln()).sum::<f64>() / self.spot.len() as f64).exp()
     }
 
     /// Directly simulate the asset price using the geometric Brownian motion formula
@@ -412,7 +538,7 @@ impl Instrument {
     ) -> f64 {
         let normal = Normal::new(0.0, 1.0).unwrap();
         let dt = time_to_maturity / steps as f64;
-        let mut price = self.spot;
+        let mut price = self.spot();
         for _ in 0..steps {
             let z = normal.sample(rng);
             price *= ((risk_free_rate - self.continuous_dividend_yield - 0.5 * volatility.powi(2))

src/options/models/binomial_tree.rs

@@ -127,7 +127,7 @@ impl OptionPricing for BinomialTreeModel {
         let mut option_values: Vec<f64> = (0..=self.steps)
             .map(|i| {
                 option.payoff(Some(
-                    option.instrument().spot * u.powi(i as i32) * d.powi((self.steps - i) as i32),
+                    option.instrument().spot() * u.powi(i as i32) * d.powi((self.steps - i) as i32),
                 ))
             })
             .collect();
@@ -140,7 +140,7 @@ impl OptionPricing for BinomialTreeModel {
 
                 if matches!(option.style(), OptionStyle::American) {
                     let early_exercise = option.payoff(Some(
-                        option.instrument().spot * u.powi(i as i32) * d.powi((step - i) as i32),
+                        option.instrument().spot() * u.powi(i as i32) * d.powi((step - i) as i32),
                     ));
                     option_values[i] = expected_value.max(early_exercise);
                 } else {
@@ -150,7 +150,7 @@ impl OptionPricing for BinomialTreeModel {
         }
 
         if matches!(option.style(), OptionStyle::American) {
-            option_values[0].max(option.strike() - option.instrument().spot) // TODO: Change to max(0.0, self.payoff(Some(self.spot)))
+            option_values[0].max(option.strike() - option.instrument().spot()) // TODO: Change to max(0.0, self.payoff(Some(self.spot)))
         } else {
             option_values[0] // Return the root node value
         }

src/options/models/black_76.rs

@@ -80,7 +80,7 @@ impl Black76Model {
         let (d1, d2) = self.calculate_d1_d2(instrument, strike, ttm);
 
         (-self.risk_free_rate * ttm).exp()
-            * (instrument.spot * normal.cdf(d1) - strike * normal.cdf(d2))
+            * (instrument.spot() * normal.cdf(d1) - strike * normal.cdf(d2))
     }
 
     /// Calculate the price of a European put option using the Black-76 formula.
@@ -105,7 +105,7 @@ impl Black76Model {
         let (d1, d2) = self.calculate_d1_d2(instrument, strike, ttm);
 
         (-self.risk_free_rate * ttm).exp()
-            * (strike * normal.cdf(-d2) - instrument.spot * normal.cdf(-d1))
+            * (strike * normal.cdf(-d2) - instrument.spot() * normal.cdf(-d1))
     }
 }