1 /* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ 2 3 /* 4 Copyright (C) 2001, 2002, 2003 Sadruddin Rejeb 5 6 This file is part of QuantLib, a free-software/open-source library 7 for financial quantitative analysts and developers - http://quantlib.org/ 8 9 QuantLib is free software: you can redistribute it and/or modify it 10 under the terms of the QuantLib license. You should have received a 11 copy of the license along with this program; if not, please email 12 <quantlib-dev@lists.sf.net>. The license is also available online at 13 <http://quantlib.org/license.shtml>. 14 15 This program is distributed in the hope that it will be useful, but WITHOUT 16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 17 FOR A PARTICULAR PURPOSE. See the license for more details. 18 */ 19 20 #include <ql/models/shortrate/onefactormodels/coxingersollross.hpp> 21 #include <ql/methods/lattices/trinomialtree.hpp> 22 #include <ql/math/distributions/chisquaredistribution.hpp> 23 24 namespace QuantLib { 25 26 class CoxIngersollRoss::VolatilityConstraint : public Constraint { 27 private: 28 class Impl : public Constraint::Impl { 29 Real k_, theta_; 30 public: Impl(Real k,Real theta)31 Impl(Real k, Real theta) : k_(k), theta_(theta) {} test(const Array & params) const32 bool test(const Array& params) const { 33 Real sigma = params[0]; 34 return (sigma > 0.0) && (sigma*sigma < 2.0*k_*theta_); 35 } 36 }; 37 public: VolatilityConstraint(Real k,Real theta)38 VolatilityConstraint(Real k, Real theta) 39 : Constraint(ext::shared_ptr<Constraint::Impl>( 40 new VolatilityConstraint::Impl(k, theta))) {} 41 }; 42 CoxIngersollRoss(Rate r0,Real theta,Real k,Real sigma,bool withFellerConstraint)43 CoxIngersollRoss::CoxIngersollRoss(Rate r0, Real theta, 44 Real k, Real sigma, 45 bool withFellerConstraint) 46 : OneFactorAffineModel(4), 47 theta_(arguments_[0]), k_(arguments_[1]), 48 sigma_(arguments_[2]), r0_(arguments_[3]) { 49 theta_ = ConstantParameter(theta, PositiveConstraint()); 50 k_ = ConstantParameter(k, PositiveConstraint()); 51 if (withFellerConstraint) 52 sigma_ = ConstantParameter(sigma, VolatilityConstraint(k,theta)); 53 else 54 sigma_ = ConstantParameter(sigma, PositiveConstraint()); 55 r0_ = ConstantParameter(r0, PositiveConstraint()); 56 } 57 58 ext::shared_ptr<OneFactorModel::ShortRateDynamics> dynamics() const59 CoxIngersollRoss::dynamics() const { 60 return ext::shared_ptr<ShortRateDynamics>( 61 new Dynamics(theta(), k() , sigma(), x0())); 62 } 63 A(Time t,Time T) const64 Real CoxIngersollRoss::A(Time t, Time T) const { 65 Real sigma2 = sigma()*sigma(); 66 Real h = std::sqrt(k()*k() + 2.0*sigma2); 67 Real numerator = 2.0*h*std::exp(0.5*(k()+h)*(T-t)); 68 Real denominator = 2.0*h + (k()+h)*(std::exp((T-t)*h) - 1.0); 69 Real value = std::log(numerator/denominator)* 70 2.0*k()*theta()/sigma2; 71 return std::exp(value); 72 } 73 B(Time t,Time T) const74 Real CoxIngersollRoss::B(Time t, Time T) const { 75 Real h = std::sqrt(k()*k() + 2.0*sigma()*sigma()); 76 Real temp = std::exp((T-t)*h) - 1.0; 77 Real numerator = 2.0*temp; 78 Real denominator = 2.0*h + (k()+h)*temp; 79 Real value = numerator/denominator; 80 return value; 81 } 82 discountBondOption(Option::Type type,Real strike,Time t,Time s) const83 Real CoxIngersollRoss::discountBondOption(Option::Type type, 84 Real strike, 85 Time t, Time s) const { 86 87 QL_REQUIRE(strike>0.0, "strike must be positive"); 88 DiscountFactor discountT = discountBond(0.0, t, x0()); 89 DiscountFactor discountS = discountBond(0.0, s, x0()); 90 91 if (t < QL_EPSILON) { 92 switch(type) { 93 case Option::Call: 94 return std::max<Real>(discountS - strike, 0.0); 95 case Option::Put: 96 return std::max<Real>(strike - discountS, 0.0); 97 default: QL_FAIL("unsupported option type"); 98 } 99 } 100 101 Real sigma2 = sigma()*sigma(); 102 Real h = std::sqrt(k()*k() + 2.0*sigma2); 103 Real b = B(t,s); 104 105 Real rho = 2.0*h/(sigma2*(std::exp(h*t) - 1.0)); 106 Real psi = (k() + h)/sigma2; 107 108 Real df = 4.0*k()*theta()/sigma2; 109 Real ncps = 2.0*rho*rho*x0()*std::exp(h*t)/(rho+psi+b); 110 Real ncpt = 2.0*rho*rho*x0()*std::exp(h*t)/(rho+psi); 111 112 NonCentralCumulativeChiSquareDistribution chis(df, ncps); 113 NonCentralCumulativeChiSquareDistribution chit(df, ncpt); 114 115 Real z = std::log(A(t,s)/strike)/b; 116 Real call = discountS*chis(2.0*z*(rho+psi+b)) - 117 strike*discountT*chit(2.0*z*(rho+psi)); 118 119 if (type == Option::Call) 120 return call; 121 else 122 return call - discountS + strike*discountT; 123 } 124 125 ext::shared_ptr<Lattice> tree(const TimeGrid & grid) const126 CoxIngersollRoss::tree(const TimeGrid& grid) const { 127 ext::shared_ptr<TrinomialTree> trinomial( 128 new TrinomialTree(dynamics()->process(), grid, true)); 129 return ext::shared_ptr<Lattice>( 130 new ShortRateTree(trinomial, dynamics(), grid)); 131 } 132 133 } 134 135