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 // Moments of normal distribution. 22 23 /* Here we calculate the higher order moments of normally distributed 24 random vector u with means equal to zero and given 25 variance-covariance matrix V, i.e. u↝��(0,V). The moment 26 generating function for such distribution is f(t)=e^½tᵀVt. If 27 we derivate it w.r.t. t and unfold the higher dimensional tensors 28 row-wise, we obtain terms like: 29 30 ∂f(t)/∂t = f(t)·Vt 31 ∂²f(t)/∂t² = f(t)·(Vt⊗v) 32 ∂³f(t)/∂t³ = f(t)·(Vt⊗Vt⊗Vt + P_?(v⊗Vt) + P_?(Vt⊗v) + v⊗Vt) 33 ∂⁴f(t)/∂t⁴ = f(t)·(Vt⊗Vt⊗Vt⊗Vt + S_?(v⊗Vt⊗Vt) + S_?(Vt⊗v⊗Vt) + S_?(Vt⊗Vt⊗v) + S_?(v⊗v)) 34 35 where v is vectorized V (v=vec(V)), and P_? is a suitable row permutation 36 (corresponds to permutation of multidimensional indices) which permutes the 37 tensor data, so that the index of a variable being derived would be the 38 last. This ensures that all (permuted) tensors can be summed yielding a 39 tensor whose indices have some order (in here we chose the order that more 40 recent derivating variables are to the right). Finally, S_? is a suitable sum 41 of various P_?. 42 43 We are interested in S_? multiplying the Kronecker powers ⊗ⁿv. The S_? is a 44 (possibly) multi-set of permutations of even order. Note that we know the 45 number of permutations in S_?. The above formulas for f(t) derivatives are 46 valid also for monomial u, and from literature we know that 2n-th moment is 47 (2n!)/(n!2ⁿ)·σ². So there are (2n!)/(n!2ⁿ) permutations in S_?. 48 49 In order to find the S_? permutation we need to define a couple of things. 50 First we define a sort of equivalence between the permutations applicable to 51 even number of indices. We write P₁ ≡ P₂ whenever P₁⁻¹∘P₂ permutes only whole 52 pairs, or items within pairs, but not indices across the pairs. For instance 53 the permutations (0,1,2,3) and (3,2,0,1) are equivalent, but (0,2,1,3) is 54 not equivalent with the two. Clearly, the ≡ relationship is an equivalence. 55 56 This allows to define a relation ⊑ between the permutation multi-sets S, 57 which is basically the subset relation ⊆ but with respect to the equivalence 58 relation ≡, more formally: 59 60 S₁ ⊑ S₂ iff P ∈ S₁ ⇒ ∃Q ∈ S₂ : P ≡ Q 61 62 This induces an equivalence S₁ ≡ S₂. 63 64 Now let Fₙ denote a set of permutations on 2n indices which is maximal with 65 respect to ⊑, and minimal with respect to ≡ (in other words, it contains 66 everything up to the equivalence ≡). It is straightforward to calculate a 67 number of permutations in Fₙ. This is a total number of all permutations 68 of 2n divided by permutations of pairs divided by permutations within the 69 pairs. This is (2n!)/(n!2ⁿ). 70 71 We now prove that S_? ≡ Fₙ. Clearly S_? ⊑ Fₙ, since Fₙ is maximal. In order 72 to prove that Fₙ ⊑ S_?, let us assert that for any permutation P and for 73 any (semi-)positive definite matrix V we have 74 P·S_?(⊗ⁿv)=S_?(⊗ⁿv). Below we show that there is a positive 75 definite matrix V of some dimension that for any two permutation 76 multi-sets S₁, S₂, we have 77 78 S₁ ≢ S₂ ⇒ S₁(⊗ⁿv) ≠ S₂(⊗ⁿv) 79 80 So it follows that for any permutation P, we have P·S_? ≡ S_?. For a purpose 81 of contradiction let P ∈ Fₙ be a permutation which is not equivalent to any 82 permutation from S_?. Since S_? is non-empty, let us pick P₀ ∈ S_?. Now 83 assert that P₀⁻¹S_? ≢ P⁻¹S_? since the first contains an identity and the 84 second does not contain a permutation equivalent to identity. Thus we have 85 (P∘P₀⁻¹)S_? ≢ S_? which gives the contradiction and we have proved that 86 Fₙ ⊑ S_?. Thus Fₙ ≡ S_?. Moreover, we know that S_? and Fₙ have the same 87 number of permutations, hence the minimality of S_? with respect to ≡. 88 89 Now it suffices to prove that there exists a positive definite V such that 90 for any two permutation multi-sets S₁ and S₂ holds S₁ ≢ S₂ ⇒ S₁(⊗ⁿv) ≠ S₂⊗ⁿv. 91 If V is a n×n matrix, then S₁ ≢ S₂ implies that there is 92 identically nonzero polynomial of elements from V of order n over 93 integers. If V=AᵀA then there is identically non-zero polynomial of 94 elements from A of order 2n. This means, that we have to find n(n+1)/2 95 tuple x of real numbers such that all identically non-zero polynomials p 96 of order 2n over integers yield p(x)≠0. 97 98 The x is constructed as follows: x_i = π^log(rᵢ), where rᵢ is i-th prime. 99 Let us consider the monom x₁^j₁·…·xₖ^jₖ. When the monom is evaluated, we get 100 101 π^{log(r₁^j₁)+…+log(rₖ^jₖ)}= π^log(r₁^j₁+…+rₖ^jₖ) 102 103 Now it is easy to see that if an integer combination of such terms is zero, 104 then the combination must be either trivial or sum to 0 and all monoms 105 must be equal. Both cases imply a polynomial identically equal to zero. So, 106 any non-trivial integer polynomial evaluated at x must be non-zero. 107 108 So, having this result in hand, now it is straightforward to calculate 109 higher moments of normal distribution. Here we define a container, which 110 does the job. In its constructor, we simply calculate Kronecker powers of v 111 and apply Fₙ to ⊗ⁿv. Fₙ is, in fact, a set of all equivalences in sense of 112 class Equivalence over 2n elements, having n classes each of them having 113 exactly 2 elements. 114 */ 115 116 #ifndef NORMAL_MOMENTS_H 117 #define NORMAL_MOMENTS_H 118 119 #include "t_container.hh" 120 121 class UNormalMoments : public TensorContainer<URSingleTensor> 122 { 123 public: 124 UNormalMoments(int maxdim, const TwoDMatrix &v); 125 private: 126 void generateMoments(int maxdim, const TwoDMatrix &v); 127 static bool selectEquiv(const Equivalence &e); 128 }; 129 130 class FNormalMoments : public TensorContainer<FRSingleTensor> 131 { 132 public: 133 FNormalMoments(const UNormalMoments &moms); 134 }; 135 136 #endif 137