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
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
|
// Copyright (C) 2016 EDF
// All Rights Reserved
// This code is published under the GNU Lesser General Public License (GNU LGPL)
#ifndef SPARSEREGRESSION_H
#define SPARSEREGRESSION_H
#include <Eigen/Dense>
#include <functional>
#ifdef _OPENMP
#include <omp.h>
#endif
#include "StOpt/core/grids/SparseSpaceGridNoBound.h"
#include "StOpt/regression/BaseRegression.h"
/** \file SparseRegression.h
* \brief Compute conditional expectation with sparse grids
* As in 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
* particle are sorted in each dimension such that the basis support have particles
* First draft is without adaptation.
* \author Xavier Warin
*/
namespace StOpt
{
/**
* \defgroup Regression with sparse grids with linear, quadratic and cubic polynomials
*@{
*/
class SparseRegression : public BaseRegression
{
private:
std::shared_ptr< SparseSpaceGridNoBound > m_spGrid ; ///< to store the grid
Eigen::ArrayXXd m_partRescaled; ///< Particles used to regress but rescaled according to transformToUnitSquare : dimension of the problem , second dimension : the number of particles
std::vector< std::vector< double> > m_mesh ; ///< For each direction, give the coordinate of the mesh
std::vector< std::vector< std::function< double(const double &) > > > m_functionScal ; ///< for each position, for the level used, defines the basis function
Eigen::LLT<Eigen::MatrixXd> m_llt; ///< llt factorization by eigen
Eigen::MatrixXd m_yReg ; ///< utilitarian to calculate pseudo inverse if the matrix is singular
bool m_bSingular ; ///< true if the regression matrix is singular
bool m_bNoRescale ; //< Don't rescale if true
/// \brief create and factorize the regression matrix
/// \param p_levelMax maximum level of the sparse grid
/// \param p_weight weight for the anisotropy : the level \f$ (l_i)_i\f$ satisfy \f$\sum_i weight[i] l_i \le NDIM + levelMax -1 \f$
/// \param p_degree Degree of the interpolation sparse grid
void createAndFactorize(const int &p_levelMax, const Eigen::ArrayXd &p_weight, const int &p_degree);
/// \brief rescale the particle values
/// \param p_aParticle The particle to rescale according the mesh
/// \param p_aPartRescaled The same particle rescaled
void rescaleAParticle(const Eigen::Ref< const Eigen::ArrayXd> p_aParticle, Eigen::Ref< Eigen::ArrayXd> p_aPartRescaled) const;
/// \brief Transform each direction to unit square : in each direction , the same number number of particles in each mesh.
/// \param p_levelMax maximum level of the sparse grid
/// \param p_weight weight for the anisotropy : the level \f$ (l_i)_i\f$ satisfy \f$\sum_i weight[i] l_i \le NDIM + levelMax -1 \f$
void transformToUnitSquare(const int &p_levelMax, const Eigen::ArrayXd &p_weight);
/// \brief Rescale simply
/// \param p_weight weight for the anisotropy : the level \f$ (l_i)_i\f$ satisfy \f$ \sum_i weight[i] l_i \le NDIM + levelMax -1 \f$
void rescale(const Eigen::ArrayXd &p_weight);
/// \brief Recursive calculation function basis at a point using "son"" knowledge
/// \param p_levelCurrent Current index of levels of the point
/// \param p_positionCurrent Current position of the point, at the given level
/// \param p_ipoint Point number
/// \param p_xMiddle Position in [0,1] of current node in each dimension
/// \param p_dx Semi mesh size
/// \param p_x Evaluation point
/// \param p_idimMin Minimal dimension search (to avoid to go twice at same node)
/// \param p_funcVal Function basis values at current node for all dimensions
/// \param p_son Son array (first dimension is the node number, second is the dimension , 0 in array corresponds to left, 1 to right)
/// \param p_nonNullFunctionValues To fill in : values of the different non null function basis
/// \param p_associatedFunctionNumber Cell number associated to each previous non num value
void recursiveFillFunctionBasisWithSon(Eigen::ArrayXc &p_levelCurrent,
Eigen::ArrayXui &p_positionCurrent,
const int &p_ipoint,
Eigen::ArrayXd &p_xMiddle,
Eigen::ArrayXd &p_dx,
const Eigen::ArrayXd &p_x,
const unsigned short int &p_idimMin,
Eigen::ArrayXd &p_funcVal,
const Eigen::Array< std::array<int, 2 >, Eigen::Dynamic, Eigen::Dynamic > &p_son,
std::vector<double> &p_nonNullFunctionValues,
std::vector<int> &p_associatedFunctionNumber) const;
/// \brief For a particle asses the value of all basis functions
/// \param p_particle particle rescaled coordinates
/// \param p_nonNullFunctionValues To fill in : values of the different non null function basis
/// \param p_associatedFunctionNumber Cell number associated to each previous non num value
void assessFuncBasis(const Eigen::ArrayXd &p_particle, std::vector<double> &p_nonNullFunctionValues, std::vector<int> &p_associatedFunctionNumber) const ;
/// \brief fill in matrix regression using the minimum of memory possible
/// Each thread fills part of the matrix for the simulations it uses, then the matrix is reconstructed
/// \param p_matReg matrix to fill in
void fillInRegressionMatrix(Eigen::MatrixXd &p_matReg);
/// \brief Fill in second member with minimal cost of storage
/// \param p_fToRegress functions to regress associated to each simulation used in optimization (number of functions to regress times the simulations)
/// \return second members to construct (number of functions to regress times the number of basis functions)
Eigen::MatrixXd fillInSecondMember(const Eigen::MatrixXd &p_fToRegress) const ;
/// \brief Reconstruct the solution using the basis coordinates
/// \param p_basisCoeff basis function coefficients (number)
/// \param p_fRegressed the result of the regression
void reconstruct(const Eigen::ArrayXXd &p_basisCoeff, Eigen::ArrayXXd &p_fRegressed) const ;
/// \brief Reconstruct the solution using the basis coordinates
/// \param p_basisCoeff basis function coefficients
/// \param p_aParticle coordinates of the point where to reconstruct the value function
/// \return regressed value
double reconstructOneParticle(const Eigen::ArrayXd &p_basisCoeff, const Eigen::ArrayXd &p_aParticle) const ;
public :
/// \brief default constructor
SparseRegression() {}
/// \brief Constructor (for constructor achieved once)
/// \param p_levelMax Level max associated to the sparse grid
/// \param p_weight Weight for anisotropic sparse grids
/// \param p_degree Degree of the interpolator
/// \param p_bNoRescale Don't use a complex rescaling
SparseRegression(const int &p_levelMax, const Eigen::ArrayXd &p_weight, const int &p_degree, bool p_bNoRescale = 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_levelMax Level max associated to the sparse grid
/// \param p_weight Weight for anisotropic sparse grids
/// \param p_degree Degree of the interpolator
SparseRegression(const bool &p_bZeroDate,
const Eigen::ArrayXXd &p_particles,
const int &p_levelMax, const Eigen::ArrayXd &p_weight,
const int &p_degree);
/// \brief constructor to recreate the dumped object (for simulation)
/// \param p_bZeroDate first date is 0?
/// \param p_spGrid sparse grid
/// \param p_mesh mesh used for rescaling
/// \param p_meanX scaled factor in each direction (average of particles values in each direction)
/// \param p_etypX scaled factor in each direction (standard deviation of particles in each direction)
/// \param p_svdMatrix svd matrix transposed used to transform particles
SparseRegression(const bool &p_bZeroDate, std::shared_ptr< SparseSpaceGridNoBound> p_spGrid, const std::vector< std::vector< double> > &p_mesh, const Eigen::ArrayXd &p_meanX,
const Eigen::ArrayXd &p_etypX, const Eigen::MatrixXd &p_svdMatrix);
/// \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 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;
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 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 Accessor
///@{
std::shared_ptr<SparseSpaceGridNoBound> getSpGrid() const
{
return m_spGrid ; ///< grid
}
const std::vector< std::vector< double> > &getMesh() const
{
return m_mesh; ///< mesh to rescale particles;
}
///@}
/// \brief get the number of basis functions
inline int getNumberOfFunction() const
{
return m_spGrid->getNbPoints();
}
/// \brief Clone the regressor
std::shared_ptr<BaseRegression> clone() const
{
return std::static_pointer_cast<BaseRegression>(std::make_shared<SparseRegression>(*this));
}
};
/**@}*/
}
#endif /* SPARSEREGRESSION */
|
