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