Commit df084f8e authored by Sébastien Villemot's avatar Sébastien Villemot

Merge branch 'identification_higher_order' into 'master'

Fix pruned state space at third order

See merge request !1711
parents f69cd149 62230166
Pipeline #3199 passed with stages
in 69 minutes and 53 seconds
% By Willi Mutschler, September 26, 2016. Email: willi@mutschler.eu
% Quadruplication Matrix as defined by
% Meijer (2005) - Matrix algebra for higher order moments. Linear Algebra and its Applications, 410,pp. 112134
function [DP6,DP6inv] = Q6_plication(p)
% Computes the 6-way duplication Matrix DP6 (and its Moore-Penrose inverse)
% such that for any p-dimensional vector x:
% y=kron(kron(kron(kron(kron(x,x),x,x),x),x)=DP6*z
% where z is of dimension np=p*(p+1)*(p+2)*(p+3)*(p+4)*(p+5)/(1*2*3*4*5*6)
% and is obtained from y by removing each second and later occurence of the
% same element. This is a generalization of the Duplication matrix.
% Reference: Meijer (2005) - Matrix algebra for higher order moments.
% Linear Algebra and its Applications, 410,pp. 112-134
% =========================================================================
% INPUTS
% * p [integer] size of vector
% -------------------------------------------------------------------------
% OUTPUTS
% * DP6 [p^6 by np] 6-way duplication matrix
% * DP6inv [np by np] Moore-Penrose inverse of DP6
% -------------------------------------------------------------------------
% This function is called by
% * pruned_state_space_system.m
% -------------------------------------------------------------------------
% This function calls
% * binom_coef (embedded)
% * mue (embedded)
% * uperm
% =========================================================================
% Copyright (C) 2020 Dynare Team
%
% Inputs:
% p: size of vector
% Outputs:
% QP: quadruplication matrix
% QPinv: Moore-Penrose inverse of QP
% This file is part of Dynare.
%
function [DP6,DP6inv] = Q6_plication(p,progress)
if nargin <2
progress =0;
end
reverseStr = ''; counti=1;
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
% =========================================================================
np = p*(p+1)*(p+2)*(p+3)*(p+4)*(p+5)/(1*2*3*4*5*6);
DP6 = spalloc(p^6,p*(p+1)*(p+2)*(p+3)*(p+4)*(p+5)/(1*2*3*4*5*6),p^6);
DP6 = spalloc(p^6,np,p^6);
counti=1;
for i1=1:p
for i2=i1:p
for i3=i2:p
for i4=i3:p
for i5=i4:p
for i6=i5:p
if progress && (rem(counti,100)== 0)
msg = sprintf(' Q6-plication Matrix Processed %d/%d', counti, np); fprintf([reverseStr, msg]); reverseStr = repmat(sprintf('\b'), 1, length(msg));
elseif progress && (counti==np)
msg = sprintf(' Q6-plication Matrix Processed %d/%d\n', counti, np); fprintf([reverseStr, msg]); reverseStr = repmat(sprintf('\b'), 1, length(msg));
end
idx = uperm([i6 i5 i4 i3 i2 i1]);
for r = 1:size(idx,1)
ii1 = idx(r,1); ii2= idx(r,2); ii3=idx(r,3); ii4=idx(r,4); ii5=idx(r,5); ii6=idx(r,6);
......@@ -45,11 +66,13 @@ end
DP6inv = (transpose(DP6)*DP6)\transpose(DP6);
function m = mue(p,i1,i2,i3,i4,i5,i6)
% Auxiliary expression, see page 122 of Meijer (2005)
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);
m = round(m);
end
function N = binom_coef(p,q,i)
% Auxiliary expression for binomial coefficients, see page 119 of Meijer (2005)
t = q; r =p+q-i;
if t==0
N=1;
......
%
% bivmom.m Date: 1/11/2004
% This Matlab program computes the product moment of X_1^{p_1}X_2^{p_2},
% where X_1 and X_2 are standard bivariate normally distributed.
function [y,dy] = bivmom(p,rho)
% Computes the product moment (and its derivative with respect to standard
% errors and correlation parameters) of X_1^{p_1}X_2^{p_2}, where X_1 and X_2
% are standard bivariate normally distributed.
% n : dimension of X
% rho: correlation coefficient between X_1 and X_2
% Reference: Kotz, Balakrishnan, and Johnson (2000), Continuous Multivariate
% Distributions, Vol. 1, p.261
% Note that there is a typo in Eq.(46.25), there should be an extra rho in front
% of the equation.
% Usage: bivmom(p,rho)
%
% Retrieved from http://www-2.rotman.utoronto.ca/~kan/papers/prodmom.zip
% This function is part of replication codes of the following paper:
% =========================================================================
% INPUTS
% p [2 by 1] powers of X_{1} and X_{2}
% rho [1 by 1] correlation coefficient between X_1 and X_2
% -------------------------------------------------------------------------
% OUTPUTS
% y [1 by 1] product moment E[X_1^{p_1}X_2^{p_2}]
% dy [1 by 1] derivative of y wrt to rho
% -------------------------------------------------------------------------
% This function is based upon bivmom.m which is part of replication codes
% of the following paper:
% Kan, R.: "From moments of sum to moments of product." Journal of
% Multivariate Analysis, 2008, vol. 99, issue 3, pages 542-554.
% bivmom.m can be retrieved from http://www-2.rotman.utoronto.ca/~kan/papers/prodmom.zip
% Further references:
% Kotz, Balakrishnan, and Johnson (2000), Continuous Multivariate Distributions, Vol. 1, p.261
% Note that there is a typo in Eq.(46.25), there should be an extra rho in front
% of the equation.
% =========================================================================
% Copyright (C) 2008-2015 Raymond Kan <kan@chass.utoronto.ca>
% Copyright (C) 2019-2020 Dynare Team
......@@ -33,7 +41,6 @@
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
% =========================================================================
function [y,dy] = bivmom(p,rho)
s1 = p(1);
s2 = p(2);
rho2 = rho^2;
......
......@@ -15,12 +15,13 @@ function k = commutation(n, m, sparseflag)
% This function is called by
% * get_perturbation_params_derivs.m (previously getH.m)
% * get_identification_jacobians.m (previously getJJ.m)
% * pruned_state_space_system.m
% -------------------------------------------------------------------------
% This function calls
% * vec (embedded)
% =========================================================================
% Copyright (C) 1997 Tom Minka <minka@microsoft.com>
% Copyright (C) 2019 Dynare Team
% Copyright (C) 2019-2020 Dynare Team
%
% This file is part of Dynare.
%
......@@ -47,24 +48,12 @@ if nargin < 3
sparseflag = 0;
end
if 0
% first method
i = 1:(n*m);
a = reshape(i, n, m);
j = vec(transpose(a));
k = zeros(n*m,n*m);
for r = i
k(r, j(r)) = 1;
end
if sparseflag
k = reshape(kron(vec(speye(n)), speye(m)), n*m, n*m);
else
% second method
k = reshape(kron(vec(eye(n)), eye(m)), n*m, n*m);
end
if sparseflag ~= 0
k = sparse(k);
end
function V = vec(A)
V = A(:);
end
......
......@@ -93,9 +93,9 @@ for jide = 1:4
elseif jide == 2
strTest = 'MINIMAL SYSTEM (Komunjer and Ng, 2011)'; strJacobian = 'Deltabar'; strMeaning = 'Jacobian of steady state and minimal system';
if options_ident.order == 2
strMeaning = 'Jacobian of second-order accurate mean and first-order minimal system';
strMeaning = 'Jacobian of first-order minimal system and second-order accurate mean';
elseif options_ident.order == 3
strMeaning = 'Jacobian of second-order accurate mean and first-order minimal system';
strMeaning = 'Jacobian of first-order minimal system and third-order accurate mean';
end
if ~no_identification_minimal
noidentification = 0; ide = ide_minimal;
......
This diff is collapsed.
......@@ -145,6 +145,12 @@ analytic_derivation_mode = options.analytic_derivation_mode;
% * -2: numerical two-sided finite difference method to compute numerically dYss, dg1, dg2, dg3, d2Yss and d2g1, the other output arguments are computed analytically as in kronflag=0
gstep = options.gstep;
order = options.order;
if isempty(options.qz_criterium)
% set default value for qz_criterium: if there are no unit roots one can use 1.0
% If they are possible, you may have have multiple unit roots and the accuracy
% decreases when computing the eigenvalues in lyapunov_symm. Hence, we normally use 1+1e-6
options = select_qz_criterium_value(options);
end
qz_criterium = options.qz_criterium;
threads_BC = options.threads.kronecker.sparse_hessian_times_B_kronecker_C;
......@@ -269,6 +275,11 @@ if d2flag
KalmanB = ghu;
end
% Store some objects
DERIVS.indpmodel = indpmodel;
DERIVS.indpstderr = indpstderr;
DERIVS.indpcorr = indpcorr;
if analytic_derivation_mode == -1
%% numerical two-sided finite difference method using function get_perturbation_params_derivs_numerical_objective.m (previously thet2tau.m in Dynare 4.5) for
% Jacobian (wrt selected stderr, corr and model parameters) of
......
......@@ -137,7 +137,7 @@ no_identification_spectrum = options_ident.no_identification_spectrum;
if info(1) == 0 %no errors in solution
% Compute parameter Jacobians for identification analysis
[MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dMINIMAL, derivatives_info] = get_identification_jacobians(estim_params_, M_, oo_, options_, options_ident, indpmodel, indpstderr, indpcorr, indvobs);
[MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dSPECTRUM_NO_MEAN, dMINIMAL, derivatives_info] = get_identification_jacobians(estim_params_, M_, oo_, options_, options_ident, indpmodel, indpstderr, indpcorr, indvobs);
if isempty(dMINIMAL)
% Komunjer and Ng is not computed if (1) minimality conditions are not fullfilled or (2) there are more shocks and measurement errors than observables, so we need to reset options
no_identification_minimal = 1;
......@@ -194,7 +194,7 @@ if info(1) == 0 %no errors in solution
options_ident.no_identification_spectrum = 1; %do not recompute dSPECTRUM
options_ident.ar = nlags; %store new lag number
options_.ar = nlags; %store new lag number
[~, ~, ~, ~, ~, ~, MOMENTS, dMOMENTS, ~, ~, ~] = get_identification_jacobians(estim_params_, M_, oo_, options_, options_ident, indpmodel, indpstderr, indpcorr, indvobs);
[~, ~, ~, ~, ~, ~, MOMENTS, dMOMENTS, ~, ~, ~, ~] = get_identification_jacobians(estim_params_, M_, oo_, options_, options_ident, indpmodel, indpstderr, indpcorr, indvobs);
ind_dMOMENTS = (find(max(abs(dMOMENTS'),[],1) > tol_deriv)); %new index with non-zero rows
end
......@@ -508,6 +508,7 @@ if info(1) == 0 %no errors in solution
ide_spectrum.norm_dSPECTRUM = norm_dSPECTRUM;
ide_spectrum.tilda_dSPECTRUM = tilda_dSPECTRUM;
ide_spectrum.dSPECTRUM = dSPECTRUM;
ide_spectrum.dSPECTRUM_NO_MEAN = dSPECTRUM_NO_MEAN;
end
%% Perform identification checks, i.e. find out which parameters are involved
......
......@@ -30,17 +30,15 @@ function out = identification_numerical_objective(params, outputflag, estim_para
% Jacobian of the dynamic model equations, and Y_t selected variables
% -------------------------------------------------------------------------
% This function is called by
% * get_first_order_solution_params_deriv.m (previously getH.m)
% * get_identification_jacobians.m (previously getJJ.m)
% -------------------------------------------------------------------------
% This function calls
% * [M.fname,'.dynamic']
% * dynare_resolve
% * dyn_vech
% * lyapunov_symm
% * resol
% * vec
% =========================================================================
% Copyright (C) 2011-2019 Dynare Team
% Copyright (C) 2011-2020 Dynare Team
%
% This file is part of Dynare.
%
......@@ -58,16 +56,6 @@ function out = identification_numerical_objective(params, outputflag, estim_para
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
% =========================================================================
if nargin < 11 || isempty(useautocorr)
useautocorr = 0;
end
if nargin < 12 || isempty(nlags)
nlags = 3;
end
if nargin < 13 || isempty(grid_nbr)
grid_nbr = 0;
end
%% Update stderr, corr and model parameters
%note that if no estimated_params_block is given, then all stderr and model parameters are selected but no corr parameters
if length(params) > length(indpmodel)
......@@ -90,50 +78,21 @@ end
%% compute Kalman transition matrices and steady state with updated parameters
[~,info,M,options,oo] = resol(0,M,options,oo);
options = rmfield(options,'options_ident');
oo.dr = pruned_state_space_system(M, options, oo.dr);
A = oo.dr.pruned.A;
B = oo.dr.pruned.B;
C = oo.dr.pruned.C(indvar,:);
D = oo.dr.pruned.D(indvar,:);
Om_z = oo.dr.pruned.Om_z;
Om_y = oo.dr.pruned.Om_y(indvar,indvar);
Varinov = oo.dr.pruned.Varinov;
obs_nbr = size(C,1);
pruned = pruned_state_space_system(M, options, oo.dr, indvar, nlags, useautocorr, 0);
%% out = [vech(cov(Y_t,Y_t)); vec(cov(Y_t,Y_{t-1}); ...; vec(cov(Y_t,Y_{t-nlags})] of indvar variables, in DR order. This is Iskrev (2010)'s J matrix.
if outputflag == 1
% Denote Ezz0 = E_t(z_t * z_t'), then the following Lyapunov equation defines the autocovariagrom: Ezz0 -A*Ezz*A' = B*Sig_e*B'
[Ezz0,u] = lyapunov_symm(A, Om_z, options.lyapunov_fixed_point_tol, options.qz_criterium, options.lyapunov_complex_threshold, 1, options.debug);
stationary_vars = (1:size(C,1))';
if ~isempty(u)
x = abs(C*u);
stationary_vars = find(all(x < options.Schur_vec_tol,2));
end
Eyy0 = NaN*ones(obs_nbr,obs_nbr);
Eyy0(stationary_vars,stationary_vars) = C(stationary_vars,:)*Ezz0*C(stationary_vars,:)' + Om_y(stationary_vars,stationary_vars);
indzeros = find(abs(Eyy0) < 1e-12); %find values that are numerical zero
Eyy0(indzeros) = 0;
if outputflag == 1
if useautocorr
sy = sqrt(diag(Ezz0)); %theoretical standard deviation
sy = sy(stationary_vars);
sy = sy*sy'; %cross products of standard deviations
sy0 = sy-diag(diag(sy))+eye(length(sy));
Eyy0corr = NaN*ones(size(C,1),size(C,1));
Eyy0corr(stationary_vars,stationary_vars) = Eyy0./sy0;
out = dyn_vech(Eyy0corr); %focus only on unique terms
out = dyn_vech(pruned.Corr_y);
else
out = dyn_vech(Eyy0); %focus only on unique terms
end
% compute autocovariances/autocorrelations of lagged observed variables
tmpEyyi = A*Ezz0*C(stationary_vars,:)' + B*Varinov*D(stationary_vars,:)';
Ai = eye(size(A,1)); %this is A^0
out = dyn_vech(pruned.Var_y);
end
for i = 1:nlags
Eyyi = NaN*ones(obs_nbr,obs_nbr);
Eyyi(stationary_vars,stationary_vars) = C(stationary_vars,:)*Ai*tmpEyyi;
if useautocorr
Eyyi = Eyyi./sy;
out = [out;vec(pruned.Corr_yi(:,:,i))];
else
out = [out;vec(pruned.Var_yi(:,:,i))];
end
out = [out;vec(Eyyi)];
Ai = Ai*A; %note that this is A^(i-1)
end
end
......@@ -142,13 +101,13 @@ if outputflag == 2
% This computes the spectral density g_omega where the interval [-pi;\pi] is discretized by grid_nbr points
freqs = (0 : pi/(grid_nbr/2):pi);% we focus only on positive values including the 0 frequency
tpos = exp( sqrt(-1)*freqs); %Fourier frequencies
IA = eye(size(A,1));
var_nbr = size(C,1);
IA = eye(size(pruned.A,1));
var_nbr = size(pruned.C,1);
out = zeros(var_nbr^2*length(freqs),1);
kk = 0;
for ig = 1:length(freqs)
Transferfct = D + C*((tpos(ig)*IA-A)\B);
g_omega = (1/(2*pi))*(Transferfct*Varinov*Transferfct'); % note that ' is the conjugate transpose
Transferfct = pruned.D + pruned.C*((tpos(ig)*IA-pruned.A)\pruned.B);
g_omega = (1/(2*pi))*(Transferfct*pruned.Varinov*Transferfct'); % note that ' is the conjugate transpose
kk = kk+1;
out(1 + (kk-1)*var_nbr^2 : kk*var_nbr^2) = g_omega(:);
end
......
function dy = prodmom_deriv(V,ii,nu,dV,dC)
% This function builds upon and extends prodmom.m to compute the
% derivatives of product moments of normally distributed variables with
% respect to standard errors and correlation parameters.
% prodmom.m is part of replication codes of the following paper:
function [y,dy] = prodmom_deriv(V,ii,nu,dV,dC)
% Computes the product moments (and its derivatives with respect to standard
% errors and correlation parameters) of normally distributed variables, i.e.
% this function computes the product moment of
% X_{i_1}^{nu_1}X_{i_2}^{nu_2}...X_{i_m}^{nu_m}, where X_{i_j} are elements
% from X ~ N(0_n,V) and V is positive semidefinite.
% Example: To get E[X_2X_4^3X_7^2], use prodmom_deriv(V,[2 4 7],[1 3 2])
% =========================================================================
% INPUTS
% V [n by n] covariance matrix of X (needs to be positive semidefinite)
% ii [m by 1] vector of i_j
% nu [nu_m by 1] power of X_{i_j}
% dV [n by n by stderrparam_nbr+corrparam_nbr] derivative of V with respect
% to selected standard error (stderr)
% and correlation (corr) parameters
% dC [n by n by stderrparam_nbr+corrparam_nbr] derivative of Correlation matrix C with respect
% to selected standard error (stderr)
% and correlation (corr) parameters
% -------------------------------------------------------------------------
% OUTPUTS
% y [1 by 1] product moment E[X_{i_1}^{nu_1}X_{i_2}^{nu_2}...X_{i_m}^{nu_m}]
% dy [1 by stderrparam_nbr+corrparam_nbr] derivatives of y wrt to selected
% standard error and corr parameters
% -------------------------------------------------------------------------
% This function is based upon prodmom.m which is part of replication codes
% of the following paper:
% Kan, R.: "From moments of sum to moments of product." Journal of
% Multivariate Analysis, 2008, vol. 99, issue 3, pages 542-554.
% prodmom.m can be retrieved from http://www-2.rotman.utoronto.ca/~kan/papers/prodmom.zip
% Further references:
% Triantafyllopoulos (2003) On the Central Moments of the Multidimensional
% Gaussian Distribution, Mathematical Scientist
% Kotz, Balakrishnan, and Johnson (2000), Continuous Multivariate
% Distributions, Vol. 1, p.261
% =========================================================================
% Copyright (C) 2008-2015 Raymond Kan <kan@chass.utoronto.ca>
% Copyright (C) 2019-2020 Dynare Team
......@@ -25,89 +52,118 @@ function dy = prodmom_deriv(V,ii,nu,dV,dC)
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
% =========================================================================
if nargin<3
nu = ones(size(ii));
nu = ones(size(ii));
end
s = sum(nu);
if s==0
dy = zeros(1,1,size(dV,3));
return
y = 1;
if nargout > 1
dy = zeros(1,1,size(dV,3));
end
return
end
if rem(s,2)==1
dy = zeros(1,1,size(dV,3));
return
y = 0;
if nargout > 1
dy = zeros(1,1,size(dV,3));
end
return
end
nuz = nu==0;
nu(nuz) = [];
ii(nuz) = [];
m = length(ii);
V = V(ii,ii);
dV = dV(ii,ii,:);
if nargout > 1
dV = dV(ii,ii,:);
end
s2 = s/2;
%
% Use univariate normal results
%
if m==1
dy = s2*V^(s2-1)*dV*prod([1:2:s-1]);
dy = reshape(dy,1,size(dV,3));
y = V^s2*prod([1:2:s-1]);
if nargout > 1
dy = s2*V^(s2-1)*dV*prod([1:2:s-1]);
dy = reshape(dy,1,size(dV,3));
end
return
end
%
% Use bivariate normal results when there are only two distinct indices
%
if m==2
rho = V(1,2)/sqrt(V(1,1)*V(2,2));
drho = dC(ii(1),ii(2),:);
[tmp,dtmp] = bivmom(nu,rho);
dy = (nu(1)/2)*V(1,1)^(nu(1)/2-1)*dV(1,1,:) * V(2,2)^(nu(2)/2) * tmp...
+ V(1,1)^(nu(1)/2) * (nu(2)/2)*V(2,2)^(nu(2)/2-1)*dV(2,2,:) * tmp...
+ V(1,1)^(nu(1)/2) * V(2,2)^(nu(2)/2) * dtmp * drho;
dy = reshape(dy,1,size(dV,3));
return
rho = V(1,2)/sqrt(V(1,1)*V(2,2));
if nargout > 1
drho = dC(ii(1),ii(2),:);
[tmp,dtmp] = bivmom(nu,rho);
dy = (nu(1)/2)*V(1,1)^(nu(1)/2-1)*dV(1,1,:) * V(2,2)^(nu(2)/2) * tmp...
+ V(1,1)^(nu(1)/2) * (nu(2)/2)*V(2,2)^(nu(2)/2-1)*dV(2,2,:) * tmp...
+ V(1,1)^(nu(1)/2) * V(2,2)^(nu(2)/2) * dtmp * drho;
dy = reshape(dy,1,size(dV,3));
else
tmp = bivmom(nu,rho);
end
y = V(1,1)^(nu(1)/2)*V(2,2)^(nu(2)/2)*tmp;
return
end
%
% Regular case
%
[nu,inu] = sort(nu,2,'descend');
V = V(inu,inu); % Extract only the relevant part of V
dV = dV(inu,inu,:); % Extract only the relevant part of dV
x = zeros(1,m);
V = V./2;
dV = dV./2;
nu2 = nu./2;
p = 2;
q = nu2*V*nu2';
%dq = nu2*dV*nu2';
%dq = multiprod(multiprod(nu2,dV),nu2');
dq = NaN(size(q,1), size(q,2), size(dV,3));
for jp = 1:size(dV,3)
dq(:,:,jp) = nu2*dV(:,:,jp)*nu2';
y = 0;
if nargout > 1
dV = dV(inu,inu,:); % Extract only the relevant part of dV
dV = dV./2;
%dq = nu2*dV*nu2';
%dq = multiprod(multiprod(nu2,dV),nu2');
dq = NaN(size(q,1), size(q,2), size(dV,3));
for jp = 1:size(dV,3)
dq(:,:,jp) = nu2*dV(:,:,jp)*nu2';
end
dy = 0;
end
dy = 0;
for i=1:fix(prod(nu+1)/2)
dy = dy+p*s2*q^(s2-1)*dq;
y = y+p*q^s2;
if nargout > 1
dy = dy+p*s2*q^(s2-1)*dq;
end
for j=1:m
if x(j)<nu(j)
x(j) = x(j)+1;
p = -round(p*(nu(j)+1-x(j))/x(j));
%dq = dq-2*(nu2-x)*dV(:,j,:)-dV(j,j,:);
%dq = dq-2*multiprod((nu2-x),dV(:,j,:))-dV(j,j,:);
for jp=1:size(dV,3)
dq(:,:,jp) = dq(:,:,jp)-2*(nu2-x)*dV(:,j,jp)-dV(j,j,jp);
end
q = q-2*(nu2-x)*V(:,j)-V(j,j);
break
x(j) = x(j)+1;
p = -round(p*(nu(j)+1-x(j))/x(j));
q = q-2*(nu2-x)*V(:,j)-V(j,j);
if nargout > 1
%dq = dq-2*(nu2-x)*dV(:,j,:)-dV(j,j,:);
%dq = dq-2*multiprod((nu2-x),dV(:,j,:))-dV(j,j,:);
for jp=1:size(dV,3)
dq(:,:,jp) = dq(:,:,jp)-2*(nu2-x)*dV(:,j,jp)-dV(j,j,jp);
end
end
break
else
x(j) = 0;
if rem(nu(j),2)==1
p = -p;
end
%dq = dq+2*nu(j)*multiprod((nu2-x),dV(:,j,:))-nu(j)^2*dV(j,j,:);
for jp=1:size(dV,3)
dq(:,:,jp) = dq(:,:,jp)+2*nu(j)*(nu2-x)*dV(:,j,jp)-nu(j)^2*dV(j,j,jp);
end
q = q+2*nu(j)*(nu2-x)*V(:,j)-nu(j)^2*V(j,j);
x(j) = 0;
if rem(nu(j),2)==1
p = -p;
end
if nargout > 1
%dq = dq+2*nu(j)*multiprod((nu2-x),dV(:,j,:))-nu(j)^2*dV(j,j,:);
for jp=1:size(dV,3)
dq(:,:,jp) = dq(:,:,jp)+2*nu(j)*(nu2-x)*dV(:,j,jp)-nu(j)^2*dV(j,j,jp);
end
end
q = q+2*nu(j)*(nu2-x)*V(:,j)-nu(j)^2*V(j,j);
end
end
end
dy = dy/prod([1:s2]);
dy = reshape(dy,1,size(dV,3));
\ No newline at end of file
y = y/prod([1:s2]);
if nargout > 1
dy = dy/prod([1:s2]);
dy = reshape(dy,1,size(dV,3));
end
\ No newline at end of file
This source diff could not be displayed because it is too large. You can view the blob instead.
% By Willi Mutschler, September 26, 2016. Email: willi@mutschler.eu
% Quadruplication Matrix as defined by
% Meijer (2005) - Matrix algebra for higher order moments. Linear Algebra and its Applications, 410,pp. 112134
function [QP,QPinv] = quadruplication(p)
% Computes the Quadruplication Matrix QP (and its Moore-Penrose inverse)
% such that for any p-dimensional vector x:
% y=kron(kron(kron(x,x),x),x)=QP*z
% where z is of dimension np=p*(p+1)*(p+2)*(p+3)/2 and is obtained from y
% by removing each second and later occurence of the same element.
% This is a generalization of the Duplication matrix.
% Reference: Meijer (2005) - Matrix algebra for higher order moments.
% Linear Algebra and its Applications, 410,pp. 112-134
% =========================================================================
% INPUTS
% * p [integer] size of vector
% -------------------------------------------------------------------------
% OUTPUTS
% * QP [p^4 by np] Quadruplication matrix
% * QPinv [np by np] Moore-Penrose inverse of QP
% -------------------------------------------------------------------------
% This function is called by
% * pruned_state_space_system.m
% -------------------------------------------------------------------------
% This function calls
% * mue (embedded)
% * uperm
% =========================================================================
% Copyright (C) 2020 Dynare Team
%
% Inputs:
% p: size of vector
% Outputs:
% QP: quadruplication matrix
% QPinv: Moore-Penrose inverse of QP
% This file is part of Dynare.
%
function [QP,QPinv] = quadruplication(p,progress,sparseflag)
if nargin <2
progress =0;
end
if nargin < 3
sparseflag = 1;
end
reverseStr = ''; counti=1;
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.