1function [Dp,DpMPinv] = duplication(p)
2% [Dp,DpMPinv] = duplication(p)
3% -------------------------------------------------------------------------
4% Duplication Matrix as defined by Magnus and Neudecker (2002), p.49
5% =========================================================================
6% INPUTS
7%   p:  [integer] length of vector
8% -------------------------------------------------------------------------
9% OUTPUTS
10%   Dp:      Duplication matrix
11%   DpMPinv: Moore-Penroze inverse of Dp
12% -------------------------------------------------------------------------
13% This function is called by
14%   * get_identification_jacobians.m (previously getJJ.m)
15% =========================================================================
16% Copyright (C) 1997 Tom Minka <minka@microsoft.com>
17% Copyright (C) 2019 Dynare Team
18%
19% This file is part of Dynare.
20%
21% Dynare is free software: you can redistribute it and/or modify
22% it under the terms of the GNU General Public License as published by
23% the Free Software Foundation, either version 3 of the License, or
24% (at your option) any later version.
25%
26% Dynare is distributed in the hope that it will be useful,
27% but WITHOUT ANY WARRANTY; without even the implied warranty of
28% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
29% GNU General Public License for more details.
30%
31% You should have received a copy of the GNU General Public License
32% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
33% =========================================================================
34% Original author: Thomas P Minka (tpminka@media.mit.edu), April 22, 2013
35
36a = tril(ones(p));
37i = find(a);
38a(i) = 1:length(i);
39a = a + transpose(tril(a,-1));
40
41j = a(:);
42
43m = p*(p+1)/2;
44Dp = spalloc(p*p,m,p^2);
45for r = 1:size(Dp,1)
46    Dp(r, j(r)) = 1;
47end
48
49if nargout > 1
50    DpMPinv = (transpose(Dp)*Dp)\transpose(Dp);
51end
52