1function [QP,QPinv] = quadruplication(p) 2% Computes the Quadruplication Matrix QP (and its Moore-Penrose inverse) 3% such that for any p-dimensional vector x: 4% y=kron(kron(kron(x,x),x),x)=QP*z 5% where z is of dimension np=p*(p+1)*(p+2)*(p+3)/2 and is obtained from y 6% by removing each second and later occurence of the same element. 7% 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% * QP [p^4 by np] Quadruplication matrix 16% * QPinv [np by np] Moore-Penrose inverse of QP 17% ------------------------------------------------------------------------- 18% This function is called by 19% * pruned_state_space_system.m 20% ------------------------------------------------------------------------- 21% This function calls 22% * mue (embedded) 23% * uperm 24% ========================================================================= 25% Copyright (C) 2020 Dynare Team 26% 27% This file is part of Dynare. 28% 29% Dynare is free software: you can redistribute it and/or modify 30% it under the terms of the GNU General Public License as published by 31% the Free Software Foundation, either version 3 of the License, or 32% (at your option) any later version. 33% 34% Dynare is distributed in the hope that it will be useful, 35% but WITHOUT ANY WARRANTY; without even the implied warranty of 36% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 37% GNU General Public License for more details. 38% 39% You should have received a copy of the GNU General Public License 40% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 41% ========================================================================= 42np = p*(p+1)*(p+2)*(p+3)/24; 43QP = spalloc(p^4,np,p^4); 44if nargout > 1 45 QPinv = spalloc(np,np,p^4); 46end 47counti = 1; 48for l=1:p 49 for k=l:p 50 for j=k:p 51 for i=j:p 52 idx = uperm([i j k l]); 53 for r = 1:size(idx,1) 54 ii = idx(r,1); jj= idx(r,2); kk=idx(r,3); ll=idx(r,4); 55 n = ii + (jj-1)*p + (kk-1)*p^2 + (ll-1)*p^3; 56 m = mue(p,i,j,k,l); 57 QP(n,m)=1; 58 if nargout > 1 59 if i==j && j==k && k==l 60 QPinv(m,n)=1; 61 elseif i==j && j==k && k>l 62 QPinv(m,n)=1/4; 63 elseif i>j && j==k && k==l 64 QPinv(m,n)=1/4; 65 elseif i==j && j>k && k==l 66 QPinv(m,n) = 1/6; 67 elseif i>j && j>k && k==l 68 QPinv(m,n) = 1/12; 69 elseif i>j && j==k && k>l 70 QPinv(m,n) = 1/12; 71 elseif i==j && j>k && k>l 72 QPinv(m,n) = 1/12; 73 elseif i>j && j>k && k>l 74 QPinv(m,n) = 1/24; 75 end 76 end 77 end 78 counti = counti+1; 79 end 80 end 81 end 82end 83%QPinv = (transpose(QP)*QP)\transpose(QP); 84 85function m = mue(p,i,j,k,l) 86 % Auxiliary expression, see page 118 of Meijer (2005) 87 m = i + (j-1)*p + 1/2*(k-1)*p^2 + 1/6*(l-1)*p^3 - 1/2*j*(j-1) + 1/6*k*(k-1)*(k-2) - 1/24*l*(l-1)*(l-2)*(l-3) - 1/2*(k-1)^2*p + 1/6*(l-1)^3*p - 1/4*(l-1)*(l-2)*p^2 - 1/4*l*(l-1)*p + 1/6*(l-1)*p; 88 m = round(m); 89end 90 91 92end