1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
|
// Copyright (C) 2016, 2017 EDF
// All Rights Reserved
// This code is published under the GNU Lesser General Public License (GNU LGPL)
#ifndef LOCALLINEARREGRESSION_H
#define LOCALLINEARREGRESSION_H
#include <vector>
#include <memory>
#include <array>
#include <Eigen/Dense>
#include "StOpt/regression/LocalRegression.h"
#include "StOpt/core/grids/InterpolatorSpectral.h"
/** \file LocalLinearRegression.h
* \brief Compute conditional expectations by using local regressions.
* See article "Monte-Carlo valorisation of American options: facts and new algorithms to improve existing methods"
* by Bouchard, Warin in "Numerical methods in finance", Springer,2012
* \author Xavier Warin
*/
namespace StOpt
{
/**
* \defgroup LocalLinear Piecewise linear regression
* \brief It implements local linear functions for performing local regression
* using P-1 finite elements.
*@{
*/
/// \class LocalLinearRegression LocalLinearRegression.h
/// To be used in Monte Carlo methods. Linear regression on each cell which is constructed such that each cell has
/// roughly the same number of particles
class LocalLinearRegression : public LocalRegression
{
protected :
Eigen::ArrayXXd m_matReg ; ///< Regression matrix (rank \f$ (n_d+1)^2 \f$ by the number of meshes ). Also store lower part of factorized matrix.
Eigen::ArrayXXd m_diagReg ; ///< Diagonal part of the factorized matrix during Choleski factorization
public :
/// \brief Default constructor
LocalLinearRegression() {}
/// \brief Constructor
/// \param p_nbMesh discretization in each direction
/// \param p_bRotationAndRecale do we use SVD
LocalLinearRegression(const Eigen::ArrayXi &p_nbMesh, bool p_bRotationAndRecale = false);
/// \brief Constructor for object constructed at each time step
/// \param p_bZeroDate first date is 0?
/// \param p_particles particles used for the meshes.
/// First dimension : dimension of the problem,
/// second dimension : the number of particles
/// \param p_nbMesh discretization in each direction
/// \param p_bRotationAndRecale do we use SVD
LocalLinearRegression(const bool &p_bZeroDate,
const Eigen::ArrayXXd &p_particles,
const Eigen::ArrayXi &p_nbMesh,
bool p_bRotationAndRecale = false);
/// \brief Second constructor , only to be used in simulation
/// \param p_bRotationAndRecale do we use SVD
LocalLinearRegression(const bool &p_bZeroDate,
const Eigen::ArrayXi &p_nbMesh,
const Eigen::Array< std::array< double, 2>, Eigen::Dynamic, Eigen::Dynamic > &p_mesh,
const std::vector< std::shared_ptr< Eigen::ArrayXd > > &p_mesh1D, const Eigen::ArrayXd &p_meanX,
const Eigen::ArrayXd &p_etypX, const Eigen::MatrixXd &p_svdMatrix,
const bool &p_bRotationAndRecale) : LocalRegression(p_bZeroDate, p_nbMesh, p_mesh, p_mesh1D, p_meanX, p_etypX, p_svdMatrix, p_bRotationAndRecale)
{
}
/// \brief Copy constructor
/// \param p_objet object to copy
LocalLinearRegression(const LocalLinearRegression &p_objet);
/// \brief update the particles used in regression and construct the matrices
/// \param p_bZeroDate first date is 0?
/// \param p_particles particles used for the meshes.
/// First dimension : dimension of the problem,
/// second dimension : the number of particles
void updateSimulations(const bool &p_bZeroDate, const Eigen::ArrayXXd &p_particles);
/// \brief Get some local accessors
///@{
inline const Eigen::ArrayXXd &getMatReg() const
{
return m_matReg;
}
inline const Eigen::ArrayXXd &getDiagReg() const
{
return m_diagReg;
}
///@}
/// \brief conditional expectation basis function coefficient calculation
/// \param p_fToRegress function to regress associated to each simulation used in optimization
/// \return regression coordinates on the basis (size : number of meshes multiplied by the dimension plus one)
/// @{
Eigen::ArrayXd getCoordBasisFunction(const Eigen::ArrayXd &p_fToRegress) const;
///@}
/// \brief conditional expectation basis function coefficient calculation for multiple functions to regress
/// \param p_fToRegress function to regress associated to each simulation used in optimization (size : number of functions to regress \times the number of Monte Carlo simulations)
/// \return regression coordinates on the basis (size : number of function to regress \times number of meshes multiplied by the dimension plus one)
/// @{
Eigen::ArrayXXd getCoordBasisFunctionMultiple(const Eigen::ArrayXXd &p_fToRegress) const ;
///@}
/// \brief conditional expectation calculation
/// \param p_fToRegress simulations to regress used in optimization
/// \return regressed value function
/// @{
Eigen::ArrayXd getAllSimulations(const Eigen::ArrayXd &p_fToRegress) const ;
Eigen::ArrayXXd getAllSimulationsMultiple(const Eigen::ArrayXXd &p_fToRegress) const;
///@}
/// \brief modify cell to fit data and get a convex mapping for a function to regress
/// Use "Convex piecewise-linear fitting" by Alessandro Magnani ยท Stephen P. Boyd, Optim Eng (2009)
/// http://web.stanford.edu/~boyd/papers/pdf/cvx_pwl_fit.pdf
/// The number of iterations should be small (< 10) as the methodology improves the convexity but doesn't give a full convex function
/// \param p_fToRegress function to regress associated to each simulation used in optimization
/// \param p_nbIterMax maximal number of iterations
/// Return the regression values
Eigen::ArrayXd getAllSimulationsConvex(const Eigen::ArrayXd &p_fToRegress, const int &p_nbIterMax) ;
/// \brief Use basis functions to reconstruct the solution
/// \param p_basisCoefficients basis coefficients
///@{
Eigen::ArrayXd reconstruction(const Eigen::ArrayXd &p_basisCoefficients) const;
Eigen::ArrayXXd reconstructionMultiple(const Eigen::ArrayXXd &p_basisCoefficients) const;
/// @}
/// \brief use basis function to reconstruct a given simulation
/// \param p_isim simulation number
/// \param p_basisCoefficients basis coefficients to reconstruct a given conditional expectation
double reconstructionASim(const int &p_isim, const Eigen::ArrayXd &p_basisCoefficients) const ;
/// \brief conditional expectation reconstruction
/// \param p_coordinates coordinates to interpolate (uncertainty sample)
/// \param p_coordBasisFunction regression coordinates on the basis (size: number of meshes multiplied by the dimension plus one)
/// \return regressed value function reconstructed for each simulation
double getValue(const Eigen::ArrayXd &p_coordinates,
const Eigen::ArrayXd &p_coordBasisFunction) const;
/// \brief permits to reconstruct a function with basis functions coefficients values given on a grid
/// \param p_coordinates coordinates (uncertainty sample)
/// \param p_ptOfStock grid point
/// \param p_interpFuncBasis spectral interpolator to interpolate the basis functions coefficients used in regression on the grid (given for each basis function)
double getAValue(const Eigen::ArrayXd &p_coordinates, const Eigen::ArrayXd &p_ptOfStock,
const std::vector< std::shared_ptr<InterpolatorSpectral> > &p_interpFuncBasis) const;
/// \brief get the number of basis functions
inline int getNumberOfFunction() const
{
if (m_bZeroDate)
return 1;
else
return m_nbMesh.prod() * (m_nbMesh.size() + 1);
}
/// \brief Clone the regressor
std::shared_ptr<BaseRegression> clone() const
{
return std::static_pointer_cast<BaseRegression>(std::make_shared<LocalLinearRegression>(*this));
}
/// \brief conditional expectation basis function coefficient calculation for a special cell
/// \param p_iCell cell number
/// \param p_fToRegress function to regress associated to each simulation used in optimization and corresponding to the cell
/// \return regression coordinates on the basis (size : the dimension of the problem plus one)
/// @{
Eigen::ArrayXd getCoordBasisFunctionOneCell(const int &p_iCell, const Eigen::ArrayXd &p_fToRegress) const ;
Eigen::ArrayXXd getCoordBasisFunctionMultipleOneCell(const int &p_iCell, const Eigen::ArrayXXd &p_fToRegress) const ;
///@}
/// \brief Given a particle and the coordinates of the mesh it belongs to, get back the conditional expectation
/// \param p_oneParticle The particle generated
/// \param p_cell the cell it belongs to
/// \param p_foncBasisCoef function basis on the current cell (the row is the number of reconstruction to achieve, the columns the number of function basis)
/// \return send back the array of conditional expectations
Eigen::ArrayXd getValuesOneCell(const Eigen::ArrayXd &p_oneParticle, const int &p_cell, const Eigen::ArrayXXd &p_foncBasisCoef) const;
};
/**@}*/
}
#endif /*LOCALLINEARREGRESSION_H*/
|
