1function [DP6,DP6inv] = Q6_plication(p) 2% Computes the 6-way duplication Matrix DP6 (and its Moore-Penrose inverse) 3% such that for any p-dimensional vector x: 4% y=kron(kron(kron(kron(kron(x,x),x,x),x),x)=DP6*z 5% where z is of dimension np=p*(p+1)*(p+2)*(p+3)*(p+4)*(p+5)/(1*2*3*4*5*6) 6% and is obtained from y by removing each second and later occurence of the 7% same element. This is a generalization of the Duplication matrix. 8% Reference: Meijer (2005) - Matrix algebra for higher order moments. 9% Linear Algebra and its Applications, 410,pp. 112-134 10% ========================================================================= 11% INPUTS 12% * p [integer] size of vector 13% ------------------------------------------------------------------------- 14% OUTPUTS 15% * DP6 [p^6 by np] 6-way duplication matrix 16% * DP6inv [np by np] Moore-Penrose inverse of DP6 17% ------------------------------------------------------------------------- 18% This function is called by 19% * pruned_state_space_system.m 20% ------------------------------------------------------------------------- 21% This function calls 22% * binom_coef (embedded) 23% * mue (embedded) 24% * uperm 25% ========================================================================= 26% Copyright (C) 2020 Dynare Team 27% 28% This file is part of Dynare. 29% 30% Dynare is free software: you can redistribute it and/or modify 31% it under the terms of the GNU General Public License as published by 32% the Free Software Foundation, either version 3 of the License, or 33% (at your option) any later version. 34% 35% Dynare is distributed in the hope that it will be useful, 36% but WITHOUT ANY WARRANTY; without even the implied warranty of 37% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 38% GNU General Public License for more details. 39% 40% You should have received a copy of the GNU General Public License 41% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 42% ========================================================================= 43np = p*(p+1)*(p+2)*(p+3)*(p+4)*(p+5)/(1*2*3*4*5*6); 44DP6 = spalloc(p^6,np,p^6); 45counti=1; 46for i1=1:p 47 for i2=i1:p 48 for i3=i2:p 49 for i4=i3:p 50 for i5=i4:p 51 for i6=i5:p 52 idx = uperm([i6 i5 i4 i3 i2 i1]); 53 for r = 1:size(idx,1) 54 ii1 = idx(r,1); ii2= idx(r,2); ii3=idx(r,3); ii4=idx(r,4); ii5=idx(r,5); ii6=idx(r,6); 55 n = ii1 + (ii2-1)*p + (ii3-1)*p^2 + (ii4-1)*p^3 + (ii5-1)*p^4 + (ii6-1)*p^5; 56 m = mue(p,i6,i5,i4,i3,i2,i1); 57 DP6(n,m)=1; 58 end 59 counti = counti+1; 60 end 61 end 62 end 63 end 64 end 65end 66DP6inv = (transpose(DP6)*DP6)\transpose(DP6); 67 68function m = mue(p,i1,i2,i3,i4,i5,i6) 69% Auxiliary expression, see page 122 of Meijer (2005) 70 m = binom_coef(p,6,1) - binom_coef(p,1,i1+1) - binom_coef(p,2,i2+1) - binom_coef(p,3,i3+1) - binom_coef(p,4,i4+1) - binom_coef(p,5,i5+1) - binom_coef(p,6,i6+1); 71 m = round(m); 72end 73 74function N = binom_coef(p,q,i) 75% Auxiliary expression for binomial coefficients, see page 119 of Meijer (2005) 76 t = q; r =p+q-i; 77 if t==0 78 N=1; 79 else 80 N=1; 81 for h = 0:(t-1) 82 N = N*(r-h); 83 end 84 N=N/factorial(t); 85 end 86end 87end