QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
fdmndimsolver.hpp
Go to the documentation of this file.
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 fdmndimsolver.hpp
21*/
22
23#ifndef quantlib_fdm_n_dim_solver_hpp
24#define quantlib_fdm_n_dim_solver_hpp
25
36
37#include <numeric>
38
39namespace QuantLib {
40
41 template <Size N>
42 class FdmNdimSolver : public LazyObject {
43 public:
44 FdmNdimSolver(const FdmSolverDesc& solverDesc,
45 const FdmSchemeDesc& schemeDesc,
46 ext::shared_ptr<FdmLinearOpComposite> op);
47
48 void performCalculations() const override;
49
50 Real interpolateAt(const std::vector<Real>& x) const;
51 Real thetaAt(const std::vector<Real>& x) const;
52
53 // template meta programming
55 void static setValue(data_table& f,
56 const std::vector<Size>& x, Real value);
57
58 private:
61 const ext::shared_ptr<FdmLinearOpComposite> op_;
62
63 const ext::shared_ptr<FdmSnapshotCondition> thetaCondition_;
64 const ext::shared_ptr<FdmStepConditionComposite> conditions_;
65
66 std::vector<std::vector<Real> > x_;
67 std::vector<Real> initialValues_;
68 const std::vector<bool> extrapolation_;
69
70 mutable ext::shared_ptr<data_table> f_;
71 mutable ext::shared_ptr<MultiCubicSpline<N> > interp_;
72 };
73
74
75 template <Size N>
77 const FdmSchemeDesc& schemeDesc,
78 ext::shared_ptr<FdmLinearOpComposite> op)
79 : solverDesc_(solverDesc), schemeDesc_(schemeDesc), op_(std::move(op)),
80 thetaCondition_(new FdmSnapshotCondition(
81 0.99 * std::min(1.0 / 365.0,
82 solverDesc.condition->stoppingTimes().empty() ?
83 solverDesc.maturity :
84 solverDesc.condition->stoppingTimes().front()))),
85 conditions_(FdmStepConditionComposite::joinConditions(thetaCondition_, solverDesc.condition)),
86 x_(solverDesc.mesher->layout()->dim().size()),
87 initialValues_(solverDesc.mesher->layout()->size()),
88 extrapolation_(std::vector<bool>(N, false)) {
89
90 QL_REQUIRE(solverDesc.mesher->layout()->dim().size() == N, "solver dim " << N
91 << "does not fit to layout dim " << solverDesc.mesher->layout()->size());
92
93 for (Size i=0; i < N; ++i) {
94 x_[i].reserve(solverDesc.mesher->layout()->dim()[i]);
95 }
96
97 for (const auto& iter : *solverDesc.mesher->layout()) {
99 ->avgInnerValue(iter, solverDesc.maturity);
100
101 const std::vector<Size>& c = iter.coordinates();
102 for (Size i=0; i < N; ++i) {
103 if ((std::accumulate(c.begin(), c.end(), 0UL) - c[i]) == 0U) {
104 x_[i].push_back(solverDesc.mesher->location(iter, i));
105 }
106 }
107 }
108
109 f_ = ext::shared_ptr<data_table>(new data_table(x_));
110 }
111
112
113 template <Size N> inline
115 Array rhs(initialValues_.size());
116 std::copy(initialValues_.begin(), initialValues_.end(), rhs.begin());
117
118 FdmBackwardSolver(op_, solverDesc_.bcSet, conditions_, schemeDesc_)
119 .rollback(rhs, solverDesc_.maturity, 0.0,
120 solverDesc_.timeSteps, solverDesc_.dampingSteps);
121
122 for (const auto& iter : *solverDesc_.mesher->layout()) {
123 setValue(*f_, iter.coordinates(), rhs[iter.index()]);
124 }
125
126 interp_ = ext::shared_ptr<MultiCubicSpline<N> >(
127 new MultiCubicSpline<N>(x_, *f_, extrapolation_));
128 }
129
130
131 template <Size N> inline
132 Real FdmNdimSolver<N>::thetaAt(const std::vector<Real>& x) const {
133 if (conditions_->stoppingTimes().front() == 0.0)
134 return Null<Real>();
135
136 calculate();
137 const Array& rhs = thetaCondition_->getValues();
138
139 data_table f(x_);
140
141 for (const auto& iter : *solverDesc_.mesher->layout()) {
142 setValue(f, iter.coordinates(), rhs[iter.index()]);
143 }
144
145 return (MultiCubicSpline<N>(x_, f)(x)
146 - interpolateAt(x)) / thetaCondition_->getTime();
147 }
148
149 template <Size N> inline
150 Real FdmNdimSolver<N>::interpolateAt(const std::vector<Real>& x) const {
151 calculate();
152
153 return (*interp_)(x);
154 }
155
156 template <Size N> inline
158 const std::vector<Size>& x, Real value) {
159 FdmNdimSolver<N-1>::setValue(f[x[x.size()-N]], x, value);
160 }
161
162 template <> inline
164 const std::vector<Size>& x, Real value) {
165 f[x.back()] = value;
166 }
167}
168
169#endif
1-D array used in linear algebra.
Definition: array.hpp:52
const_iterator begin() const
Definition: array.hpp:503
void rollback(array_type &a, Time from, Time to, Size steps, Size dampingSteps)
void performCalculations() const override
std::vector< Real > initialValues_
static void setValue(data_table &f, const std::vector< Size > &x, Real value)
const ext::shared_ptr< FdmStepConditionComposite > conditions_
MultiCubicSpline< N >::data_table data_table
const ext::shared_ptr< FdmSnapshotCondition > thetaCondition_
ext::shared_ptr< data_table > f_
const FdmSolverDesc solverDesc_
Real interpolateAt(const std::vector< Real > &x) const
FdmNdimSolver(const FdmSolverDesc &solverDesc, const FdmSchemeDesc &schemeDesc, ext::shared_ptr< FdmLinearOpComposite > op)
Real thetaAt(const std::vector< Real > &x) const
std::vector< std::vector< Real > > x_
const std::vector< bool > extrapolation_
ext::shared_ptr< MultiCubicSpline< N > > interp_
const ext::shared_ptr< FdmLinearOpComposite > op_
const FdmSchemeDesc schemeDesc_
Framework for calculation on demand and result caching.
Definition: lazyobject.hpp:35
N-dimensional cubic spline interpolation between discrete points.
template class providing a null value for a given type.
Definition: null.hpp:76
ext::function< Real(Real)> f_
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
layer of abstraction to calculate the inner value
memory layout of a fdm linear operator
mesher for a fdm grid
step condition for value inspection
composite of fdm step conditions
generic finite difference model
QL_REAL Real
real number
Definition: types.hpp:50
std::size_t Size
size of a container
Definition: types.hpp:58
framework for calculation on demand and result caching
N-dimensional cubic spline interpolation between discrete points.
Definition: any.hpp:35
STL namespace.
const ext::shared_ptr< FdmInnerValueCalculator > calculator
const ext::shared_ptr< FdmMesher > mesher
std::uint64_t x_