Skip to content

Commit

Permalink
feat: update heston calibrator
Browse files Browse the repository at this point in the history
  • Loading branch information
@dancixx
dancixx committed Oct 9, 2024
1 parent d48d52e commit 03dc80e
Show file tree
Hide file tree
Showing 2 changed files with 106 additions and 89 deletions.
183 changes: 96 additions & 87 deletions src/quant/heston/calibrator.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@ use std::cell::RefCell;
use impl_new_derive::ImplNew;
use levenberg_marquardt::{LeastSquaresProblem, LevenbergMarquardt};
use nalgebra::{DMatrix, DVector, Dyn, Owned};
use ndarray::{Array, Array1};
use ndarray::Array1;

use crate::{
quant::{r#trait::Pricer, OptionType},
Expand All @@ -12,6 +12,40 @@ use crate::{

use super::pricer::HestonPricer;

/// Heston model parameters
#[derive(Clone, Debug)]
pub struct HestonParams {
pub v0: f64,
pub theta: f64,
pub rho: f64,
pub kappa: f64,
pub sigma: f64,
}

impl From<HestonParams> for DVector<f64> {
fn from(params: HestonParams) -> Self {
DVector::from_vec(vec![
params.v0,
params.theta,
params.rho,
params.kappa,
params.sigma,
])
}
}

impl From<DVector<f64>> for HestonParams {
fn from(params: DVector<f64>) -> Self {
HestonParams {
v0: params[0],
theta: params[1],
rho: params[2],
kappa: params[3],
sigma: params[4],
}
}
}

/// Heston calibrator
#[derive(ImplNew)]
pub struct HestonCalibrator {
Expand All @@ -30,15 +64,11 @@ pub struct HestonCalibrator {
/// Option type
pub option_type: OptionType,
/// Initial guess for the calibration from the NMLE method
pub initial_params: Option<DVector<f64>>,
pub initial_params: Option<HestonParams>,
}

impl HestonCalibrator {
pub fn calibrate(&mut self) {
if self.initial_params.is_none() {
panic!("Initial guess for the calibration is required. \n Use the initial_params method to set the initial guess \n or use the initial_params argument in the constructor.");
}

println!("Initial guess: {:?}", self.initial_params.as_ref().unwrap());

let (result, report) = LevenbergMarquardt::new().minimize(HestonCalibrationProblem::new(
Expand All @@ -47,36 +77,40 @@ impl HestonCalibrator {
self.s.clone().into(),
self.k.clone().into(),
self.tau,
self.r,
self.q,
&self.option_type,
));

// Print the c_market
println!("Market prices: {:?}", self.c_market);

let residuals = result.residuals().unwrap();

// Print the c_model
println!("Model prices: {:?}", result.residuals().unwrap());
println!(
"Model prices: {:?}",
DVector::from_vec(self.c_market.clone()) + residuals
);

// Print the result of the calibration
println!("Calibration report: {:?}", result.params);

// Print the result of the calibration
println!("Calibration report: {:?}", report);
}

/// Initial guess for the calibration
/// http://scis.scichina.com/en/2018/042202.pdf
///
/// Using NMLE (Normal Maximum Likelihood Estimation) method
pub fn initial_params(&mut self, s: Array1<f64>, v: Array1<f64>, r: f64) {
self.initial_params = Some(DVector::from_vec(nmle_heston(s, v, r)));
self.initial_params = Some(nmle_heston(s, v, r));
}
}

/// A calibrator.
#[derive(ImplNew)]
pub(crate) struct HestonCalibrationProblem<'a> {
/// Params to calibrate.
pub params: DVector<f64>,
pub params: HestonParams,
/// Option prices from the market.
pub c_market: DVector<f64>,
/// Asset price vector.
Expand All @@ -85,6 +119,10 @@ pub(crate) struct HestonCalibrationProblem<'a> {
pub k: DVector<f64>,
/// Time to maturity.
pub tau: f64,
/// Risk-free rate.
pub r: f64,
/// Dividend yield.
pub q: Option<f64>,
/// Option type
pub option_type: &'a OptionType,
/// Derivate matrix.
Expand All @@ -97,11 +135,11 @@ impl<'a> LeastSquaresProblem<f64, Dyn, Dyn> for HestonCalibrationProblem<'a> {
type ResidualStorage = Owned<f64, Dyn>;

fn set_params(&mut self, params: &DVector<f64>) {
self.params.copy_from(params);
self.params = HestonParams::from(params.clone());
}

fn params(&self) -> DVector<f64> {
self.params.clone()
self.params.clone().into()
}

fn residuals(&self) -> Option<DVector<f64>> {
Expand All @@ -111,14 +149,14 @@ impl<'a> LeastSquaresProblem<f64, Dyn, Dyn> for HestonCalibrationProblem<'a> {
for (idx, _) in self.c_market.iter().enumerate() {
let pricer = HestonPricer::new(
self.s[idx],
self.params[0],
self.params.v0,
self.k[idx],
0.5,
None,
self.params[2],
self.params[3],
self.params[1],
self.params[4],
self.r,
self.q,
self.params.rho,
self.params.kappa,
self.params.theta,
self.params.sigma,
None,
self.tau,
None,
Expand All @@ -135,8 +173,6 @@ impl<'a> LeastSquaresProblem<f64, Dyn, Dyn> for HestonCalibrationProblem<'a> {
}

let _ = std::mem::replace(&mut *self.derivates.borrow_mut(), derivates);
// println!("c_model: {:?}", c_model);
// println!("c_market: {:?}", self.c_market.clone());
Some(c_model - self.c_market.clone())
}

Expand All @@ -146,81 +182,54 @@ impl<'a> LeastSquaresProblem<f64, Dyn, Dyn> for HestonCalibrationProblem<'a> {

// The Jacobian matrix is a matrix of partial derivatives
// of the residuals with respect to the parameters.
let jacobian = DMatrix::from_vec(
derivates.len() / self.params.len(),
self.params.len(),
derivates,
);
let jacobian = DMatrix::from_vec(derivates.len() / 5, 5, derivates);

Some(jacobian)
}
}

#[cfg(test)]
mod tests {

use polars::series::ChunkCompare;

use crate::{quant::yahoo::Yahoo, stochastic::N};

use super::*;

#[test]
fn test_heston_calibrate() {
let mut yahoo = Yahoo::default();
yahoo.set_symbol("AAPL");
yahoo.get_options_chain(&OptionType::Call);
yahoo.get_price_history();
let options = yahoo.options.as_ref().unwrap();
// remove where last price less than 0.01
let mask = options.column("last_price").unwrap().gt(10.00).unwrap();
let options = options.filter(&mask).unwrap();

// Get Price history
let price_history = yahoo.price_history.as_ref().unwrap();
let s = price_history.select(["close"]).unwrap();
let s = s
.select_at_idx(0)
.unwrap()
.f64()
.unwrap()
.into_no_null_iter()
.collect::<Vec<f64>>();
// convert to years the epoch time
let tau = (yahoo.options_chain.as_ref().unwrap().option_chain.result[0].options[0]
.expiration_date as f64
- chrono::Local::now().timestamp() as f64)
/ 31536000.0;
let tau = 24.0 / 365.0;
println!("Time to maturity: {}", tau);
let c_market = options.select(["last_price"]).unwrap();
let c_market = c_market
.select_at_idx(0)
.unwrap()
.f64()
.unwrap()
.into_no_null_iter()
.collect::<Vec<f64>>();

let k = options.select(["strike"]).unwrap();
let k = k
.select_at_idx(0)
.unwrap()
.f64()
.unwrap()
.into_no_null_iter()
.collect::<Vec<f64>>();

let mut calibrator = HestonCalibrator::new(
s.clone(),
k.clone(),
0.5,
None,
c_market,
0.5,
OptionType::Call,
None,
);
calibrator.initial_params(Array1::from(s), Array1::default(N), 0.05);
calibrator.calibrate();

let s = vec![
425.73, 425.73, 425.73, 425.67, 425.68, 425.65, 425.65, 425.68, 425.65, 425.16, 424.78,
425.19,
];

let k = vec![
395.0, 400.0, 405.0, 410.0, 415.0, 420.0, 425.0, 430.0, 435.0, 440.0, 445.0, 450.0,
];

let c_market = vec![
30.75, 25.88, 21.00, 16.50, 11.88, 7.69, 4.44, 2.10, 0.78, 0.25, 0.10, 0.10,
];

let v0 = Array1::linspace(0.0, 0.01, 10);

for v in v0.iter() {
let mut calibrator = HestonCalibrator::new(
s.clone(),
k.clone(),
6.40e-4,
None,
c_market.clone(),
tau,
OptionType::Call,
Some(HestonParams {
v0: *v,
theta: 6.47e-5,
rho: -1.98e-3,
kappa: 6.57e-3,
sigma: 5.09e-4,
}),
);
calibrator.calibrate();
}
}
}
12 changes: 10 additions & 2 deletions src/stats/mle.rs
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,7 @@
use ndarray::Array1;

use crate::quant::heston::calibrator::HestonParams;

/// Maximum likelihood estimation for Heston model
/// http://scis.scichina.com/en/2018/042202.pdf
///
Expand All @@ -10,7 +12,7 @@ use ndarray::Array1;
///
/// # Returns
/// Vec<f64> - estimated parameters
pub fn nmle_heston(s: Array1<f64>, v: Array1<f64>, r: f64) -> Vec<f64> {
pub fn nmle_heston(s: Array1<f64>, v: Array1<f64>, r: f64) -> HestonParams {
let n = v.len();
let delta = 1.0 / n as f64;
let mut sum = [0.0; 4];
Expand Down Expand Up @@ -58,5 +60,11 @@ pub fn nmle_heston(s: Array1<f64>, v: Array1<f64>, r: f64) -> Vec<f64> {

let rho_hat = sum_dw1dw2 / (n as f64 * delta);

vec![v[0], theta_hat, rho_hat, kappa_hat, sigma_hat]
HestonParams {
v0: v[0],
theta: theta_hat,
rho: rho_hat,
kappa: kappa_hat,
sigma: sigma_hat,
}
}

0 comments on commit 03dc80e

Please sign in to comment.