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