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