1function [y,dy] = bivmom(p,rho) 2% Computes the product moment (and its derivative with respect to standard 3% errors and correlation parameters) of X_1^{p_1}X_2^{p_2}, where X_1 and X_2 4% are standard bivariate normally distributed. 5% n : dimension of X 6% rho: correlation coefficient between X_1 and X_2 7% ========================================================================= 8% INPUTS 9% p [2 by 1] powers of X_{1} and X_{2} 10% rho [1 by 1] correlation coefficient between X_1 and X_2 11% ------------------------------------------------------------------------- 12% OUTPUTS 13% y [1 by 1] product moment E[X_1^{p_1}X_2^{p_2}] 14% dy [1 by 1] derivative of y wrt to rho 15% ------------------------------------------------------------------------- 16% This function is based upon bivmom.m which is part of replication codes 17% of the following paper: 18% Kan, R.: "From moments of sum to moments of product." Journal of 19% Multivariate Analysis, 2008, vol. 99, issue 3, pages 542-554. 20% bivmom.m can be retrieved from http://www-2.rotman.utoronto.ca/~kan/papers/prodmom.zip 21% Further references: 22% Kotz, Balakrishnan, and Johnson (2000), Continuous Multivariate Distributions, Vol. 1, p.261 23% Note that there is a typo in Eq.(46.25), there should be an extra rho in front 24% of the equation. 25% ========================================================================= 26% Copyright (C) 2008-2015 Raymond Kan <kan@chass.utoronto.ca> 27% Copyright (C) 2019-2020 Dynare Team 28% 29% This file is part of Dynare. 30% 31% Dynare is free software: you can redistribute it and/or modify 32% it under the terms of the GNU General Public License as published by 33% the Free Software Foundation, either version 3 of the License, or 34% (at your option) any later version. 35% 36% Dynare is distributed in the hope that it will be useful, 37% but WITHOUT ANY WARRANTY; without even the implied warranty of 38% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 39% GNU General Public License for more details. 40% 41% You should have received a copy of the GNU General Public License 42% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 43% ========================================================================= 44s1 = p(1); 45s2 = p(2); 46rho2 = rho^2; 47if nargout > 1 48 drho2 = 2*rho; 49end 50if rem(s1+s2,2)==1 51 y = 0; 52 return 53end 54r = fix(s1/2); 55s = fix(s2/2); 56y = 1; 57c = 1; 58if nargout > 1 59 dy = 0; 60 dc = 0; 61end 62odd = 2*rem(s1,2); 63for j=1:min(r,s) 64 if nargout > 1 65 dc = 2*dc*(r+1-j)*(s+1-j)*rho2/(j*(2*j-1+odd)) + 2*c*(r+1-j)*(s+1-j)*drho2/(j*(2*j-1+odd)); 66 end 67 c = 2*c*(r+1-j)*(s+1-j)*rho2/(j*(2*j-1+odd)); 68 y = y+c; 69 if nargout > 1 70 dy = dy + dc; 71 end 72end 73if odd 74 if nargout > 1 75 dy = y + dy*rho; 76 end 77 y = y*rho; 78end 79y = prod([1:2:s1])*prod([1:2:s2])*y; 80if nargout > 1 81 dy = prod([1:2:s1])*prod([1:2:s2])*dy; 82end 83 84