1 /* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2 
3 /*
4  Copyright (C) 2007 François du Vignaud
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/math/matrixutilities/tapcorrelations.hpp>
21 #include <cmath>
22 
23 namespace QuantLib {
24 
triangularAnglesParametrization(const Array & angles,Size matrixSize,Size rank)25     Disposable<Matrix> triangularAnglesParametrization(const Array& angles,
26                                                        Size matrixSize,
27                                                        Size rank) {
28 
29         // what if rank == 1?
30         QL_REQUIRE((rank-1) * (2*matrixSize - rank) == 2*angles.size(),
31                     "rank-1) * (matrixSize - rank/2) == angles.size()");
32         Matrix m(matrixSize, matrixSize);
33 
34         // first row filling
35         m[0][0] = 1.0;
36         for (Size j=1; j<matrixSize; ++j)
37             m[0][j] = 0.0;
38 
39         // next ones...
40         Size k = 0; //angles index
41         for (Size i=1; i<m.rows(); ++i) {
42             Real sinProduct = 1.0;
43             Size bound = std::min(i,rank-1);
44             for (Size j=0; j<bound; ++j) {
45                 m[i][j] = std::cos(angles[k]);
46                 m[i][j] *= sinProduct;
47                 sinProduct *= std::sin(angles[k]);
48                 ++k;
49             }
50             m[i][bound] = sinProduct;
51             for (Size j=bound+1; j<m.rows(); ++j)
52                 m[i][j] = 0;
53         }
54         return m;
55     }
56 
lmmTriangularAnglesParametrization(const Array & angles,Size matrixSize,Size)57     Disposable<Matrix> lmmTriangularAnglesParametrization(const Array& angles,
58                                                           Size matrixSize,
59                                                           Size) {
60         Matrix m(matrixSize, matrixSize);
61         for (Size i=0; i<m.rows(); ++i) {
62             Real cosPhi, sinPhi;
63             if (i>0) {
64                 cosPhi = std::cos(angles[i-1]);
65                 sinPhi = std::sin(angles[i-1]);
66             } else {
67                 cosPhi = 1.0;
68                 sinPhi = 0.0;
69             }
70 
71             for (Size j=0; j<i; ++j)
72                 m[i][j] = sinPhi * m[i-1][j];
73 
74             m[i][i] = cosPhi;
75 
76             for (Size j=i+1; j<m.rows(); ++j)
77                 m[i][j] = 0.0;
78         }
79         return m;
80     }
81 
triangularAnglesParametrizationUnconstrained(const Array & x,Size matrixSize,Size rank)82     Disposable<Matrix> triangularAnglesParametrizationUnconstrained(
83                                                             const Array& x,
84                                                             Size matrixSize,
85                                                             Size rank) {
86         Array angles(x.size());
87         //we convert the unconstrained parameters in angles
88         for (Size i = 0; i < x.size(); ++i)
89             angles[i] = M_PI*.5 - std::atan(x[i]);
90         return triangularAnglesParametrization(angles, matrixSize, rank);
91     }
92 
lmmTriangularAnglesParametrizationUnconstrained(const Array & x,Size matrixSize,Size rank)93     Disposable<Matrix> lmmTriangularAnglesParametrizationUnconstrained(
94                                                             const Array& x,
95                                                             Size matrixSize,
96                                                             Size rank) {
97         Array angles(x.size());
98         //we convert the unconstrained parameters in angles
99         for (Size i = 0; i < x.size(); ++i)
100             angles[i] = M_PI*.5 - std::atan(x[i]);
101         return lmmTriangularAnglesParametrization(angles, matrixSize, rank);
102     }
103 
triangularAnglesParametrizationRankThree(Real alpha,Real t0,Real epsilon,Size matrixSize)104     Disposable<Matrix> triangularAnglesParametrizationRankThree(
105                                             Real alpha, Real t0,
106                                             Real epsilon, Size matrixSize) {
107         Matrix m(matrixSize, 3);
108         for (Size i=0; i<m.rows(); ++i) {
109             Real t = t0 * (1 - std::exp(epsilon*Real(i)));
110             Real phi = std::atan(alpha * t);
111             m[i][0] = std::cos(t)*std::cos(phi);
112             m[i][1] = std::sin(t)*std::cos(phi);
113             m[i][2] = -std::sin(phi);
114         }
115         return m;
116     }
117 
triangularAnglesParametrizationRankThreeVectorial(const Array & parameters,Size nbRows)118     Disposable<Matrix> triangularAnglesParametrizationRankThreeVectorial(
119                                                     const Array& parameters,
120                                                     Size nbRows) {
121         QL_REQUIRE(parameters.size() == 3,
122                    "the parameter array must contain exactly 3 values" );
123         return triangularAnglesParametrizationRankThree(parameters[0],
124                                                         parameters[1],
125                                                         parameters[2],
126                                                         nbRows);
127 
128     }
129 
value(const Array & x) const130     Real FrobeniusCostFunction::value(const Array& x) const {
131         Array temp = values(x);
132         return DotProduct(temp, temp);
133     }
134 
values(const Array & x) const135     Disposable<Array> FrobeniusCostFunction::values(const Array& x) const {
136         Array result((target_.rows()*(target_.columns()-1))/2);
137         Matrix pseudoRoot = f_(x, matrixSize_, rank_);
138         Matrix differences = pseudoRoot * transpose(pseudoRoot) - target_;
139         Size k = 0;
140         // then we store the elementwise differences in a vector.
141         for (Size i=0; i<target_.rows(); ++i) {
142             for (Size j=0; j<i; ++j){
143                 result[k] = differences[i][j];
144                 ++k;
145             }
146         }
147         return result;
148     }
149 }
150