1 /*
2 * Copyright © 2004 Ondra Kamenik
3 * Copyright © 2019 Dynare Team
4 *
5 * This file is part of Dynare.
6 *
7 * Dynare is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * Dynare is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with Dynare. If not, see <http://www.gnu.org/licenses/>.
19 */
20
21 #include <memory>
22
23 #include "normal_moments.hh"
24 #include "permutation.hh"
25 #include "kron_prod.hh"
26 #include "tl_static.hh"
27
UNormalMoments(int maxdim,const TwoDMatrix & v)28 UNormalMoments::UNormalMoments(int maxdim, const TwoDMatrix &v)
29 : TensorContainer<URSingleTensor>(1)
30 {
31 if (maxdim >= 2)
32 generateMoments(maxdim, v);
33 }
34
35 /* Here we fill up the container with the tensors for d=2,4,6,… up to the given
36 dimension. Each tensor of moments is equal to Fₙ(⊗ⁿv). This has a dimension
37 equal to 2n. See the header file for proof and details.
38
39 Here we sequentially construct the Kronecker power ⊗ⁿv and apply Fₙ.
40 */
41 void
generateMoments(int maxdim,const TwoDMatrix & v)42 UNormalMoments::generateMoments(int maxdim, const TwoDMatrix &v)
43 {
44 TL_RAISE_IF(v.nrows() != v.ncols(),
45 "Variance-covariance matrix is not square in UNormalMoments constructor");
46
47 int nv = v.nrows();
48 auto mom2 = std::make_unique<URSingleTensor>(nv, 2);
49 mom2->getData() = v.getData();
50 insert(std::move(mom2));
51 auto kronv = std::make_unique<URSingleTensor>(nv, 2);
52 kronv->getData() = v.getData();
53 for (int d = 4; d <= maxdim; d += 2)
54 {
55 auto newkronv = std::make_unique<URSingleTensor>(nv, d);
56 KronProd::kronMult(ConstVector(v.getData()),
57 ConstVector(kronv->getData()),
58 newkronv->getData());
59 kronv = std::move(newkronv);
60 auto mom = std::make_unique<URSingleTensor>(nv, d);
61 // apply Fₙ to ‘kronv’
62 /* Here we go through all equivalences, select only those having 2
63 elements in each class, then go through all elements in ‘kronv’ and
64 add to permuted location of ‘mom’.
65
66 The permutation must be taken as inverse of the permutation implied by
67 the equivalence, since we need a permutation which after application
68 to identity of indices yileds indices in the equivalence classes. Note
69 how the Equivalence::apply() method works. */
70 mom->zeros();
71 const EquivalenceSet eset = TLStatic::getEquiv(d);
72 for (const auto &cit : eset)
73 if (selectEquiv(cit))
74 {
75 Permutation per(cit);
76 per.inverse();
77 for (Tensor::index it = kronv->begin(); it != kronv->end(); ++it)
78 {
79 IntSequence ind(kronv->dimen());
80 per.apply(it.getCoor(), ind);
81 Tensor::index it2(*mom, ind);
82 mom->get(*it2, 0) += kronv->get(*it, 0);
83 }
84 }
85 insert(std::move(mom));
86 }
87 }
88
89 /* We return true for an equivalence whose each class has 2 elements. */
90
91 bool
selectEquiv(const Equivalence & e)92 UNormalMoments::selectEquiv(const Equivalence &e)
93 {
94 if (2*e.numClasses() != e.getN())
95 return false;
96 for (const auto &si : e)
97 if (si.length() != 2)
98 return false;
99 return true;
100 }
101
102 /* Here we go through all the unfolded container, fold each tensor and
103 insert it. */
FNormalMoments(const UNormalMoments & moms)104 FNormalMoments::FNormalMoments(const UNormalMoments &moms)
105 : TensorContainer<FRSingleTensor>(1)
106 {
107 for (const auto &mom : moms)
108 insert(std::make_unique<FRSingleTensor>(*(mom.second)));
109 }
110