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