1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2006 Mark Joshi
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/methods/montecarlo/genericlsregression.hpp>
21#include <ql/math/statistics/sequencestatistics.hpp>
22#include <ql/math/matrixutilities/svd.hpp>
23
24namespace QuantLib {
25
26 Real genericLongstaffSchwartzRegression(
27 std::vector<std::vector<NodeData> >& simulationData,
28 std::vector<std::vector<Real> >& basisCoefficients) {
29
30 Size steps = simulationData.size();
31 basisCoefficients.resize(new_size: steps-1);
32
33 for (Size i=steps-1; i!=0; --i) {
34
35 std::vector<NodeData>& exerciseData = simulationData[i];
36
37 // 1) find the covariance matrix of basis function values and
38 // deflated cash-flows
39 Size N = exerciseData.front().values.size();
40 std::vector<Real> temp(N+1);
41 SequenceStatistics stats(N+1);
42
43 Size j;
44 for (j=0; j<exerciseData.size(); ++j) {
45 if (exerciseData[j].isValid) {
46 std::copy(first: exerciseData[j].values.begin(),
47 last: exerciseData[j].values.end(),
48 result: temp.begin());
49 temp.back() = exerciseData[j].cumulatedCashFlows
50 - exerciseData[j].controlValue;
51
52 stats.add(sample: temp);
53 }
54 }
55
56 std::vector<Real> means = stats.mean();
57 Matrix covariance = stats.covariance();
58
59 Matrix C(N,N);
60 Array target(N);
61 for (Size k=0; k<N; ++k) {
62 target[k] = covariance[k][N] + means[k]*means[N];
63 for (Size l=0; l<=k; ++l)
64 C[k][l] = C[l][k] = covariance[k][l] + means[k]*means[l];
65 }
66
67 // 2) solve for least squares regression
68 Array alphas = SVD(C).solveFor(target);
69 basisCoefficients[i-1].resize(new_size: N);
70 std::copy(first: alphas.begin(), last: alphas.end(),
71 result: basisCoefficients[i-1].begin());
72
73 // 3) use exercise strategy to divide paths into exercise and
74 // non-exercise domains
75 for (j=0; j<exerciseData.size(); ++j) {
76 if (exerciseData[j].isValid) {
77 Real exerciseValue = exerciseData[j].exerciseValue;
78 Real continuationValue =
79 exerciseData[j].cumulatedCashFlows;
80 Real estimatedContinuationValue =
81 std::inner_product(
82 first1: exerciseData[j].values.begin(),
83 last1: exerciseData[j].values.end(),
84 first2: alphas.begin(),
85 init: exerciseData[j].controlValue);
86
87 // for exercise paths, add deflated rebate to
88 // deflated cash-flows at previous time frame;
89 // for non-exercise paths, add deflated cash-flows to
90 // deflated cash-flows at previous time frame
91 Real value = estimatedContinuationValue <= exerciseValue ?
92 exerciseValue :
93 continuationValue;
94
95 simulationData[i-1][j].cumulatedCashFlows += value;
96 }
97 }
98 }
99
100 // the value of the product can now be estimated by averaging
101 // over all paths
102 Statistics estimate;
103 std::vector<NodeData>& estimatedData = simulationData[0];
104 for (auto& j : estimatedData)
105 estimate.add(value: j.cumulatedCashFlows);
106
107 return estimate.mean();
108 }
109
110}
111
112

source code of quantlib/ql/methods/montecarlo/genericlsregression.cpp