1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2011 Klaus Spanderen
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/*! \file gemanroncoroniprocess.cpp
21 \brief Geman-Roncoroni process
22*/
23
24#include <ql/math/functional.hpp>
25#include <ql/processes/eulerdiscretization.hpp>
26#include <ql/experimental/processes/gemanroncoroniprocess.hpp>
27
28
29namespace QuantLib {
30
31 GemanRoncoroniProcess::GemanRoncoroniProcess(
32 Real x0,
33 Real alpha, Real beta,
34 Real gamma, Real delta,
35 Real eps, Real zeta, Real d,
36 Real k, Real tau,
37 Real sig2, Real a, Real b,
38 Real theta1, Real theta2, Real theta3,
39 Real psi)
40 : StochasticProcess1D(ext::shared_ptr<discretization>(
41 new EulerDiscretization)),
42 x0_(x0),
43 alpha_(alpha), beta_(beta),
44 gamma_(gamma), delta_(delta),
45 eps_(eps), zeta_(zeta), d_(d),
46 k_(k), tau_(tau),
47 sig2_(sig2), a_(a), b_(b),
48 theta1_(theta1), theta2_(theta2), theta3_(theta3),
49 psi_(psi) {
50 }
51
52 Real GemanRoncoroniProcess::x0() const {
53 return x0_;
54 }
55
56 Real GemanRoncoroniProcess::drift(Time t, Real x) const {
57 const Real mu = alpha_ + beta_*t + gamma_*std::cos(x: eps_+2*M_PI*t)
58 + delta_*std::cos(x: zeta_+4*M_PI*t);
59 const Real muPrime = beta_ - gamma_*2*M_PI*std::sin(x: eps_+2*M_PI*t)
60 - delta_*4*M_PI*std::sin(x: zeta_+4*M_PI*t);
61
62 return muPrime + theta1_*(mu-x);
63 }
64
65 Real GemanRoncoroniProcess::diffusion(Time t, Real /*x*/) const {
66 return std::sqrt(x: sig2_ + a_*squared(x: std::cos(M_PI*t+b_)));
67 }
68
69 Real GemanRoncoroniProcess::stdDeviation(Time t0, Real /*x0*/, Time dt) const {
70 const Volatility sig2t = sig2_+a_*squared(x: std::cos(M_PI*t0+b_));
71
72 return std::sqrt(x: sig2t/(2*theta1_)*(1.0-std::exp(x: -2*theta1_*dt)));
73 }
74
75 Real GemanRoncoroniProcess::evolve(Time t0, Real x0,
76 Time dt, Real dw) const {
77 // random number generator for the jump part
78 if (!urng_) {
79 typedef PseudoRandom::urng_type urng_type;
80 urng_ = ext::make_shared<urng_type>(args: (unsigned long)(1234UL * dw + 56789UL));
81 }
82 Array du(3);
83 du[0] = urng_->next().value;
84 du[1] = urng_->next().value;
85
86 return evolve(t0, x0, dt, dw, du);
87 }
88
89 Real GemanRoncoroniProcess::evolve(Time t0, Real x0, Time dt,
90 Real dw, const Array& du) const {
91 Real retVal;
92 const Time t = t0 + 0.5*dt;
93 const Real mu = alpha_ + beta_*t + gamma_*std::cos(x: eps_ +2*M_PI*t)
94 + delta_*std::cos(x: zeta_+4*M_PI*t);
95
96 const Real j = -1.0/theta3_
97 *std::log(x: 1.0+du[1]*(std::exp(x: -theta3_*psi_)-1.0));
98
99 if (x0 <= mu+d_) {
100 retVal = StochasticProcess1D::evolve(t0: t, x0, dt, dw);
101
102 const Real jumpIntensity
103 = theta2_*(2.0/(1+std::fabs(x: std::sin(M_PI*(t-tau_)/k_)))-1);
104 const Time interarrival = -1.0/jumpIntensity*std::log(x: du[0]);
105
106 if (interarrival < dt) {
107 retVal += j;
108 }
109 }
110 else {
111 retVal = x0-j;
112 }
113
114 return retVal;
115 }
116}
117

source code of quantlib/ql/experimental/processes/gemanroncoroniprocess.cpp