1% 2% prodmom.m Date: 4/29/2006 3% This Matlab program computes the product moment of X_{i_1}^{nu_1}X_{i_2}^{nu_2}...X_{i_m}^{nu_m}, 4% where X_{i_j} are elements from X ~ N(0_n,V). 5% V only needs to be positive semidefinite. 6% V: variance-covariance matrix of X 7% ii: vector of i_j 8% nu: power of X_{i_j} 9% Reference: Triantafyllopoulos (2003) On the Central Moments of the Multidimensional 10% Gaussian Distribution, Mathematical Scientist 11% Kotz, Balakrishnan, and Johnson (2000), Continuous Multivariate 12% Distributions, Vol. 1, p.261 13% Note that there is a typo in Eq.(46.25), there should be an extra rho in front 14% of the equation. 15% Usage: prodmom(V,[i1 i2 ... ir],[nu1 nu2 ... nur]) 16% Example: To get E[X_2X_4^3X_7^2], use prodmom(V,[2 4 7],[1 3 2]) 17% 18% Retrieved from http://www-2.rotman.utoronto.ca/~kan/papers/prodmom.zip 19% This function is part of replication codes of the following paper: 20% Kan, R.: "From moments of sum to moments of product." Journal of 21% Multivariate Analysis, 2008, vol. 99, issue 3, pages 542-554. 22% ========================================================================= 23% Copyright (C) 2008-2015 Raymond Kan <kan@chass.utoronto.ca> 24% Copyright (C) 2019-2020 Dynare Team 25% 26% This file is part of Dynare. 27% 28% Dynare is free software: you can redistribute it and/or modify 29% it under the terms of the GNU General Public License as published by 30% the Free Software Foundation, either version 3 of the License, or 31% (at your option) any later version. 32% 33% Dynare is distributed in the hope that it will be useful, 34% but WITHOUT ANY WARRANTY; without even the implied warranty of 35% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 36% GNU General Public License for more details. 37% 38% You should have received a copy of the GNU General Public License 39% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 40% ========================================================================= 41function y = prodmom(V,ii,nu); 42if nargin<3 43 nu = ones(size(ii)); 44end 45s = sum(nu); 46if s==0 47 y = 1; 48 return 49end 50if rem(s,2)==1 51 y = 0; 52 return 53end 54nuz = nu==0; 55nu(nuz) = []; 56ii(nuz) = []; 57m = length(ii); 58V = V(ii,ii); 59s2 = s/2; 60% 61% Use univariate normal results 62% 63if m==1 64 y = V^s2*prod([1:2:s-1]); 65 return 66end 67% 68% Use bivariate normal results when there are only two distinct indices 69% 70if m==2 71 rho = V(1,2)/sqrt(V(1,1)*V(2,2)); 72 y = V(1,1)^(nu(1)/2)*V(2,2)^(nu(2)/2)*bivmom(nu,rho); 73 return 74end 75% 76% Regular case 77% 78[nu,inu] = sort(nu,2,'descend'); 79V = V(inu,inu); % Extract only the relevant part of V 80x = zeros(1,m); 81V = V./2; 82nu2 = nu./2; 83p = 2; 84q = nu2*V*nu2'; 85y = 0; 86for i=1:fix(prod(nu+1)/2) 87 y = y+p*q^s2; 88 for j=1:m 89 if x(j)<nu(j) 90 x(j) = x(j)+1; 91 p = -round(p*(nu(j)+1-x(j))/x(j)); 92 q = q-2*(nu2-x)*V(:,j)-V(j,j); 93 break 94 else 95 x(j) = 0; 96 if rem(nu(j),2)==1 97 p = -p; 98 end 99 q = q+2*nu(j)*(nu2-x)*V(:,j)-nu(j)^2*V(j,j); 100 end 101 end 102end 103y = y/prod([1:s2]);