diff --git a/matlab/Q6_plication.m b/matlab/Q6_plication.m index c08e861f3..0692364f1 100644 --- a/matlab/Q6_plication.m +++ b/matlab/Q6_plication.m @@ -1,33 +1,54 @@ -% 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. 112–134 +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 . +% ========================================================================= 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; diff --git a/matlab/bivmom.m b/matlab/bivmom.m index f759a7f9d..f850fe1c2 100644 --- a/matlab/bivmom.m +++ b/matlab/bivmom.m @@ -1,19 +1,27 @@ -% -% 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 % 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 . % ========================================================================= -function [y,dy] = bivmom(p,rho) s1 = p(1); s2 = p(2); rho2 = rho^2; diff --git a/matlab/commutation.m b/matlab/commutation.m index c03bf9a0e..a87017e9b 100644 --- a/matlab/commutation.m +++ b/matlab/commutation.m @@ -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 -% 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 diff --git a/matlab/disp_identification.m b/matlab/disp_identification.m index 716d52f63..29142c2fa 100644 --- a/matlab/disp_identification.m +++ b/matlab/disp_identification.m @@ -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; diff --git a/matlab/get_identification_jacobians.m b/matlab/get_identification_jacobians.m index 0620b7f48..7e7443bfd 100644 --- a/matlab/get_identification_jacobians.m +++ b/matlab/get_identification_jacobians.m @@ -1,4 +1,4 @@ -function [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) +function [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) % function [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) % previously getJJ.m in Dynare 4.5 % Sets up the Jacobians needed for identification analysis @@ -54,7 +54,14 @@ function [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dM % * order==1: corresponds to Iskrev (2010)'s J matrix % * order==2: corresponds to Mutschler (2015)'s \bar{M}_2 matrix, i.e. theoretical moments from the pruned state space system % -% dSPECTRUM: [totparam_nbr by totparam_nbr] in DR order. Gram matrix of Jacobian (wrt all params) of spectral density for VAROBS variables, where +% dSPECTRUM: [totparam_nbr by totparam_nbr] in DR order. Gram matrix of Jacobian (wrt all params) of mean and of spectral density for VAROBS variables, where +% spectral density at frequency w: f(w) = (2*pi)^(-1)*H(exp(-i*w))*E[Inov*Inov']*ctranspose(H(exp(-i*w)) with H being the Transfer function +% dSPECTRUM = dMEAN*dMEAN + int_{-\pi}^\pi transpose(df(w)/dp')*(df(w)/dp') dw +% * order==1: corresponds to Qu and Tkachenko (2012)'s G matrix, where Inov and H are computed from linear state space system +% * order==2: corresponds to Mutschler (2015)'s G_2 matrix, where Inov and H are computed from second-order pruned state space system +% * order==3: Inov and H are computed from third-order pruned state space system +% +% dSPECTRUM_NO_MEAN:[totparam_nbr by totparam_nbr] in DR order. Gram matrix of Jacobian (wrt all params) of spectral density for VAROBS variables, where % spectral density at frequency w: f(w) = (2*pi)^(-1)*H(exp(-i*w))*E[Inov*Inov']*ctranspose(H(exp(-i*w)) with H being the Transfer function % dSPECTRUM = int_{-\pi}^\pi transpose(df(w)/dp')*(df(w)/dp') dw % * order==1: corresponds to Qu and Tkachenko (2012)'s G matrix, where Inov and H are computed from linear state space system @@ -81,11 +88,11 @@ function [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dM % * get_minimal_state_representation % * duplication % * dyn_vech +% * fjaco % * get_perturbation_params_derivs (previously getH) % * get_all_parameters -% * fjaco -% * lyapunov_symm % * identification_numerical_objective (previously thet2tau) +% * pruned_state_space_system % * vec % ========================================================================= % Copyright (C) 2010-2020 Dynare Team @@ -220,22 +227,9 @@ elseif order == 3 end % Get (pruned) state space representation: -options.options_ident.indvobs = indvobs; -options.options_ident.indpmodel = indpmodel; -options.options_ident.indpstderr = indpstderr; -options.options_ident.indpcorr = indpcorr; -oo.dr = pruned_state_space_system(M, options, oo.dr); -MEAN = oo.dr.pruned.E_y; dMEAN = oo.dr.pruned.dE_y; -A = oo.dr.pruned.A; dA = oo.dr.pruned.dA; -B = oo.dr.pruned.B; dB = oo.dr.pruned.dB; -C = oo.dr.pruned.C; dC = oo.dr.pruned.dC; -D = oo.dr.pruned.D; dD = oo.dr.pruned.dD; -c = oo.dr.pruned.c; dc = oo.dr.pruned.dc; -d = oo.dr.pruned.d; dd = oo.dr.pruned.dd; -Varinov = oo.dr.pruned.Varinov; dVarinov = oo.dr.pruned.dVarinov; -Om_z = oo.dr.pruned.Om_z; dOm_z = oo.dr.pruned.dOm_z; -Om_y = oo.dr.pruned.Om_y; dOm_y = oo.dr.pruned.dOm_y; - +pruned = pruned_state_space_system(M, options, oo.dr, indvobs, nlags, useautocorr, 1); +MEAN = pruned.E_y; +dMEAN = pruned.dE_y; %storage for Jacobians used in dsge_likelihood.m for analytical Gradient and Hession of likelihood (only at order=1) derivatives_info = struct(); if order == 1 @@ -247,12 +241,19 @@ if order == 1 end %% Compute dMOMENTS -%Note that our state space system for computing second moments is the following: -% zhat = A*zhat(-1) + B*xi, where zhat = z - E(z) -% yhat = C*zhat(-1) + D*xi, where yhat = y - E(y) -if ~no_identification_moments - MOMENTS = identification_numerical_objective(xparam1, 1, estim_params, M, oo, options, indpmodel, indpstderr, indpcorr, indvobs, useautocorr, nlags, grid_nbr); %[outputflag=1] - MOMENTS = [MEAN; MOMENTS]; +if ~no_identification_moments + if useautocorr + E_yy = pruned.Corr_y; dE_yy = pruned.dCorr_y; + E_yyi = pruned.Corr_yi; dE_yyi = pruned.dCorr_yi; + else + E_yy = pruned.Var_y; dE_yy = pruned.dVar_y; + E_yyi = pruned.Var_yi; dE_yyi = pruned.dVar_yi; + end + MOMENTS = [MEAN; dyn_vech(E_yy)]; + for i=1:nlags + MOMENTS = [MOMENTS; vec(E_yyi(:,:,i))]; + end + if kronflag == -1 %numerical derivative of autocovariogram dMOMENTS = fjaco(str2func('identification_numerical_objective'), xparam1, 1, estim_params, M, oo, options, indpmodel, indpstderr, indpcorr, indvobs, useautocorr, nlags, grid_nbr); %[outputflag=1] @@ -260,58 +261,10 @@ if ~no_identification_moments else dMOMENTS = zeros(obs_nbr + obs_nbr*(obs_nbr+1)/2 + nlags*obs_nbr^2 , totparam_nbr); dMOMENTS(1:obs_nbr,:) = dMEAN; %add Jacobian of first moments of VAROBS variables - % Denote Ezz0 = E[zhat*zhat'], then the following Lyapunov equation defines the autocovariagram: Ezz0 -A*Ezz0*A' = B*Sig_xi*B' = Om_z - [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:length(indvobs))'; - 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); - %here method=1 is used, whereas all other calls of lyapunov_symm use method=2. The reason is that T and U are persistent, and input matrix A is the same, so using option 2 for all the rest of iterations spares a lot of computing time while not repeating Schur every time - indzeros = find(abs(Eyy0) < 1e-12); %find values that are numerical zero - Eyy0(indzeros) = 0; - if useautocorr - sdy = sqrt(diag(Eyy0)); %theoretical standard deviation - sdy = sdy(stationary_vars); - sy = sdy*sdy'; %cross products of standard deviations - end - for jp = 1:totparam_nbr - if jp <= (stderrparam_nbr+corrparam_nbr) - %Note that for stderr and corr parameters, the derivatives of the system matrices are zero, i.e. dA=dB=dC=dD=0 - dEzz0 = lyapunov_symm(A,dOm_z(:,:,jp),options.lyapunov_fixed_point_tol,options.qz_criterium,options.lyapunov_complex_threshold,2,options.debug); %here method=2 is used to spare a lot of computing time while not repeating Schur every time - dEyy0 = C*dEzz0*C' + dOm_y(:,:,jp); - else %model parameters - dEzz0 = lyapunov_symm(A,dOm_z(:,:,jp)+dA(:,:,jp)*Ezz0*A'+A*Ezz0*dA(:,:,jp)',options.lyapunov_fixed_point_tol,options.qz_criterium,options.lyapunov_complex_threshold,2,options.debug); %here method=2 is used to spare a lot of computing time while not repeating Schur every time - dEyy0 = dC(:,:,jp)*Ezz0*C' + C*dEzz0*C' + C*Ezz0*dC(:,:,jp)' + dOm_y(:,:,jp); - end - %indzeros = find(abs(dEyy0) < 1e-12); - %dEyy0(indzeros) = 0; - if useautocorr - dsy = 1/2./sdy.*diag(dEyy0); - dsy = dsy*sdy'+sdy*dsy'; - dEyy0corr = (dEyy0.*sy-dsy.*Eyy0)./(sy.*sy); - dEyy0corr = dEyy0corr-diag(diag(dEyy0corr))+diag(diag(dEyy0)); - dMOMENTS(obs_nbr+1 : obs_nbr+obs_nbr*(obs_nbr+1)/2 , jp) = dyn_vech(dEyy0corr); %focus only on VAROBS variables - else - dMOMENTS(obs_nbr+1 : obs_nbr+obs_nbr*(obs_nbr+1)/2 , jp) = dyn_vech(dEyy0); %focus only on VAROBS variables - end - tmpEyyi = A*Ezz0*C(stationary_vars,:)' + B*Varinov*D(stationary_vars,:)'; - %we could distinguish between stderr and corr params, but this has no real speed effect as we multipliy with zeros - dtmpEyyi = dA(:,:,jp)*Ezz0*C' + A*dEzz0*C' + A*Ezz0*dC(:,:,jp)' + dB(:,:,jp)*Varinov*D' + B*dVarinov(:,:,jp)*D' + B*Varinov*dD(:,:,jp)'; - Ai = eye(size(A,1)); %this is A^0 - dAi = zeros(size(A,1),size(A,1)); %this is d(A^0) + dMOMENTS(obs_nbr+1 : obs_nbr+obs_nbr*(obs_nbr+1)/2 , jp) = dyn_vech(dE_yy(:,:,jp)); for i = 1:nlags - Eyyi = C(stationary_vars,:)*Ai*tmpEyyi; - dEyyi = dC(:,:,jp)*Ai*tmpEyyi + C*dAi*tmpEyyi + C*Ai*dtmpEyyi; - if useautocorr - dEyyi = (dEyyi.*sy-dsy.*Eyyi)./(sy.*sy); - end - dMOMENTS(obs_nbr + obs_nbr*(obs_nbr+1)/2 + (i-1)*obs_nbr^2 + 1 : obs_nbr + obs_nbr*(obs_nbr+1)/2 + i*obs_nbr^2, jp) = vec(dEyyi); %focus only on VAROBS variables - dAi = dAi*A + Ai*dA(:,:,jp); %note that this is d(A^(i-1)) - Ai = Ai*A; %note that this is A^(i-1) + dMOMENTS(obs_nbr + obs_nbr*(obs_nbr+1)/2 + (i-1)*obs_nbr^2 + 1 : obs_nbr + obs_nbr*(obs_nbr+1)/2 + i*obs_nbr^2, jp) = vec(dE_yyi(:,:,i,jp)); end end end @@ -357,7 +310,7 @@ if ~no_identification_spectrum freqs = (0 : pi/(grid_nbr/2):pi); % we focus only on positive frequencies tpos = exp( sqrt(-1)*freqs); %positive Fourier frequencies tneg = exp(-sqrt(-1)*freqs); %negative Fourier frequencies - IA = eye(size(A,1)); + IA = eye(size(pruned.A,1)); if kronflag == -1 %numerical derivative of spectral density dOmega_tmp = fjaco(str2func('identification_numerical_objective'), xparam1, 2, estim_params, M, oo, options, indpmodel, indpstderr, indpcorr, indvobs, useautocorr, nlags, grid_nbr); %[outputflag=2] @@ -366,68 +319,63 @@ if ~no_identification_spectrum kk = kk+1; dOmega = dOmega_tmp(1 + (kk-1)*obs_nbr^2 : kk*obs_nbr^2,:); if ig == 1 % add zero frequency once - dSPECTRUM = dOmega'*dOmega; + dSPECTRUM_NO_MEAN = dOmega'*dOmega; else % due to symmetry to negative frequencies we can add positive frequencies twice - dSPECTRUM = dSPECTRUM + 2*(dOmega'*dOmega); + dSPECTRUM_NO_MEAN = dSPECTRUM_NO_MEAN + 2*(dOmega'*dOmega); end end elseif kronflag == 1 %use Kronecker products - dA = reshape(dA,size(dA,1)*size(dA,2),size(dA,3)); - dB = reshape(dB,size(dB,1)*size(dB,2),size(dB,3)); - dC = reshape(dC,size(dC,1)*size(dC,2),size(dC,3)); - dD = reshape(dD,size(dD,1)*size(dD,2),size(dD,3)); - dVarinov = reshape(dVarinov,size(dVarinov,1)*size(dVarinov,2),size(dVarinov,3)); - K_obs_exo = commutation(obs_nbr,size(Varinov,1)); + dA = reshape(pruned.dA,size(pruned.dA,1)*size(pruned.dA,2),size(pruned.dA,3)); + dB = reshape(pruned.dB,size(pruned.dB,1)*size(pruned.dB,2),size(pruned.dB,3)); + dC = reshape(pruned.dC,size(pruned.dC,1)*size(pruned.dC,2),size(pruned.dC,3)); + dD = reshape(pruned.dD,size(pruned.dD,1)*size(pruned.dD,2),size(pruned.dD,3)); + dVarinov = reshape(pruned.dVarinov,size(pruned.dVarinov,1)*size(pruned.dVarinov,2),size(pruned.dVarinov,3)); + K_obs_exo = commutation(obs_nbr,size(pruned.Varinov,1)); for ig=1:length(freqs) z = tneg(ig); - zIminusA = (z*IA - A); + zIminusA = (z*IA - pruned.A); zIminusAinv = zIminusA\IA; - Transferfct = D + C*zIminusAinv*B; % Transfer function + Transferfct = pruned.D + pruned.C*zIminusAinv*pruned.B; % Transfer function dzIminusA = -dA; dzIminusAinv = kron(-(transpose(zIminusA)\IA),zIminusAinv)*dzIminusA; %this takes long - dTransferfct = dD + DerivABCD(C,dC,zIminusAinv,dzIminusAinv,B,dB); %also long + dTransferfct = dD + DerivABCD(pruned.C,dC,zIminusAinv,dzIminusAinv,pruned.B,dB); %this takes long dTransferfct_conjt = K_obs_exo*conj(dTransferfct); - dOmega = (1/(2*pi))*DerivABCD(Transferfct,dTransferfct,Varinov,dVarinov,Transferfct',dTransferfct_conjt); %also long + dOmega = (1/(2*pi))*DerivABCD(Transferfct,dTransferfct,pruned.Varinov,dVarinov,Transferfct',dTransferfct_conjt); %also long if ig == 1 % add zero frequency once - dSPECTRUM = dOmega'*dOmega; + dSPECTRUM_NO_MEAN = dOmega'*dOmega; else % due to symmetry to negative frequencies we can add positive frequencies twice - dSPECTRUM = dSPECTRUM + 2*(dOmega'*dOmega); + dSPECTRUM_NO_MEAN = dSPECTRUM_NO_MEAN + 2*(dOmega'*dOmega); end end - %put back into tensor notation - dA = reshape(dA,size(A,1),size(A,2),totparam_nbr); - dB = reshape(dB,size(B,1),size(B,2),totparam_nbr); - dC = reshape(dC,size(C,1),size(C,2),totparam_nbr); - dD = reshape(dD,size(D,1),size(D,2),totparam_nbr); - dVarinov = reshape(dVarinov,size(Varinov,1),size(Varinov,2),totparam_nbr); elseif (kronflag==0) || (kronflag==-2) for ig = 1:length(freqs) - IzminusA = tpos(ig)*IA - A; - invIzminusA = IzminusA\eye(size(A,1)); - Transferfct = D + C*invIzminusA*B; + IzminusA = tpos(ig)*IA - pruned.A; + invIzminusA = IzminusA\eye(size(pruned.A,1)); + Transferfct = pruned.D + pruned.C*invIzminusA*pruned.B; dOmega = zeros(obs_nbr^2,totparam_nbr); for j = 1:totparam_nbr if j <= stderrparam_nbr+corrparam_nbr %stderr and corr parameters: only dSig is nonzero - dOmega_tmp = Transferfct*dVarinov(:,:,j)*Transferfct'; + dOmega_tmp = Transferfct*pruned.dVarinov(:,:,j)*Transferfct'; else %model parameters - dinvIzminusA = -invIzminusA*(-dA(:,:,j))*invIzminusA; - dTransferfct = dD(:,:,j) + dC(:,:,j)*invIzminusA*B + C*dinvIzminusA*B + C*invIzminusA*dB(:,:,j); - dOmega_tmp = dTransferfct*Varinov*Transferfct' + Transferfct*dVarinov(:,:,j)*Transferfct' + Transferfct*Varinov*dTransferfct'; + dinvIzminusA = -invIzminusA*(-pruned.dA(:,:,j))*invIzminusA; + dTransferfct = pruned.dD(:,:,j) + pruned.dC(:,:,j)*invIzminusA*pruned.B + pruned.C*dinvIzminusA*pruned.B + pruned.C*invIzminusA*pruned.dB(:,:,j); + dOmega_tmp = dTransferfct*pruned.Varinov*Transferfct' + Transferfct*pruned.dVarinov(:,:,j)*Transferfct' + Transferfct*pruned.Varinov*dTransferfct'; end dOmega(:,j) = (1/(2*pi))*dOmega_tmp(:); end if ig == 1 % add zero frequency once - dSPECTRUM = dOmega'*dOmega; + dSPECTRUM_NO_MEAN = dOmega'*dOmega; else % due to symmetry to negative frequencies we can add positive frequencies twice - dSPECTRUM = dSPECTRUM + 2*(dOmega'*dOmega); + dSPECTRUM_NO_MEAN = dSPECTRUM_NO_MEAN + 2*(dOmega'*dOmega); end end end % Normalize Matrix and add steady state Jacobian, note that G is real and symmetric by construction - dSPECTRUM = 2*pi*dSPECTRUM./(2*length(freqs)-1) + dMEAN'*dMEAN; - dSPECTRUM = real(dSPECTRUM); + dSPECTRUM_NO_MEAN = real(2*pi*dSPECTRUM_NO_MEAN./(2*length(freqs)-1)); + dSPECTRUM = dSPECTRUM_NO_MEAN + dMEAN'*dMEAN; else + dSPECTRUM_NO_MEAN = []; dSPECTRUM = []; end @@ -442,14 +390,14 @@ if ~no_identification_minimal dMINIMAL = []; else % Derive and check minimal state vector of first-order - SYS.A = oo.dr.ghx(oo.dr.pruned.indx,:); - SYS.dA = oo.dr.derivs.dghx(oo.dr.pruned.indx,:,:); - SYS.B = oo.dr.ghu(oo.dr.pruned.indx,:); - SYS.dB = oo.dr.derivs.dghu(oo.dr.pruned.indx,:,:); - SYS.C = oo.dr.ghx(oo.dr.pruned.indy,:); - SYS.dC = oo.dr.derivs.dghx(oo.dr.pruned.indy,:,:); - SYS.D = oo.dr.ghu(oo.dr.pruned.indy,:); - SYS.dD = oo.dr.derivs.dghu(oo.dr.pruned.indy,:,:); + SYS.A = oo.dr.ghx(pruned.indx,:); + SYS.dA = oo.dr.derivs.dghx(pruned.indx,:,:); + SYS.B = oo.dr.ghu(pruned.indx,:); + SYS.dB = oo.dr.derivs.dghu(pruned.indx,:,:); + SYS.C = oo.dr.ghx(pruned.indy,:); + SYS.dC = oo.dr.derivs.dghx(pruned.indy,:,:); + SYS.D = oo.dr.ghu(pruned.indy,:); + SYS.dD = oo.dr.derivs.dghu(pruned.indy,:,:); [CheckCO,minnx,SYS] = get_minimal_state_representation(SYS,1); if CheckCO == 0 diff --git a/matlab/get_perturbation_params_derivs.m b/matlab/get_perturbation_params_derivs.m index 08827e9b5..2deb29da7 100644 --- a/matlab/get_perturbation_params_derivs.m +++ b/matlab/get_perturbation_params_derivs.m @@ -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 diff --git a/matlab/identification_analysis.m b/matlab/identification_analysis.m index 5c31d0135..6395a51b3 100644 --- a/matlab/identification_analysis.m +++ b/matlab/identification_analysis.m @@ -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 diff --git a/matlab/identification_numerical_objective.m b/matlab/identification_numerical_objective.m index a55d4d8ee..48ba28bfc 100644 --- a/matlab/identification_numerical_objective.m +++ b/matlab/identification_numerical_objective.m @@ -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 . % ========================================================================= -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 diff --git a/matlab/prodmom_deriv.m b/matlab/prodmom_deriv.m index b4cd75f66..c154df8a7 100644 --- a/matlab/prodmom_deriv.m +++ b/matlab/prodmom_deriv.m @@ -1,10 +1,37 @@ -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 % 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 . % ========================================================================= 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) 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 diff --git a/matlab/pruned_state_space_system.m b/matlab/pruned_state_space_system.m index b0c5c4f04..df5cf5ab4 100644 --- a/matlab/pruned_state_space_system.m +++ b/matlab/pruned_state_space_system.m @@ -1,763 +1,1219 @@ -function dr = pruned_state_space_system(M, options, dr) -% This can be set up into the following ABCD representation: -% z = c + A*z(-1) + B*xi -% y = ybar + d + C*z(-1) + D*xi +function pruned_state_space = pruned_state_space_system(M, options, dr, indy, nlags, useautocorr, compute_derivs) +% Set up the pruned state space ABCD representation: +% z = c + A*z(-1) + B*inov +% y = ys + d + C*z(-1) + D*inov +% References: +% - Andreasen, Martin M., Jesús Fernández-Villaverde and Juan F. Rubio-Ramírez (2018): +% "The Pruned State-Space System for Non-Linear DSGE Models: Theory and Empirical Applications", +% Review of Economic Studies, Volume 85, Issue 1, Pages 1–49. +% - Mutschler, Willi (2018): "Higher-order statistics for DSGE models", +% Econometrics and Statistics, Volume 6, Pages 44-56. +% ========================================================================= +% INPUTS +% M: [structure] storing the model information +% options: [structure] storing the options +% dr: [structure] storing the results from perturbation approximation +% indy: [vector] index of control variables in DR order +% nlags: [integer] number of lags in autocovariances and autocorrelations +% useautocorr: [boolean] true: compute autocorrelations +% ------------------------------------------------------------------------- +% OUTPUTS +% pruned_state_space: [structure] with the following fields: +% indx: [x_nbr by 1] +% index of state variables +% indy: [y_nbr by 1] +% index of control variables +% A: [z_nbr by z_nbr] +% state space transition matrix A mapping previous states to current states +% B: [z_nbr by inov_nbr] +% state space transition matrix B mapping current inovations to current states +% c: [z_nbr by 1] +% state space transition matrix c mapping constants to current states +% C: [y_nbr by z_nbr] +% state space measurement matrix C mapping previous states to current controls +% D: [y_nbr by inov_nbr] +% state space measurement matrix D mapping current inovations to current controls +% d: [y_nbr by 1] +% state space measurement matrix d mapping constants to current controls +% Var_inov [inov_nbr by inov_nbr] +% contemporenous covariance matrix of innovations, i.e. E[inov*inov'] +% Var_z [z_nbr by z_nbr] +% contemporenous covariance matrix of states z +% Var_y [y_nbr by y_nbr] +% contemporenous covariance matrix of controls y +% Var_yi [y_nbr by y_nbr by nlags] +% autocovariance matrix of controls y +% Corr_y [y_nbr by y_nbr] +% contemporenous correlation matrix of controls y +% Corr_yi [y_nbr by y_nbr by nlags] +% autocorrelation matrix of controls y +% E_y [y_nbr by 1] +% unconditional theoretical mean of control variables y +% +% if compute_derivs == 1, then the following additional fields are outputed: +% dA: [z_nbr by z_nbr by totparam_nbr] +% parameter Jacobian of A +% dB: [z_nbr by inov_nbr by totparam_nbr] +% parameter Jacobian of B +% dc: [z_nbr by totparam_nbr] +% parameter Jacobian of c +% dC: [y_nbr by z_nbr by totparam_nbr] +% parameter Jacobian of C +% dD: [y_nbr by inov_nbr by totparam_nbr] +% parameter Jacobian of D +% dd: [y_nbr by totparam_nbr] +% parameter Jacobian of d +% dVar_inov [inov_nbr by inov_nbr by totparam_nbr] +% parameter Jacobian of Var_inov +% dVar_z [z_nbr by z_nbr by totparam_nbr] +% parameter Jacobian of Var_z +% dVar_y [y_nbr by y_nbr by totparam_nbr] +% parameter Jacobian of Var_y +% dVar_yi [y_nbr by y_nbr by nlags by totparam_nbr] +% parameter Jacobian of Var_yi +% dCorr_y [y_nbr by y_nbr by totparam_nbr] +% parameter Jacobian of Corr_y +% dCorr_yi [y_nbr by y_nbr by nlags by totparam_nbr] +% parameter Jacobian of Corr_yi +% dE_y [y_nbr by totparam_nbr] +% parameter Jacobian of E_y +% ------------------------------------------------------------------------- +% This function is called by +% * get_identification_jacobians.m +% * identification_numerical_objective.m +% ------------------------------------------------------------------------- +% This function calls +% * allVL1.m +% * commutation.m +% * disclyap_fast.m +% * duplication.m +% * lyapunov_symm.m +% * prodmom +% * prodmom_deriv +% * Q6_plication +% * quadruplication.m +% * vec.m +% ========================================================================= +% Copyright (C) 2019-2020 Dynare Team +% +% This file is part of Dynare. +% +% 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 . +% ========================================================================= -% Denote -% hx = dr.ghx(indx,:); hu = dr.ghu(indx,:); -% hxx = dr.ghxx(indx,:); hxu = dr.ghxu(indx,:); huu = dr.ghuu(indx,:); hs2 = dr.ghs2(indx,:); -% gx = dr.ghx(indy,:); gu = dr.ghu(indy,:); -% gxx = dr.ghxx(indy,:); gxu = dr.ghxu(indy,:); guu = dr.ghuu(indy,:); gs2 = dr.ghs2(indy,:); -% hxss = dr.ghxss(indx,:); huss = dr.ghuss(indx,:); hxxx = dr.ghxxx(indx,:); huuu = dr.ghuuu(indx,:); hxxu = dr.ghxxu(indx,:); hxuu = dr.ghxxu(indx,:); -% gxss = dr.ghxss(indy,:); guss = dr.ghuss(indy,:); gxxx = dr.ghxxx(indy,:); guuu = dr.ghuuu(indy,:); gxxu = dr.ghxxu(indy,:); gxuu = dr.ghxxu(indy,:); -% Law of motion for first-order effects states xf and selected endogenous first-order effects variables yf: -% xf = hx*xf(-1) + hu*u -% yf = gx*xf(-1) + gu*u -% Law of motion for second-order effects states xs and selected endogenous second-order effects variables ys: -% xs = hx*xs(-1) + 1/2*hxx*kron(xf(-1),xf(-1)) + 1/2*huu*kron(u,u) + hxu*kron(xf(-1),u) + 1/2*hs2 -% ys = gx*xs(-1) + 1/2*gxx*kron(xf(-1),xf(-1)) + 1/2*guu*kron(u,u) + gxu*kron(xf(-1),u) + 1/2*gs2 -% Law of motion for third-order effects states xrd and selected endogenous second-order effects variables yrd: -% xrd = hx*xrd(-1) + hxx*kron(xf(-1),xs(-1)) + hxu*kron(xs(-1),u) + 3/6*hxss*xf(-1) + 3/6*huss*u + 1/6*hxxx*kron(xf(-1),kron(xf(-1),xf(-1))) + 1/6*huuu*kron(u,kron(u,u)) + 3/6*hxxu*kron(xf(-1),kron(xf(-1),u)) + 3/6*hxuu*kron(xf(-1),kron(u,u)) -% yrd = gx*xrd(-1) + gxx*kron(xf(-1),xs(-1)) + gxu*kron(xs(-1),u) + 3/6*gxss*xf(-1) + 3/6*guss*u + 1/6*gxxx*kron(xf(-1),kron(xf(-1),xf(-1))) + 1/6*guuu*kron(u,kron(u,u)) + 3/6*gxxu*kron(xf(-1),kron(xf(-1),u)) + 3/6*gxuu*kron(xf(-1),kron(u,u)) - % Selected variables: y = ybar + yf + ys - % Pruned state vector: z = [xf; xs; kron(xf,xf)] - % Pruned innovation vector: xi = [u; kron(u,u)-vec(Sigma_u); kron(xf(-1),u); kron(u,xf(-1))] - % Second moments of xi - % Sig_xi = E[u*u' , u*kron(u,u)' , zeros(u_nbr,x_nbr*u_nbr), zeros(u_nbr,u_nbr*x_nbr); - % kron(u,u)*u', kron(u,u)*kron(u,u)'-Sig_u(:)*Sig_u(:)', zeros(u_nbr^2,x_nbr*u_nbr), zeros(u_nbr^2,u_nbr*x_nbr); - % zeros(x_nbr*u_nbr,u_nbr), zeros(x_nbr*u_nbr,u_nbr^2), kron(xf(-1),u)*kron(xf(-1),u)', kron(xf(-1),u)*kron(u,xf(-1))'; - % zeros(u_nbr*x_nbr,u_nbr), zeros(u_nbr*x_nbr,u_nbr^2), kron(u,xf(-1))*kron(xf(-1),u)', kron(u,xf(-1))*kron(u,xf(-1))';] - % That is, we only need to compute: - % - Sig_xi_11 = E[u*u']=Sig_u - % - Sig_xi_21 = E[kron(u,u)*u']=0 for symmetric distributions - % - Sig_xi_12 = transpose(Sig_xi_21)= 0 for symmetric distributions - % - Sig_xi_22 = E[kron(u,u)*kron(u,u)']-Sig_u(:)*Sig_u(:)' which requires the fourth-order product moments of u - % - Sig_xi_34 = E[kron(xf(-1),u)*kron(xf(-1),u)'] which requires the second-order moments of xf and second-order moments of u - % - Sig_xi_43 = transpose(Sig_xi_34) - % - Sig_xi_33 = - % - Sig_xi_44 = +%% MAIN IDEA: +% Decompose the state vector x into first-order effects xf, second-order +% effects xs, and third-order effects xrd, i.e. x=xf+xs+xrd. Then, Dynare's +% perturbation approximation for the state vector up to third order +% (with Gaussian innovations u, i.e. no odd moments, hxxs=huus=hxus=hsss=0) is: +% x = hx*( xf(-1)+xs(-1)+xrd(-1) ) +% + hu*u +% + 1/2*hxx*kron( xf(-1)+xs(-1)+xrd(-1) , xf(-1)+xs(-1)+xrd(-1) ) +% + hxu*kron( xf(-1)+xs(-1)+xrd(-1) , u ) +% + 1/2*huu*kron( u , u ) +% + 1/2*hss*sig^2 +% + 1/6*hxxx*kron( kron( xf(-1)+xs(-1)+xrd(-1) , xf(-1)+xs(-1)+xrd(-1) ) , xf(-1)+xs(-1)+xrd(-1) ) +% + 1/6*huuu*kron( kron( u , u ) , u ) +% + 3/6*hxxu*kron( kron( xf(-1)+xs(-1)+xrd(-1) , xf(-1)+xs(-1)+xrd(-1) ) , u ) +% + 3/6*hxuu*kron( kron( xf(-1)+xs(-1)+xrd(-1) , u ) , u) +% + 3/6*hxss*( xf(-1)+xs(-1)+xrd(-1) )*sig^2 +% + 3/6*huss*u*sig^2 +% where: +% hx = dr.ghx(indx,:); hu = dr.ghu(indx,:); +% hxx = dr.ghxx(indx,:); hxu = dr.ghxu(indx,:); huu = dr.ghuu(indx,:); hss = dr.ghs2(indx,:); +% hxxx = dr.ghxxx(indx,:); huuu = dr.ghuuu(indx,:); hxxu = dr.ghxxu(indx,:); hxuu = dr.ghxxu(indx,:); hxss = dr.ghxss(indx,:); huss = dr.ghuss(indx,:); +% and similarly for control variables: +% y = gx*( xf(-1)+xs(-1)+xrd(-1) ) +% + gu*u +% + 1/2*gxx*kron( xf(-1)+xs(-1)+xrd(-1) , xf(-1)+xs(-1)+xrd(-1) ) +% + gxu*kron( xf(-1)+xs(-1)+xrd(-1) , u ) +% + 1/2*guu*kron( u , u ) +% + 1/2*gss*sig^2 +% + 1/6*gxxx*kron( kron( xf(-1)+xs(-1)+xrd(-1) , xf(-1)+xs(-1)+xrd(-1) ) , xf(-1)+xs(-1)+xrd(-1) ) +% + 1/6*guuu*kron( kron( u , u ) , u ) +% + 3/6*gxxu*kron( kron( xf(-1)+xs(-1)+xrd(-1) , xf(-1)+xs(-1)+xrd(-1) ) , u ) +% + 3/6*gxuu*kron( kron( xf(-1)+xs(-1)+xrd(-1) , u ) , u) +% + 3/6*gxss*( xf(-1)+xs(-1)+xrd(-1) )*sig^2 +% + 3/6*guss*u*sig^2 +% where: +% gx = dr.ghx(indy,:); gu = dr.ghu(indy,:); +% gxx = dr.ghxx(indy,:); gxu = dr.ghxu(indy,:); guu = dr.ghuu(indy,:); gss = dr.ghs2(indy,:); +% gxxx = dr.ghxxx(indy,:); guuu = dr.ghuuu(indy,:); gxxu = dr.ghxxu(indy,:); gxuu = dr.ghxxu(indy,:); gxss = dr.ghxss(indy,:); guss = dr.ghuss(indy,:); +% +% PRUNING means getting rid of terms higher than the approximation order, i.e. +% - involving fourth-order effects: kron(xf,xrd), kron(xs,xs), kron(xrd,xf), kron(xrd,u), +% kron(kron(xf,xf),xs), kron(kron(xf,xs),xf), kron(kron(xs,xf),xf) +% kron(kron(xf,xs),u), kron(kron(xs,xf),u) +% kron(kron(xs,u),u) +% xs*sig^2 +% - involving fifth-order effects: kron(xs,xrd), kron(xrd,xs), +% kron(kron(xf,xf),xrd), kron(kron(xf,xs),xs), kron(kron(xf,xrd),xf), kron(kron(xs,xf),xs), kron(kron(xs,xs),xf), kron(kron(xrd,xf),xf) +% kron(kron(xf,xrd),u), kron(kron(xs,xs),u), kron(kron(xrd,xf),u) +% kron(kron(xrd,u),u) +% xrd*sig^2 +% - involving sixth-order effects: kron(xrd,xrd), +% kron(kron(xf,xs),xrd), kron(kron(xf,xrd),xs), kron(kron(xs,xrd),xrd), kron(kron(xs,xs),xs), kron(kron(xs,xrd),xf), kron(kron(xf,xf),xs), kron(kron(xrd,xs),xf) +% kron(kron(xs,xrd),u), kron(kron(xrd,xs),u) +% - involving seventh-order effects: kron(kron(xf,xrd),xrd), kron(kron(xs,xs),xrd), kron(kron(xs,xrd),xs), kron(kron(xrd,xf),xrd), kron(kron(xrd,xs),xs), kron(kron(xrd,xrd),xf) +% kron(kron(xrd,xrd),u) +% - involving eighth-order effects: kron(kron(xs,xrd),xrd), kron(kron(xrd,xs),xrd), kron(kron(xrd,xrd),xs) +% - involving ninth-order effects: kron(kron(xrd,xrd),xrd) +% Note that u is treated as a first-order effect and the perturbation parameter sig as a variable. +% +% SUMMARY OF LAW OF MOTIONS: Set up the law of motions for the individual effects, but keep only effects of same order +% Notation: I_n=eye(n) and K_m_n=commutation(m,n) +% +% First-order effects: keep xf and u +% xf = hx*xf(-1) + hu*u +% Note that we +% +% Second-order effects: keep xs, kron(xf,xf), kron(u,u), kron(xf,u), and sig^2 +% xs = hx*xs(-1) + 1/2*hxx*kron(xf(-1),xf(-1)) + 1/2*huu*(kron(u,u)-Sigma_e(:)+Sigma_e(:)) + hxu*kron(xf(-1),u) + 1/2*hss*sig^2% +% +% Third-order effects: keep xrd, kron(xf,xs), kron(xs,xf), kron(xs,u), kron(kron(xf,xf),xf), kron(kron(u,u),u), kron(kron(xf,xf),u), kron(kron(xf,u),u), xf*sig^2, u*sig^2 +% xrd = hx*xrd(-1) + 1/2*hxx*(kron(xf(-1),xs(-1))+kron(xs(-1),xf(-1))) + hxu*kron(xs(-1),u) + 1/6*hxxx*kron(xf(-1),kron(xf(-1),xf(-1))) + 1/6*huuu*kron(u,kron(u,u)) + 3/6*hxxu*kron(xf(-1),kron(xf(-1),u)) + 3/6*hxuu*kron(xf(-1),kron(u,u)) + 3/6*hxss*xf(-1)*sig^2 + 3/6*huss*u*sig^2 +% Simplified (due to symmetry in hxx): +% xrd = hx*xrd(-1) + hxx*(kron(xf(-1),xs(-1)) + hxu*kron(xs(-1),u) + 1/6*hxxx*kron(xf(-1),kron(xf(-1),xf(-1))) + 1/6*huuu*kron(u,kron(u,u)) + 3/6*hxxu*kron(xf(-1),kron(xf(-1),u)) + 3/6*hxuu*kron(xf(-1),kron(u,u)) + 3/6*hxss*xf(-1)*sig^2 + 3/6*huss*u*sig^2% +% +% Auxiliary equation kron(xf,xf) to set up the VAR(1) pruned state space system +% kron(xf,xf) = kron(hx,hx)*kron(xf(-1),xf(-1)) + kron(hu,hu)*(kron(u,u)-Sigma_e(:)+Sigma_e(:)) + kron(hx,u)*kron(xf(-1),u) + kron(u,hx)*kron(u,xf(-1)) +% Simplified using commutation matrix: +% kron(xf,xf) = kron(hx,hx)*kron(xf(-1),xf(-1)) + (I_xx+K_x_x)*kron(hx,hu)*kron(xf(-1),u) + kron(hu,hu)*kron(u,u) +% +% Auxiliary equation kron(xf,xs) to set up the VAR(1) pruned state space system +% kron(xf,xs) = kron(hx,hx)*kron(xf(-1),xs(-1)) + kron(hu,hx)*kron(u,xs(-1)) +% + kron(hx,1/2*hxx)*kron(kron(xf(-1),xf(-1)),xf(-1)) + kron(hu,1/2*hxx)*kron(kron(u,xf(-1)),xf(-1)) +% + kron(hx,1/2*huu)*kron(kron(xf(-1),u),u) + kron(hu,1/2*huu)*kron(kron(u,u),u) +% + kron(hx,hxu)*kron(kron(xf(-1),xf(-1)),u) + kron(hu,hxu)*kron(kron(u,xf(-1)),u) +% + kron(hx,1/2*hss)*xf(-1)*sig^2 + kron(hu,1/2*hss)*u*sig^2 +% Simplified using commutation matrix: +% kron(xf,xs) = kron(hx,hx)*kron(xf(-1),xs(-1)) +% + K_x_x*kron(hx,hu)*kron(xs(-1),u) +% + kron(hx,1/2*hxx)*kron(kron(xf(-1),xf(-1)),xf(-1)) +% + ( kron(hx,hxu) + K_x_x*kron(1/2*hxx,hu) )*kron(kron(xf(-1),xf(-1)),u) +% + ( kron(hx,1/2*huu) + K_x_x*kron(hxu,hu) )*kron(kron(xf(-1),u),u) +% + kron(hu,1/2*huu)*kron(kron(u,u),u) +% + kron(hx,1/2*hss)*xf(-1)*sig^2 +% + kron(hu,1/2*hss)*u*sig^2 +% +% Auxiliary equation kron(kron(xf,xf),xf) to set up the VAR(1) pruned state space system +% kron(kron(xf,xf),xf) = kron(kron(hx,hx),hx)*kron(kron(xf(-1),xf(-1)),xf(-1)) +% + kron(kron(hx,hu),hx)*kron(kron(xf(-1),u),xf(-1)) +% + kron(kron(hx,hx),hu)*kron(kron(xf(-1),xf(-1)),u) +% + kron(kron(hx,hu),hu)*kron(kron(xf(-1),u),u) +% + kron(kron(hu,hx),hx)*kron(kron(u,xf(-1)),xf(-1)) +% + kron(kron(hu,hu),hx)*kron(kron(u,u),xf(-1)) +% + kron(kron(hu,hx),hu)*kron(kron(u,xf(-1)),u) +% + kron(kron(hu,hu),hu)*kron(kron(u,u),u) +% Simplified using commutation matrix: +% kron(kron(xf,xf),xf) = kron(kron(hx,hx),hx)*kron(kron(xf(-1),xf(-1)),xf(-1)) +% + ( kron(kron(hx,hx),hu) + K_xx_x*kron(hx,(I_xx+K_x_x)*kron(hx,hu)) )*kron(kron(xf(-1),xf(-1)),u) +% + ( kron((I_xx+K_x_x)*kron(hx,hu),hu) + K_xx_x*kron(kron(hx,hu),hu) )*kron(kron(xf(-1),u),u) +% + kron(kron(hu,hu),hu)*kron(kron(u,u),u) +% +% Law of motion for control variables y (either VAROBS variables or if no VAROBS statement is given then for all endogenous variables) +% y = steady_state(y) +% + gx*( xf(-1)+xs(-1)+xrd(-1) ) +% + gu*u +% + 1/2*gxx*kron(xf(-1),xf(-1)) + gxx*kron(xf(-1),xs(-1)) +% + gxu*kron(xf(-1),u) + gxu*kron(xs(-1),u) +% + 1/2*guu*(kron(u,u)-Sigma_e+Sigma_e) +% + 1/2*gss*sig^2 +% + 1/6*gxxx*kron(kron(xf(-1),xf(-1)),xf(-1)) +% + 1/6*guuu*kron(kron(u,u),u) +% + 3/6*gxxu*kron(kron(xf(-1),xf(-1),u) +% + 3/6*gxuu*kron(kron(xf(-1),u),u) +% + 3/6*gxss*xf(-1)*sig^2 +% + 3/6*guss*u*sig^2 +% +% See code below how z and inov are defined at first, second, and third order, +% and how to set up A, B, C, D and compute unconditional first and second moments of inov, z and y -%% Auxiliary indices, numbers and flags +%% Auxiliary indices and objects 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 + % Note that unit roots are only possible at first-order, at higher order we set it to 1 + options.qz_criterium = 1+1e-6; +end indx = [M.nstatic+(1:M.nspred) M.endo_nbr+(1:size(dr.ghx,2)-M.nspred)]'; -indy = (1:M.endo_nbr)'; %by default select all variables -compute_derivs = false; -if isfield(options,'options_ident') - compute_derivs = true; - stderrparam_nbr = length(options.options_ident.indpstderr); - corrparam_nbr = size(options.options_ident.indpcorr,1); - modparam_nbr = length(options.options_ident.indpmodel); +if isempty(indy) + indy = (1:M.endo_nbr)'; %by default select all variables in DR order +end +u_nbr = M.exo_nbr; +x_nbr = length(indx); +y_nbr = length(indy); +Yss = dr.ys(dr.order_var); +hx = dr.ghx(indx,:); +gx = dr.ghx(indy,:); +hu = dr.ghu(indx,:); +gu = dr.ghu(indy,:); +E_uu = M.Sigma_e; %this is E[u*u'] + +if compute_derivs + stderrparam_nbr = length(dr.derivs.indpstderr); + corrparam_nbr = size(dr.derivs.indpcorr,1); + modparam_nbr = length(dr.derivs.indpmodel); totparam_nbr = stderrparam_nbr+corrparam_nbr+modparam_nbr; - indy = options.options_ident.indvobs; -end -Yss = dr.ys(dr.order_var); if compute_derivs; dYss = dr.derivs.dYss; end -E_uu = M.Sigma_e; if compute_derivs; dE_uu = dr.derivs.dSigma_e; end -Correlation_matrix = M.Correlation_matrix; if compute_derivs; dCorrelation_matrix = dr.derivs.dCorrelation_matrix; end -ghx = dr.ghx; if compute_derivs; dghx = dr.derivs.dghx; end -ghu = dr.ghu; if compute_derivs; dghu = dr.derivs.dghu; end -Om = ghu*E_uu*transpose(ghu); if compute_derivs; dOm = dr.derivs.dOm; end -y_nbr = length(indy); -x_nbr = length(indx); -u_nbr = M.exo_nbr; - -% indices for extended state vector z and extended shock vector e -id_z1_xf = (1:x_nbr); -id_e1_u = (1:u_nbr); -if options.order > 1 - id_z2_xs = id_z1_xf(end) + (1:x_nbr); - id_z3_xf_xf = id_z2_xs(end) + (1:x_nbr*x_nbr); - id_e2_u_u = id_e1_u(end) + (1:u_nbr*u_nbr); - id_e3_xf_u = id_e2_u_u(end) + (1:x_nbr*u_nbr); - id_e4_u_xf = id_e3_xf_u(end) + (1:u_nbr*x_nbr); - - ghxx = dr.ghxx; if compute_derivs; dghxx = dr.derivs.dghxx; end - ghxu = dr.ghxu; if compute_derivs; dghxu = dr.derivs.dghxu; end - ghuu = dr.ghuu; if compute_derivs; dghuu = dr.derivs.dghuu; end - ghs2 = dr.ghs2; if compute_derivs; dghs2 = dr.derivs.dghs2; end -end -if options.order > 2 - id_z4_xrd = id_z3_xf_xf(end) + (1:x_nbr); - id_z5_xf_xs = id_z4_xrd(end) + (1:x_nbr*x_nbr); - id_z6_xf_xf_xf = id_z5_xf_xs(end) + (1:x_nbr*x_nbr*x_nbr); - id_e5_xs_u = id_e4_u_xf(end) + (1:x_nbr*u_nbr); - id_e6_u_xs = id_e5_xs_u(end) + (1:u_nbr*x_nbr); - id_e7_xf_xf_u = id_e6_u_xs(end) + (1:x_nbr*x_nbr*u_nbr); - id_e8_xf_u_xf = id_e7_xf_xf_u(end) + (1:x_nbr*u_nbr*x_nbr); - id_e9_u_xf_xf = id_e8_xf_u_xf(end) + (1:u_nbr*x_nbr*x_nbr); - id_e10_xf_u_u = id_e9_u_xf_xf(end) + (1:x_nbr*u_nbr*u_nbr); - id_e11_u_xf_u = id_e10_xf_u_u(end) + (1:u_nbr*x_nbr*u_nbr); - id_e12_u_u_xf = id_e11_u_xf_u(end) + (1:u_nbr*u_nbr*x_nbr); - id_e13_u_u_u = id_e12_u_u_xf(end) + (1:u_nbr*u_nbr*u_nbr); - - ghxxx = dr.ghxxx; if compute_derivs; dghxxx = dr.derivs.dghxxx; end - ghxxu = dr.ghxxu; if compute_derivs; dghxxu = dr.derivs.dghxxu; end - ghxuu = dr.ghxuu; if compute_derivs; dghxuu = dr.derivs.dghxuu; end - ghuuu = dr.ghuuu; if compute_derivs; dghuuu = dr.derivs.dghuuu; end - ghxss = dr.ghxss; if compute_derivs; dghxss = dr.derivs.dghxss; end - ghuss = dr.ghuss; if compute_derivs; dghuss = dr.derivs.dghuss; end + dYss = dr.derivs.dYss; + dhx = dr.derivs.dghx(indx,:,:); + dgx = dr.derivs.dghx(indy,:,:); + dhu = dr.derivs.dghu(indx,:,:); + dgu = dr.derivs.dghu(indy,:,:); + dE_uu = dr.derivs.dSigma_e; +end + +% first-order approximation indices for extended state vector z and extended innovations vector inov +id_z1_xf = (1:x_nbr); +id_inov1_u = (1:u_nbr); +if order > 1 + % second-order approximation indices for extended state vector z and extended innovations vector inov + id_z2_xs = id_z1_xf(end) + (1:x_nbr); + id_z3_xf_xf = id_z2_xs(end) + (1:x_nbr*x_nbr); + id_inov2_u_u = id_inov1_u(end) + (1:u_nbr*u_nbr); + id_inov3_xf_u = id_inov2_u_u(end) + (1:x_nbr*u_nbr); + + hxx = dr.ghxx(indx,:); + gxx = dr.ghxx(indy,:); + hxu = dr.ghxu(indx,:); + gxu = dr.ghxu(indy,:); + huu = dr.ghuu(indx,:); + guu = dr.ghuu(indy,:); + hss = dr.ghs2(indx,:); + gss = dr.ghs2(indy,:); + if compute_derivs + dhxx = dr.derivs.dghxx(indx,:,:); + dgxx = dr.derivs.dghxx(indy,:,:); + dhxu = dr.derivs.dghxu(indx,:,:); + dgxu = dr.derivs.dghxu(indy,:,:); + dhuu = dr.derivs.dghuu(indx,:,:); + dguu = dr.derivs.dghuu(indy,:,:); + dhss = dr.derivs.dghs2(indx,:); + dgss = dr.derivs.dghs2(indy,:); + end +end +if order > 2 + % third-order approximation indices for extended state vector z and extended innovations vector inov + id_z4_xrd = id_z3_xf_xf(end) + (1:x_nbr); + id_z5_xf_xs = id_z4_xrd(end) + (1:x_nbr*x_nbr); + id_z6_xf_xf_xf = id_z5_xf_xs(end) + (1:x_nbr*x_nbr*x_nbr); + id_inov4_xs_u = id_inov3_xf_u(end) + (1:x_nbr*u_nbr); + id_inov5_xf_xf_u = id_inov4_xs_u(end) + (1:x_nbr*x_nbr*u_nbr); + id_inov6_xf_u_u = id_inov5_xf_xf_u(end) + (1:x_nbr*u_nbr*u_nbr); + id_inov7_u_u_u = id_inov6_xf_u_u(end) + (1:u_nbr*u_nbr*u_nbr); + + hxxx = dr.ghxxx(indx,:); + gxxx = dr.ghxxx(indy,:); + hxxu = dr.ghxxu(indx,:); + gxxu = dr.ghxxu(indy,:); + hxuu = dr.ghxuu(indx,:); + gxuu = dr.ghxuu(indy,:); + huuu = dr.ghuuu(indx,:); + guuu = dr.ghuuu(indy,:); + hxss = dr.ghxss(indx,:); + gxss = dr.ghxss(indy,:); + huss = dr.ghuss(indx,:); + guss = dr.ghuss(indy,:); + if compute_derivs + dhxxx = dr.derivs.dghxxx(indx,:,:); + dgxxx = dr.derivs.dghxxx(indy,:,:); + dhxxu = dr.derivs.dghxxu(indx,:,:); + dgxxu = dr.derivs.dghxxu(indy,:,:); + dhxuu = dr.derivs.dghxuu(indx,:,:); + dgxuu = dr.derivs.dghxuu(indy,:,:); + dhuuu = dr.derivs.dghuuu(indx,:,:); + dguuu = dr.derivs.dghuuu(indy,:,:); + dhxss = dr.derivs.dghxss(indx,:,:); + dgxss = dr.derivs.dghxss(indy,:,:); + dhuss = dr.derivs.dghuss(indx,:,:); + dguss = dr.derivs.dghuss(indy,:,:); + end +end + +%% First-order state space system +% Auxiliary state vector z is defined by: z = [xf] +% Auxiliary innovations vector inov is defined by: inov = [u] +z_nbr = x_nbr; +inov_nbr = M.exo_nbr; +A = hx; +B = hu; +c = zeros(x_nbr,1); +C = gx; +D = gu; +d = zeros(y_nbr,1); +Varinov = E_uu; +E_inovzlag1 = zeros(inov_nbr,z_nbr); %at first-order E[inov*z(-1)'] = 0 +Om_z = B*Varinov*B'; +E_xf = zeros(x_nbr,1); + +lyapunov_symm_method = 1; %method=1 to initialize persistent variables +[Var_z,Schur_u] = lyapunov_symm(A, Om_z,... %at first-order this algorithm is well established and also used in th_autocovariances.m + options.lyapunov_fixed_point_tol, options.qz_criterium, options.lyapunov_complex_threshold,... + lyapunov_symm_method,... + options.debug); %we use Schur_u to take care of (possible) nonstationary VAROBS variables in moment computations +%find stationary vars +stationary_vars = (1:y_nbr)'; +if ~isempty(Schur_u) + %base this only on first order, because if first-order is stable so are the higher-order pruned systems + x = abs(gx*Schur_u); + stationary_vars = find(all(x < options.Schur_vec_tol,2)); end -%% First-order -z_nbr = x_nbr; -e_nbr = M.exo_nbr; -A = ghx(indx,:); -B = ghu(indx,:); -C = ghx(indy,:); -D = ghu(indy,:); -c = zeros(x_nbr,1); -d = zeros(y_nbr,1); -Varinov = E_uu; -E_z = zeros(x_nbr,1); %this is E_xf if compute_derivs == 1 - dA = dghx(indx,:,:); - dB = dghu(indx,:,:); - dC = dghx(indy,:,:); - dD = dghu(indy,:,:); - dc = zeros(x_nbr,totparam_nbr); - dd = zeros(y_nbr,totparam_nbr); - dVarinov = dE_uu; - dE_z = zeros(x_nbr,totparam_nbr); + dA = dhx; + dB = dhu; + dc = zeros(x_nbr,totparam_nbr); + dC = dgx; + dD = dgu; + dd = zeros(y_nbr,totparam_nbr); + dVarinov = dE_uu; + dE_xf = zeros(x_nbr,totparam_nbr); + dE_inovzlag1 = zeros(z_nbr,inov_nbr,totparam_nbr); + dVar_z = zeros(z_nbr,z_nbr,totparam_nbr); + lyapunov_symm_method = 2;%to spare a lot of computing time while not repeating Schur every time + for jp1 = 1:totparam_nbr + if jp1 <= stderrparam_nbr+corrparam_nbr + dOm_z_jp1 = B*dVarinov(:,:,jp1)*B'; + dVar_z(:,:,jp1) = lyapunov_symm(A, dOm_z_jp1,... + options.lyapunov_fixed_point_tol,options.qz_criterium,options.lyapunov_complex_threshold,... + lyapunov_symm_method,... + options.debug); + else + dOm_z_jp1 = dB(:,:,jp1)*Varinov*B' + B*Varinov*dB(:,:,jp1)'; + dVar_z(:,:,jp1) = lyapunov_symm(A, dA(:,:,jp1)*Var_z*A' + A*Var_z*dA(:,:,jp1)' + dOm_z_jp1,... + options.lyapunov_fixed_point_tol,options.qz_criterium,options.lyapunov_complex_threshold,... + lyapunov_symm_method,... + options.debug); + end + end end if order > 1 + options.qz_criterium = 1; %pruned state space has no unit roots % Some common and useful objects for order > 1 - % Compute E[xf*xf'] - E_xfxf = lyapunov_symm(ghx(indx,:), Om(indx,indx), options.lyapunov_fixed_point_tol, options.qz_criterium, options.lyapunov_complex_threshold, 1, options.debug); %use 1 to initialize persistent variables - K_ux = commutation(u_nbr,x_nbr); %commutation matrix - K_xu = commutation(x_nbr,u_nbr); %commutation matrix - QPu = quadruplication(u_nbr); %quadruplication matrix as in Meijer (2005), similar to Magnus-Neudecker definition of duplication matrix (i.e. dyn_vec) but for fourth-order product moments - %Compute unique product moments of E[kron(u*u',uu')] = E[kron(u,u)*kron(u,u)'] - COMBOS4 = flipud(allVL1(u_nbr, 4)); %all possible combinations of powers that sum up to four for fourth-order product moments of u - E_u_u_u_u = zeros(u_nbr*(u_nbr+1)/2*(u_nbr+2)/3*(u_nbr+3)/4,1); %only unique entries of E[kron(u,kron(u,kron(u,u)))] - if compute_derivs && (stderrparam_nbr+corrparam_nbr>0) - dE_u_u_u_u = zeros(u_nbr*(u_nbr+1)/2*(u_nbr+2)/3*(u_nbr+3)/4,stderrparam_nbr+corrparam_nbr); + E_xfxf = Var_z; + if compute_derivs + dE_xfxf = dVar_z; end - for j = 1:size(COMBOS4,1) - E_u_u_u_u(j) = prodmom(E_uu, 1:u_nbr, COMBOS4(j,:)); + hx_hx = kron(hx,hx); + hx_hu = kron(hx,hu); + hu_hu = kron(hu,hu); + I_xx = eye(x_nbr^2); + K_x_x = commutation(x_nbr,x_nbr,1); + invIx_hx = (eye(x_nbr)-hx)\eye(x_nbr); + + %Compute unique fourth order product moments of u, i.e. unique(E[kron(kron(kron(u,u),u),u)],'stable') + u_nbr4 = u_nbr*(u_nbr+1)/2*(u_nbr+2)/3*(u_nbr+3)/4; + QPu = quadruplication(u_nbr); + COMBOS4 = flipud(allVL1(u_nbr, 4)); %all possible (unique) combinations of powers that sum up to four + E_u_u_u_u = zeros(u_nbr4,1); %only unique entries + if compute_derivs && (stderrparam_nbr+corrparam_nbr>0) + dE_u_u_u_u = zeros(u_nbr4,stderrparam_nbr+corrparam_nbr); + end + for j4 = 1:size(COMBOS4,1) if compute_derivs && (stderrparam_nbr+corrparam_nbr>0) - dE_u_u_u_u(j,1:(stderrparam_nbr+corrparam_nbr)) = prodmom_deriv(E_uu, 1:u_nbr, COMBOS4(j,:), dE_uu(:,:,1:(stderrparam_nbr+corrparam_nbr)), dCorrelation_matrix(:,:,1:(stderrparam_nbr+corrparam_nbr))); + [E_u_u_u_u(j4), dE_u_u_u_u(j4,:)] = prodmom_deriv(E_uu, 1:u_nbr, COMBOS4(j4,:), dE_uu(:,:,1:(stderrparam_nbr+corrparam_nbr)), dr.derivs.dCorrelation_matrix(:,:,1:(stderrparam_nbr+corrparam_nbr))); + else + E_u_u_u_u(j4) = prodmom(E_uu, 1:u_nbr, COMBOS4(j4,:)); end end - E_u_u_u_u = QPu*E_u_u_u_u; %add duplicate product moments, i.e. this is E[kron(u,kron(u,kron(u,u)))] - E_uu_uu = reshape(E_u_u_u_u,u_nbr^2,u_nbr^2); %E[kron(u*u',uu')] = E[kron(u,u)*kron(u,u)'] + E_xfxf_uu = kron(E_xfxf,E_uu'); + +%% Second-order pruned state space system +% Auxiliary state vector z is defined by: z = [xf;xs;kron(xf,xf)] +% Auxiliary innovations vector inov is defined by: inov = [u;kron(u,u)-E_uu(:);kron(xf,u)] + z_nbr = x_nbr + x_nbr + x_nbr^2; + inov_nbr = u_nbr + u_nbr^2 + x_nbr*u_nbr; - % E[kron((xf*xf'),(u*u'))] - E_xfxf_E_uu = kron(E_xfxf,E_uu); - if compute_derivs - dE_xfxf = zeros(size(E_xfxf,1),size(E_xfxf,2),totparam_nbr); - dE_xfxf_E_uu = zeros(size(E_xfxf_E_uu,1),size(E_xfxf_E_uu,2),totparam_nbr); - for jp = 1:totparam_nbr - if jp <= (stderrparam_nbr+corrparam_nbr) - %Jacobian of E_xfxf wrt exogenous paramters: dE_xfxf(:,:,jp)-hx*dE_xfxf(:,:,jp)*hx'=dOm(:,:,jp), because dhx is zero by construction for stderr and corr parameters - dE_xfxf(:,:,jp) = lyapunov_symm(ghx(indx,:), dOm(indx,indx,jp),options.lyapunov_fixed_point_tol,options.qz_criterium,options.lyapunov_complex_threshold,2,options.debug); - else - %Jacobian of E_xfxf wrt model parameters: dE_xfxf(:,:,jp) - hx*dE_xfxf(:,:,jp)*hx' = d(hu*Sig_u*hu')(:,:,jp) + dhx(:,:,jp)*E_xfxf*hx'+ hx*E_xfxf*dhx(:,:,jp)' - dE_xfxf(:,:,jp) = lyapunov_symm(ghx(indx,:), dghx(indx,:,jp)*E_xfxf*ghx(indx,:)'+ghx(indx,:)*E_xfxf*dghx(indx,:,jp)'+dOm(indx,indx,jp),options.lyapunov_fixed_point_tol,options.qz_criterium,options.lyapunov_complex_threshold,2,options.debug); - %here method=2 is used to spare a lot of computing time while not repeating Schur every time - end - dE_xfxf_E_uu(:,:,jp) = kron(dE_xfxf(:,:,jp),E_uu) + kron(E_xfxf,dE_uu(:,:,jp)); - end - end - - % Second-order pruned state space system - z_nbr = x_nbr+x_nbr+x_nbr*x_nbr; - e_nbr = u_nbr+u_nbr*u_nbr+x_nbr*u_nbr+u_nbr*x_nbr; - A = zeros(z_nbr, z_nbr); - A(id_z1_xf , id_z1_xf ) = ghx(indx,:); - A(id_z2_xs , id_z2_xs ) = ghx(indx,:); - A(id_z2_xs , id_z3_xf_xf) = 0.5*ghxx(indx,:); - A(id_z3_xf_xf , id_z3_xf_xf) = kron(ghx(indx,:),ghx(indx,:)); - - B = zeros(z_nbr, e_nbr); - B(id_z1_xf , id_e1_u ) = ghu(indx,:); - B(id_z2_xs , id_e2_u_u ) = 0.5*ghuu(indx,:); - B(id_z2_xs , id_e3_xf_u) = ghxu(indx,:); - B(id_z3_xf_xf , id_e2_u_u ) = kron(ghu(indx,:),ghu(indx,:)); - B(id_z3_xf_xf , id_e3_xf_u) = kron(ghx(indx,:),ghu(indx,:)); - B(id_z3_xf_xf , id_e4_u_xf) = kron(ghu(indx,:),ghx(indx,:)); - - C = zeros(y_nbr, z_nbr); - C(1:y_nbr , id_z1_xf ) = ghx(indy,:); - C(1:y_nbr , id_z2_xs ) = ghx(indy,:); - C(1:y_nbr , id_z3_xf_xf) = 0.5*ghxx(indy,:); - - D = zeros(y_nbr, e_nbr); - D(1:y_nbr , id_e1_u ) = ghu(indy,:); - D(1:y_nbr , id_e2_u_u ) = 0.5*ghuu(indy,:); - D(1:y_nbr , id_e3_xf_u) = ghxu(indy,:); - + A(id_z1_xf , id_z1_xf ) = hx; + A(id_z2_xs , id_z2_xs ) = hx; + A(id_z2_xs , id_z3_xf_xf) = 1/2*hxx; + A(id_z3_xf_xf , id_z3_xf_xf) = hx_hx; + + B = zeros(z_nbr, inov_nbr); + B(id_z1_xf , id_inov1_u ) = hu; + B(id_z2_xs , id_inov2_u_u ) = 1/2*huu; + B(id_z2_xs , id_inov3_xf_u) = hxu; + B(id_z3_xf_xf , id_inov2_u_u ) = hu_hu; + B(id_z3_xf_xf , id_inov3_xf_u) = (I_xx+K_x_x)*hx_hu; + c = zeros(z_nbr, 1); - c(id_z2_xs , 1) = 0.5*ghs2(indx,:) + 0.5*ghuu(indx,:)*E_uu(:); - c(id_z3_xf_xf , 1) = kron(ghu(indx,:),ghu(indx,:))*E_uu(:); - - d = zeros(y_nbr, 1); - d(1:y_nbr, 1) = 0.5*ghs2(indy,:) + 0.5*ghuu(indy,:)*E_uu(:); + c(id_z2_xs , 1) = 1/2*hss + 1/2*huu*E_uu(:); + c(id_z3_xf_xf , 1) = hu_hu*E_uu(:); - Varinov = zeros(e_nbr,e_nbr); - Varinov(id_e1_u , id_e1_u) = E_uu; - Varinov(id_e2_u_u , id_e2_u_u) = E_uu_uu-E_uu(:)*transpose(E_uu(:)); - Varinov(id_e3_xf_u , id_e3_xf_u) = E_xfxf_E_uu; - Varinov(id_e4_u_xf , id_e3_xf_u) = K_ux*E_xfxf_E_uu; - Varinov(id_e3_xf_u , id_e4_u_xf) = transpose(Varinov(id_e4_u_xf,id_e3_xf_u)); - Varinov(id_e4_u_xf , id_e4_u_xf) = K_ux*E_xfxf_E_uu*transpose(K_ux); - - I_hx = speye(x_nbr)-ghx(indx,:); - I_hxinv = I_hx\speye(x_nbr); - E_xs = I_hxinv*(0.5*ghxx(indx,:)*E_xfxf(:) + c(x_nbr+(1:x_nbr),1)); - - E_z = [zeros(x_nbr,1); E_xs; E_xfxf(:)]; - if compute_derivs - dA = zeros(size(A,1),size(A,2),totparam_nbr); - dB = zeros(size(B,1),size(B,2),totparam_nbr); - dC = zeros(size(C,1),size(C,2),totparam_nbr); - dD = zeros(size(D,1),size(D,2),totparam_nbr); - dc = zeros(size(c,1),totparam_nbr); - dd = zeros(size(d,1),totparam_nbr); - dVarinov = zeros(size(Varinov,1),size(Varinov,2),totparam_nbr); - dE_xs = zeros(size(E_xs,1),size(E_xs,2),totparam_nbr); - dE_z = zeros(size(E_z,1),size(E_z,2),totparam_nbr); - for jp = 1:totparam_nbr - dA(id_z1_xf , id_z1_xf , jp) = dghx(indx,:,jp); - dA(id_z2_xs , id_z2_xs , jp) = dghx(indx,:,jp); - dA(id_z2_xs , id_z3_xf_xf , jp) = 0.5*dghxx(indx,:,jp); - dA(id_z3_xf_xf , id_z3_xf_xf , jp) = kron(dghx(indx,:,jp),ghx(indx,:)) + kron(ghx(indx,:),dghx(indx,:,jp)); + C = zeros(y_nbr, z_nbr); + C(: , id_z1_xf ) = gx; + C(: , id_z2_xs ) = gx; + C(: , id_z3_xf_xf) = 1/2*gxx; - dB(id_z1_xf , id_e1_u , jp) = dghu(indx,:,jp); - dB(id_z2_xs , id_e2_u_u , jp) = 0.5*dghuu(indx,:,jp); - dB(id_z2_xs , id_e3_xf_u , jp) = dghxu(indx,:,jp); - dB(id_z3_xf_xf , id_e2_u_u , jp) = kron(dghu(indx,:,jp),ghu(indx,:)) + kron(ghu(indx,:),dghu(indx,:,jp)); - dB(id_z3_xf_xf , id_e3_xf_u , jp) = kron(dghx(indx,:,jp),ghu(indx,:)) + kron(ghx(indx,:),dghu(indx,:,jp)); - dB(id_z3_xf_xf , id_e4_u_xf , jp) = kron(dghu(indx,:,jp),ghx(indx,:)) + kron(ghu(indx,:),dghx(indx,:,jp)); + D = zeros(y_nbr, inov_nbr); + D(: , id_inov1_u ) = gu; + D(: , id_inov2_u_u ) = 1/2*guu; + D(: , id_inov3_xf_u) = gxu; - dC(1:y_nbr , id_z1_xf , jp) = dghx(indy,:,jp); - dC(1:y_nbr , id_z2_xs , jp) = dghx(indy,:,jp); - dC(1:y_nbr , id_z3_xf_xf , jp) = 0.5*dghxx(indy,:,jp); + d = 1/2*gss + 1/2*guu*E_uu(:); - dD(1:y_nbr , id_e1_u , jp) = dghu(indy,:,jp); - dD(1:y_nbr , id_e2_u_u , jp) = 0.5*dghuu(indy,:,jp); - dD(1:y_nbr , id_e3_xf_u , jp) = dghxu(indy,:,jp); + Varinov = zeros(inov_nbr,inov_nbr); + Varinov(id_inov1_u , id_inov1_u) = E_uu; + %Varinov(id_inov1_u , id_inov2_u_u ) = zeros(u_nbr,u_nbr^2); + %Varinov(id_inov1_u , id_inov3_xf_u) = zeros(u_nbr,x_nbr*u_nbr); + %Varinov(id_inov2_u_u , id_inov1_u ) = zeros(u_nbr^2,u_nbr); + Varinov(id_inov2_u_u , id_inov2_u_u ) = reshape(QPu*E_u_u_u_u,u_nbr^2,u_nbr^2)-E_uu(:)*E_uu(:)'; + %Varinov(id_inov2_u_u , id_inov3_xf_u) = zeros(u_nbr^2,x_nbr*u_nbr); + %Varinov(id_inov3_xf_u , id_inov1_u ) = zeros(x_nbr*u_nbr,u_nbr); + %Varinov(id_inov3_xf_u , id_inov2_u_u ) = zeros(x_nbr*u_nbr,u_nbr^2); + Varinov(id_inov3_xf_u , id_inov3_xf_u) = E_xfxf_uu; - dc(id_z2_xs , jp) = 0.5*dghs2(indx,jp) + 0.5*(dghuu(indx,:,jp)*E_uu(:) + ghuu(indx,:)*vec(dE_uu(:,:,jp))); - dc(id_z3_xf_xf , jp) = (kron(dghu(indx,:,jp),ghu(indx,:)) + kron(ghu(indx,:),dghu(indx,:,jp)))*E_uu(:) + kron(ghu(indx,:),ghu(indx,:))*vec(dE_uu(:,:,jp)); + E_xs = invIx_hx*(1/2*hxx*E_xfxf(:) + c(id_z2_xs,1)); + E_inovzlag1 = zeros(inov_nbr,z_nbr); %at second-order E[z(-1)*inov'] = 0 + Om_z = B*Varinov*transpose(B); - dd(1:y_nbr , jp) = 0.5*dghs2(indy,jp) + 0.5*(dghuu(indy,:,jp)*E_uu(:) + ghuu(indy,:)*vec(dE_uu(:,:,jp))); + lyapunov_symm_method = 1; %method=1 to initialize persistent variables (if errorflag) + [Var_z, errorflag] = disclyap_fast(A,Om_z,options.lyapunov_doubling_tol); + if errorflag %use Schur-based method + fprintf('PRUNED_STATE_SPACE_SYSTEM: error flag in disclyap_fast at order=2, use lyapunov_symm\n'); + Var_z = lyapunov_symm(A,Om_z,... + options.lyapunov_fixed_point_tol,options.qz_criterium,options.lyapunov_complex_threshold,... + lyapunov_symm_method,... + options.debug); + lyapunov_symm_method = 2; %in the following we can reuse persistent variables + end + % Make sure some stuff is zero due to Gaussianity + Var_z(id_z1_xf , id_z2_xs ) = zeros(x_nbr,x_nbr); + Var_z(id_z1_xf , id_z3_xf_xf) = zeros(x_nbr,x_nbr^2); + Var_z(id_z2_xs , id_z1_xf ) = zeros(x_nbr,x_nbr); + Var_z(id_z3_xf_xf , id_z1_xf ) = zeros(x_nbr^2,x_nbr); - if jp <= (stderrparam_nbr+corrparam_nbr) - dVarinov(id_e1_u , id_e1_u , jp) = dE_uu(:,:,jp); - dVarinov(id_e2_u_u , id_e2_u_u , jp) = reshape(QPu*dE_u_u_u_u(:,jp), u_nbr^2, u_nbr^2) - (reshape(dE_uu(:,:,jp),u_nbr^2,1)*transpose(E_uu(:)) + E_uu(:)*transpose(reshape(dE_uu(:,:,jp),u_nbr^2,1))); + if compute_derivs + dA = zeros(z_nbr,z_nbr,totparam_nbr); + dB = zeros(z_nbr,inov_nbr,totparam_nbr); + dc = zeros(z_nbr,totparam_nbr); + dC = zeros(y_nbr,z_nbr,totparam_nbr); + dD = zeros(y_nbr,inov_nbr,totparam_nbr); + dd = zeros(y_nbr,totparam_nbr); + dVarinov = zeros(inov_nbr,inov_nbr,totparam_nbr); + dE_xs = zeros(x_nbr,totparam_nbr); + dE_inovzlag1 = zeros(inov_nbr,z_nbr,totparam_nbr); + dVar_z = zeros(z_nbr,z_nbr,totparam_nbr); + + for jp2 = 1:totparam_nbr + if jp2 <= (stderrparam_nbr+corrparam_nbr) + dE_uu_jp2 = dE_uu(:,:,jp2); + dE_u_u_u_u_jp2 = QPu*dE_u_u_u_u(:,jp2); + else + dE_uu_jp2 = zeros(u_nbr,u_nbr); + dE_u_u_u_u_jp2 = zeros(u_nbr^4,1); end - dVarinov(id_e3_xf_u , id_e3_xf_u , jp) = dE_xfxf_E_uu(:,:,jp); - dVarinov(id_e4_u_xf , id_e3_xf_u , jp) = K_ux*dE_xfxf_E_uu(:,:,jp); - dVarinov(id_e3_xf_u , id_e4_u_xf , jp) = transpose(dVarinov(id_e4_u_xf,id_e3_xf_u,jp)); - dVarinov(id_e4_u_xf , id_e4_u_xf , jp) = dVarinov(id_e4_u_xf , id_e3_xf_u , jp)*transpose(K_ux); + dhx_jp2 = dhx(:,:,jp2); + dhu_jp2 = dhu(:,:,jp2); + dhxx_jp2 = dhxx(:,:,jp2); + dhxu_jp2 = dhxu(:,:,jp2); + dhuu_jp2 = dhuu(:,:,jp2); + dhss_jp2 = dhss(:,jp2); + dgx_jp2 = dgx(:,:,jp2); + dgu_jp2 = dgu(:,:,jp2); + dgxx_jp2 = dgxx(:,:,jp2); + dgxu_jp2 = dgxu(:,:,jp2); + dguu_jp2 = dguu(:,:,jp2); + dgss_jp2 = dgss(:,jp2); + dhx_hx_jp2 = kron(dhx_jp2,hx) + kron(hx,dhx_jp2); + dhu_hu_jp2 = kron(dhu_jp2,hu) + kron(hu,dhu_jp2); + dhx_hu_jp2 = kron(dhx_jp2,hu) + kron(hx,dhu_jp2); + dE_xfxf_jp2 = dE_xfxf(:,:,jp2); + dE_xfxf_uu_jp2 = kron(dE_xfxf_jp2,E_uu) + kron(E_xfxf,dE_uu_jp2); + + dA(id_z1_xf , id_z1_xf , jp2) = dhx_jp2; + dA(id_z2_xs , id_z2_xs , jp2) = dhx_jp2; + dA(id_z2_xs , id_z3_xf_xf , jp2) = 1/2*dhxx_jp2; + dA(id_z3_xf_xf , id_z3_xf_xf , jp2) = dhx_hx_jp2; + + dB(id_z1_xf , id_inov1_u , jp2) = dhu_jp2; + dB(id_z2_xs , id_inov2_u_u , jp2) = 1/2*dhuu_jp2; + dB(id_z2_xs , id_inov3_xf_u , jp2) = dhxu_jp2; + dB(id_z3_xf_xf , id_inov2_u_u , jp2) = dhu_hu_jp2; + dB(id_z3_xf_xf , id_inov3_xf_u , jp2) = (I_xx+K_x_x)*dhx_hu_jp2; + + dc(id_z2_xs , jp2) = 1/2*dhss_jp2 + 1/2*dhuu_jp2*E_uu(:) + 1/2*huu*dE_uu_jp2(:); + dc(id_z3_xf_xf , jp2) = dhu_hu_jp2*E_uu(:) + hu_hu*dE_uu_jp2(:); + + dC(: , id_z1_xf , jp2) = dgx_jp2; + dC(: , id_z2_xs , jp2) = dgx_jp2; + dC(: , id_z3_xf_xf , jp2) = 1/2*dgxx_jp2; + + dD(: , id_inov1_u , jp2) = dgu_jp2; + dD(: , id_inov2_u_u , jp2) = 1/2*dguu_jp2; + dD(: , id_inov3_xf_u , jp2) = dgxu_jp2; + + dd(:,jp2) = 1/2*dgss_jp2 + 1/2*guu*dE_uu_jp2(:) + 1/2*dguu_jp2*E_uu(:); + + dVarinov(id_inov1_u , id_inov1_u , jp2) = dE_uu_jp2; + dVarinov(id_inov2_u_u , id_inov2_u_u , jp2) = reshape(dE_u_u_u_u_jp2,u_nbr^2,u_nbr^2) - dE_uu_jp2(:)*E_uu(:)' - E_uu(:)*dE_uu_jp2(:)'; + dVarinov(id_inov3_xf_u , id_inov3_xf_u , jp2) = dE_xfxf_uu_jp2; - dE_xs(:,jp) = I_hxinv*( dghx(indx,:,jp)*E_xs + 1/2*(dghxx(indx,:,jp)*E_xfxf(:) + ghxx(indx,:)*vec(dE_xfxf(:,:,jp)) + dc(x_nbr+(1:x_nbr),jp)) ); - dE_z(:,jp) = [zeros(x_nbr,1); dE_xs(:,jp); vec(dE_xfxf(:,:,jp))]; + dE_xs(:,jp2) = invIx_hx*( dhx_jp2*E_xs + 1/2*dhxx_jp2*E_xfxf(:) + 1/2*hxx*dE_xfxf_jp2(:) + dc(id_z2_xs,jp2) ); + dOm_z_jp2 = dB(:,:,jp2)*Varinov*B' + B*dVarinov(:,:,jp2)*B' + B*Varinov*dB(:,:,jp2)'; + + [dVar_z(:,:,jp2), errorflag] = disclyap_fast(A, dA(:,:,jp2)*Var_z*A' + A*Var_z*dA(:,:,jp2)' + dOm_z_jp2, options.lyapunov_doubling_tol); + if errorflag + dVar_z(:,:,jp2) = lyapunov_symm(A, dA(:,:,jp2)*Var_z*A' + A*Var_z*dA(:,:,jp2)' + dOm_z_jp2,... + options.lyapunov_fixed_point_tol,options.qz_criterium,options.lyapunov_complex_threshold,... + lyapunov_symm_method,... + options.debug); + if lyapunov_symm_method == 1 + lyapunov_symm_method = 2; %now we can reuse persistent schur + end + end + % Make sure some stuff is zero due to Gaussianity + dVar_z(id_z1_xf , id_z2_xs , jp2) = zeros(x_nbr,x_nbr); + dVar_z(id_z1_xf , id_z3_xf_xf , jp2) = zeros(x_nbr,x_nbr^2); + dVar_z(id_z2_xs , id_z1_xf , jp2) = zeros(x_nbr,x_nbr); + dVar_z(id_z3_xf_xf , id_z1_xf , jp2) = zeros(x_nbr^2,x_nbr); end end if order > 2 - Q6Pu = Q6_plication(u_nbr); %quadruplication matrix as in Meijer (2005), similar to Magnus-Neudecker definition of duplication matrix (i.e. dyn_vec) but for fourth-order product moments - %Compute unique product moments of E[kron(u*u',uu')] = E[kron(u,u)*kron(u,u)'] - COMBOS6 = flipud(allVL1(u_nbr, 6)); %all possible combinations of powers that sum up to four for fourth-order product moments of u - E_u_u_u_u_u_u = zeros(u_nbr*(u_nbr+1)/2*(u_nbr+2)/3*(u_nbr+3)/4*(u_nbr+4)/5*(u_nbr+5)/6,1); %only unique entries of E[kron(u,kron(u,kron(u,kron(u,kron(u,u)))))] - if compute_derivs && (stderrparam_nbr+corrparam_nbr>0) - dE_u_u_u_u_u_u = zeros(size(E_u_u_u_u_u_u,1),stderrparam_nbr+corrparam_nbr); - end - for j = 1:size(COMBOS6,1) - E_u_u_u_u_u_u(j) = prodmom(E_uu, 1:u_nbr, COMBOS6(j,:)); - if compute_derivs && (stderrparam_nbr+corrparam_nbr>0) - dE_u_u_u_u_u_u(j,1:(stderrparam_nbr+corrparam_nbr)) = prodmom_deriv(E_uu, 1:u_nbr, COMBOS6(j,:), dE_uu(:,:,1:(stderrparam_nbr+corrparam_nbr)), dCorrelation_matrix(:,:,1:(stderrparam_nbr+corrparam_nbr))); - end - end + % Some common and useful objects for order > 2 + K_u_xx = commutation(u_nbr,x_nbr^2,1); + K_u_ux = commutation(u_nbr,u_nbr*x_nbr,1); + hx_hss2 = kron(hx,1/2*hss); + hu_hss2 = kron(hu,1/2*hss); + hx_hxx2 = kron(hx,1/2*hxx); + hxx2_hu = kron(1/2*hxx,hu); + hx_hxu = kron(hx,hxu); + hxu_hu = kron(hxu,hu); + hx_huu2 = kron(hx,1/2*huu); + hu_huu2 = kron(hu,1/2*huu); + hx_hx_hx = kron(hx,hx_hx); + hx_hx_hu = kron(hx_hx,hu); + hu_hx_hx = kron(hu,hx_hx); + hu_hu_hu = kron(hu_hu,hu); + hx_hu_hu = kron(hx,hu_hu); + hu_hx_hu = kron(hu,hx_hu); + invIxx_hx_hx = (eye(x_nbr^2)-hx_hx)\eye(x_nbr^2); - Var_z = lyapunov_symm(A, B*Varinov*transpose(B), options.lyapunov_fixed_point_tol, options.qz_criterium, options.lyapunov_complex_threshold, 1, options.debug); - %E_xfxf = Var_z(id_z1_xf,id_z1_xf); %this is E_xfxf, as E_xf=0, we already have that - E_xfxs = Var_z(id_z1_xf,id_z2_xs); %as E_xf=0 - E_xf_xf_xf = vec(Var_z(id_z1_xf,id_z3_xf_xf)); %as E_xf=0, this is vec(E[xf*kron(xf,xf)']) - E_xsxf = Var_z(id_z2_xs,id_z1_xf); %as E_xf=0 - E_xsxs = Var_z(id_z2_xs,id_z2_xs) + E_xs*transpose(E_xs); - E_xs_xf_xf = vec(Var_z(id_z2_xs,id_z3_xf_xf)+E_xs*E_xfxf(:)'); %this is vec(E[xs*kron(xf,xf)']) - %E_xf_xf_xf = vec(Var_z(id_z3_xf_xf,id_z1_xf)); - E_xf_xf_xs = vec(Var_z(id_z3_xf_xf,id_z2_xs) + E_xfxf(:)*E_xs'); - E_xf_xf_xf_xf = vec(Var_z(id_z3_xf_xf,id_z3_xf_xf) + E_xfxf(:)*E_xfxf(:)'); - %E_xs_E_uu = kron(E_xs,E_uu); - E_xfxs_E_uu = kron(E_xfxs,E_uu); - - z_nbr = x_nbr+x_nbr+x_nbr^2+x_nbr+x_nbr*x_nbr+x_nbr^3; - e_nbr = u_nbr+u_nbr^2+x_nbr*u_nbr+u_nbr*x_nbr+x_nbr*u_nbr+u_nbr*x_nbr+x_nbr*x_nbr*u_nbr+x_nbr*u_nbr*x_nbr+u_nbr*x_nbr*x_nbr+u_nbr*u_nbr*x_nbr+u_nbr*x_nbr*u_nbr+x_nbr*u_nbr*u_nbr+u_nbr*u_nbr*u_nbr; - hx_hx = kron(ghx(indx,:),ghx(indx,:)); - hu_hu = kron(ghu(indx,:),ghu(indx,:)); - hx_hu = kron(ghx(indx,:),ghu(indx,:)); - hu_hx = kron(ghu(indx,:),ghx(indx,:)); - + % Reuse second-order results + %E_xfxf = Var_z(id_z1_xf, id_z1_xf ); %this is E[xf*xf'], we already have that + %E_xfxs = Var_z(id_z1_xf, id_z2_xs ); %this is E[xf*xs']=0 due to gaussianity + %E_xfxf_xf = Var_z(id_z1_xf, id_z3_xf_xf ); %this is E[xf*kron(xf_xf)']=0 due to gaussianity + %E_xsxf = Var_z(id_z2_xs, id_z1_xf ); %this is E[xs*xf']=0 due to gaussianity + E_xsxs = Var_z(id_z2_xs, id_z2_xs ) + E_xs*transpose(E_xs); %this is E[xs*xs'] + E_xsxf_xf = Var_z(id_z2_xs, id_z3_xf_xf ) + E_xs*E_xfxf(:)'; %this is E[xs*kron(xf,xf)'] + %E_xf_xfxf = Var_z(id_z3_xf_xf, id_z1_xf ); %this is E[kron(xf,xf)*xf']=0 due to gaussianity + E_xf_xfxs = Var_z(id_z3_xf_xf, id_z2_xs ) + E_xfxf(:)*E_xs'; %this is E[kron(xf,xf)*xs'] + E_xf_xfxf_xf = Var_z(id_z3_xf_xf, id_z3_xf_xf) + E_xfxf(:)*E_xfxf(:)'; %this is E[kron(xf,xf)*kron(xf,xf)'] + E_xrdxf = reshape(invIxx_hx_hx*vec(... + hxx*reshape( commutation(x_nbr^2,x_nbr,1)*E_xf_xfxs(:), x_nbr^2,x_nbr)*hx'... + + hxu*kron(E_xs,E_uu)*hu'... + + 1/6*hxxx*reshape(E_xf_xfxf_xf,x_nbr^3,x_nbr)*hx'... + + 1/6*huuu*reshape(QPu*E_u_u_u_u,u_nbr^3,u_nbr)*hu'... + + 3/6*hxxu*kron(E_xfxf(:),E_uu)*hu'... + + 3/6*hxuu*kron(E_xfxf,E_uu(:))*hx'... + + 3/6*hxss*E_xfxf*hx'... + + 3/6*huss*E_uu*hu'... + ),... + x_nbr,x_nbr); %this is E[xrd*xf'] if compute_derivs - dVar_z = zeros(size(Var_z,1),size(Var_z,2),totparam_nbr); - dE_xfxs = zeros(size(E_xfxs,1),size(E_xfxs,2),totparam_nbr); - dE_xf_xf_xf = zeros(size(E_xf_xf_xf,1),totparam_nbr); - dE_xsxf = zeros(size(E_xsxf,1),size(E_xsxf,2),totparam_nbr); - dE_xsxs = zeros(size(E_xsxs,1),size(E_xsxs,2),totparam_nbr); - dE_xs_xf_xf = zeros(size(E_xs_xf_xf,1),totparam_nbr); - dE_xf_xf_xs = zeros(size(E_xf_xf_xs,1),totparam_nbr); - dE_xf_xf_xf_xf = zeros(size(E_xf_xf_xf_xf,1),totparam_nbr); - dE_xfxs_E_uu = zeros(size(E_xfxs_E_uu,1),size(E_xfxs_E_uu,2),totparam_nbr); - dhx_hx = zeros(size(hx_hx,1),size(hx_hx,2),totparam_nbr); - dhu_hu = zeros(size(hu_hu,1),size(hu_hu,2),totparam_nbr); - dhx_hu = zeros(size(hx_hu,1),size(hx_hu,2),totparam_nbr); - dhu_hx = zeros(size(hu_hx,1),size(hu_hx,2),totparam_nbr); - for jp=1:totparam_nbr - dVar_z(:,:,jp) = lyapunov_symm(A, dB(:,:,jp)*Varinov*transpose(B) + B*dVarinov(:,:,jp)*transpose(B) +B*Varinov*transpose(dB(:,:,jp)) + dA(:,:,jp)*Var_z*A' + A*Var_z*transpose(dA(:,:,jp)),... ,... - options.lyapunov_fixed_point_tol, options.qz_criterium, options.lyapunov_complex_threshold, 2, options.debug); %2 is used as we use 1 above - dE_xfxs(:,:,jp) = dVar_z(id_z1_xf,id_z2_xs,jp); - dE_xf_xf_xf(:,jp) = vec(dVar_z(id_z1_xf,id_z3_xf_xf,jp)); - dE_xsxf(:,:,jp) = dVar_z(id_z2_xs,id_z1_xf,jp); - dE_xsxs(:,:,jp) = dVar_z(id_z2_xs,id_z2_xs,jp) + dE_xs(:,jp)*transpose(E_xs) + E_xs*transpose(dE_xs(:,jp)); - dE_xs_xf_xf(:,jp) = vec(dVar_z(id_z2_xs,id_z3_xf_xf,jp) + dE_xs(:,jp)*E_xfxf(:)' + E_xs*vec(dE_xfxf(:,:,jp))'); - dE_xf_xf_xs(:,jp) = vec(dVar_z(id_z3_xf_xf,id_z2_xs,jp) + vec(dE_xfxf(:,:,jp))*E_xs' + E_xfxf(:)*dE_xs(:,jp)'); - dE_xf_xf_xf_xf(:,jp) = vec(dVar_z(id_z3_xf_xf,id_z3_xf_xf,jp) + vec(dE_xfxf(:,:,jp))*E_xfxf(:)' + E_xfxf(:)*vec(dE_xfxf(:,:,jp))'); - dE_xfxs_E_uu(:,:,jp) = kron(dE_xfxs(:,:,jp),E_uu) + kron(E_xfxs,dE_uu(:,:,jp)); - - dhx_hx(:,:,jp) = kron(dghx(indx,:,jp),ghx(indx,:)) + kron(ghx(indx,:),dghx(indx,:,jp)); - dhu_hu(:,:,jp) = kron(dghu(indx,:,jp),ghu(indx,:)) + kron(ghu(indx,:),dghu(indx,:,jp)); - dhx_hu(:,:,jp) = kron(dghx(indx,:,jp),ghu(indx,:)) + kron(ghx(indx,:),dghu(indx,:,jp)); - dhu_hx(:,:,jp) = kron(dghu(indx,:,jp),ghx(indx,:)) + kron(ghu(indx,:),dghx(indx,:,jp)); + dE_xsxs = zeros(x_nbr,x_nbr,totparam_nbr); + dE_xsxf_xf = zeros(x_nbr,x_nbr^2,totparam_nbr); + dE_xf_xfxs = zeros(x_nbr^2,x_nbr,totparam_nbr); + dE_xf_xfxf_xf = zeros(x_nbr^2,x_nbr^2,totparam_nbr); + dE_xrdxf = zeros(x_nbr,x_nbr,totparam_nbr); + for jp2 = 1:totparam_nbr + if jp2 < (stderrparam_nbr+corrparam_nbr) + dE_u_u_u_u_jp2 = QPu*dE_u_u_u_u(:,jp2); + else + dE_u_u_u_u_jp2 = zeros(u_nbr^4,1); + end + dE_xsxs(:,:,jp2) = dVar_z(id_z2_xs , id_z2_xs , jp2) + dE_xs(:,jp2)*transpose(E_xs) + E_xs*transpose(dE_xs(:,jp2)); + dE_xsxf_xf(:,:,jp2) = dVar_z(id_z2_xs , id_z3_xf_xf , jp2) + dE_xs(:,jp2)*E_xfxf(:)' + E_xs*vec(dE_xfxf(:,:,jp2))'; + dE_xf_xfxs(:,:,jp2) = dVar_z(id_z3_xf_xf , id_z2_xs , jp2) + vec(dE_xfxf(:,:,jp2))*E_xs' + E_xfxf(:)*dE_xs(:,jp2)'; + dE_xf_xfxf_xf(:,:,jp2) = dVar_z(id_z3_xf_xf , id_z3_xf_xf , jp2) + vec(dE_xfxf(:,:,jp2))*E_xfxf(:)' + E_xfxf(:)*vec(dE_xfxf(:,:,jp2))'; + dE_xrdxf(:,:,jp2) = reshape(invIxx_hx_hx*vec(... + dhx(:,:,jp2)*E_xrdxf*hx' + hx*E_xrdxf*dhx(:,:,jp2)'... + + dhxx(:,:,jp2)*reshape( commutation(x_nbr^2,x_nbr,1)*E_xf_xfxs(:), x_nbr^2,x_nbr)*hx' + hxx*reshape( commutation(x_nbr^2,x_nbr,1)*vec(dE_xf_xfxs(:,:,jp2)), x_nbr^2,x_nbr)*hx' + hxx*reshape( commutation(x_nbr^2,x_nbr,1)*E_xf_xfxs(:), x_nbr^2,x_nbr)*dhx(:,:,jp2)'... + + dhxu(:,:,jp2)*kron(E_xs,E_uu)*hu' + hxu*kron(dE_xs(:,jp2),E_uu)*hu' + hxu*kron(E_xs,dE_uu(:,:,jp2))*hu' + hxu*kron(E_xs,E_uu)*dhu(:,:,jp2)'... + + 1/6*dhxxx(:,:,jp2)*reshape(E_xf_xfxf_xf,x_nbr^3,x_nbr)*hx' + 1/6*hxxx*reshape(dE_xf_xfxf_xf(:,:,jp2),x_nbr^3,x_nbr)*hx' + 1/6*hxxx*reshape(E_xf_xfxf_xf,x_nbr^3,x_nbr)*dhx(:,:,jp2)'... + + 1/6*dhuuu(:,:,jp2)*reshape(QPu*E_u_u_u_u,u_nbr^3,u_nbr)*hu' + 1/6*huuu*reshape(dE_u_u_u_u_jp2,u_nbr^3,u_nbr)*hu' + 1/6*huuu*reshape(QPu*E_u_u_u_u,u_nbr^3,u_nbr)*dhu(:,:,jp2)'... + + 3/6*dhxxu(:,:,jp2)*kron(E_xfxf(:),E_uu)*hu' + 3/6*hxxu*kron(vec(dE_xfxf(:,:,jp2)),E_uu)*hu' + 3/6*hxxu*kron(E_xfxf(:),dE_uu(:,:,jp2))*hu' + 3/6*hxxu*kron(E_xfxf(:),E_uu)*dhu(:,:,jp2)'... + + 3/6*dhxuu(:,:,jp2)*kron(E_xfxf,E_uu(:))*hx' + 3/6*hxuu*kron(dE_xfxf(:,:,jp2),E_uu(:))*hx' + 3/6*hxuu*kron(E_xfxf,vec(dE_uu(:,:,jp2)))*hx' + 3/6*hxuu*kron(E_xfxf,E_uu(:))*dhx(:,:,jp2)'... + + 3/6*dhxss(:,:,jp2)*E_xfxf*hx' + 3/6*hxss*dE_xfxf(:,:,jp2)*hx' + 3/6*hxss*E_xfxf*dhx(:,:,jp2)'... + + 3/6*dhuss(:,:,jp2)*E_uu*hu' + 3/6*huss*dE_uu(:,:,jp2)*hu' + 3/6*huss*E_uu*dhu(:,:,jp2)'... + ), x_nbr, x_nbr); end end + % Compute unique sixth-order product moments of u, i.e. unique(E[kron(kron(kron(kron(kron(u,u),u),u),u),u)],'stable') + u_nbr6 = u_nbr*(u_nbr+1)/2*(u_nbr+2)/3*(u_nbr+3)/4*(u_nbr+4)/5*(u_nbr+5)/6; + Q6Pu = Q6_plication(u_nbr); + COMBOS6 = flipud(allVL1(u_nbr, 6)); %all possible (unique) combinations of powers that sum up to six + E_u_u_u_u_u_u = zeros(u_nbr6,1); %only unique entries + if compute_derivs && (stderrparam_nbr+corrparam_nbr>0) + dE_u_u_u_u_u_u = zeros(u_nbr6,stderrparam_nbr+corrparam_nbr); + end + for j6 = 1:size(COMBOS6,1) + if compute_derivs && (stderrparam_nbr+corrparam_nbr>0) + [E_u_u_u_u_u_u(j6), dE_u_u_u_u_u_u(j6,:)] = prodmom_deriv(E_uu, 1:u_nbr, COMBOS6(j6,:), dE_uu(:,:,1:(stderrparam_nbr+corrparam_nbr)), dr.derivs.dCorrelation_matrix(:,:,1:(stderrparam_nbr+corrparam_nbr))); + else + E_u_u_u_u_u_u(j6) = prodmom(E_uu, 1:u_nbr, COMBOS6(j6,:)); + end + end + +%% Third-order pruned state space system +% Auxiliary state vector z is defined by: z = [xf; xs; kron(xf,xf); xrd; kron(xf,xs); kron(kron(xf,xf),xf)] +% Auxiliary innovations vector inov is defined by: inov = [u; kron(u,u)-E_uu(:); kron(xf,u); kron(xs,u); kron(kron(xf,xf),u); kron(kron(xf,u),u); kron(kron(u,u),u))] + z_nbr = x_nbr + x_nbr + x_nbr^2 + x_nbr + x_nbr^2 + x_nbr^3; + inov_nbr = u_nbr + u_nbr^2 + x_nbr*u_nbr + x_nbr*u_nbr + x_nbr^2*u_nbr + x_nbr*u_nbr^2 + u_nbr^3; + A = zeros(z_nbr,z_nbr); - A(id_z1_xf , id_z1_xf ) = ghx(indx,:); - A(id_z2_xs , id_z2_xs ) = ghx(indx,:); - A(id_z2_xs , id_z3_xf_xf ) = 1/2*ghxx(indx,:); + A(id_z1_xf , id_z1_xf ) = hx; + A(id_z2_xs , id_z2_xs ) = hx; + A(id_z2_xs , id_z3_xf_xf ) = 1/2*hxx; A(id_z3_xf_xf , id_z3_xf_xf ) = hx_hx; - A(id_z4_xrd , id_z1_xf ) = 3/6*ghxss(indx,:); - A(id_z4_xrd , id_z4_xrd ) = ghx(indx,:); - A(id_z4_xrd , id_z5_xf_xs ) = ghxx(indx,:); - A(id_z4_xrd , id_z6_xf_xf_xf) = 1/6*ghxxx(indx,:); - A(id_z5_xf_xs , id_z1_xf ) = kron(ghx(indx,:),1/2*ghs2(indx,:)); + A(id_z4_xrd , id_z1_xf ) = 3/6*hxss; + A(id_z4_xrd , id_z4_xrd ) = hx; + A(id_z4_xrd , id_z5_xf_xs ) = hxx; + A(id_z4_xrd , id_z6_xf_xf_xf) = 1/6*hxxx; + A(id_z5_xf_xs , id_z1_xf ) = hx_hss2; A(id_z5_xf_xs , id_z5_xf_xs ) = hx_hx; - A(id_z5_xf_xs , id_z6_xf_xf_xf) = kron(ghx(indx,:),1/2*ghxx(indx,:)); - A(id_z6_xf_xf_xf , id_z6_xf_xf_xf) = kron(ghx(indx,:),hx_hx); + A(id_z5_xf_xs , id_z6_xf_xf_xf) = hx_hxx2; + A(id_z6_xf_xf_xf , id_z6_xf_xf_xf) = hx_hx_hx; - B = zeros(z_nbr,e_nbr); - B(id_z1_xf , id_e1_u ) = ghu(indx,:); - B(id_z2_xs , id_e2_u_u ) = 1/2*ghuu(indx,:); - B(id_z2_xs , id_e3_xf_u ) = ghxu(indx,:); - B(id_z3_xf_xf , id_e2_u_u ) = hu_hu; - B(id_z3_xf_xf , id_e3_xf_u ) = hx_hu; - B(id_z3_xf_xf , id_e4_u_xf ) = hu_hx; - B(id_z4_xrd , id_e1_u ) = 3/6*ghuss(indx,:); - B(id_z4_xrd , id_e5_xs_u ) = ghxu(indx,:); - B(id_z4_xrd , id_e7_xf_xf_u) = 3/6*ghxxu(indx,:); - B(id_z4_xrd , id_e10_xf_u_u) = 3/6*ghxuu(indx,:); - B(id_z4_xrd , id_e13_u_u_u ) = 1/6*ghuuu(indx,:); - B(id_z5_xf_xs , id_e1_u ) = kron(ghu(indx,:),1/2*ghs2(indx,:)); - B(id_z5_xf_xs , id_e6_u_xs ) = hu_hx; - B(id_z5_xf_xs , id_e7_xf_xf_u) = kron(ghx(indx,:),ghxu(indx,:)); - B(id_z5_xf_xs , id_e9_u_xf_xf) = kron(ghu(indx,:),1/2*ghxx(indx,:)); - B(id_z5_xf_xs , id_e10_xf_u_u) = kron(ghx(indx,:),1/2*ghuu(indx,:)); - B(id_z5_xf_xs , id_e11_u_xf_u) = kron(ghu(indx,:),ghxu(indx,:)); - B(id_z5_xf_xs , id_e13_u_u_u ) = kron(ghu(indx,:),1/2*ghuu(indx,:)); - B(id_z6_xf_xf_xf , id_e7_xf_xf_u) = kron(hx_hx,ghu(indx,:)); - B(id_z6_xf_xf_xf , id_e8_xf_u_xf) = kron(ghx(indx,:),hu_hx); - B(id_z6_xf_xf_xf , id_e9_u_xf_xf) = kron(ghu(indx,:),hx_hx); - B(id_z6_xf_xf_xf , id_e10_xf_u_u) = kron(hx_hu,ghu(indx,:)); - B(id_z6_xf_xf_xf , id_e11_u_xf_u) = kron(ghu(indx,:),hx_hu); - B(id_z6_xf_xf_xf , id_e12_u_u_xf) = kron(hu_hu,ghx(indx,:)); - B(id_z6_xf_xf_xf , id_e13_u_u_u ) = kron(ghu(indx,:),hu_hu); + B = zeros(z_nbr,inov_nbr); + B(id_z1_xf , id_inov1_u ) = hu; + B(id_z2_xs , id_inov2_u_u ) = 1/2*huu; + B(id_z2_xs , id_inov3_xf_u ) = hxu; + B(id_z3_xf_xf , id_inov2_u_u ) = hu_hu; + B(id_z3_xf_xf , id_inov3_xf_u ) = (I_xx+K_x_x)*hx_hu; + B(id_z4_xrd , id_inov1_u ) = 3/6*huss; + B(id_z4_xrd , id_inov4_xs_u ) = hxu; + B(id_z4_xrd , id_inov5_xf_xf_u) = 3/6*hxxu; + B(id_z4_xrd , id_inov6_xf_u_u ) = 3/6*hxuu; + B(id_z4_xrd , id_inov7_u_u_u ) = 1/6*huuu; + B(id_z5_xf_xs , id_inov1_u ) = hu_hss2; + B(id_z5_xf_xs , id_inov4_xs_u ) = K_x_x*hx_hu; + B(id_z5_xf_xs , id_inov5_xf_xf_u) = hx_hxu + K_x_x*hxx2_hu; + B(id_z5_xf_xs , id_inov6_xf_u_u ) = hx_huu2 + K_x_x*hxu_hu; + B(id_z5_xf_xs , id_inov7_u_u_u ) = hu_huu2; + B(id_z6_xf_xf_xf , id_inov5_xf_xf_u) = hx_hx_hu + kron(hx,K_x_x*hx_hu) + hu_hx_hx*K_u_xx; + B(id_z6_xf_xf_xf , id_inov6_xf_u_u ) = hx_hu_hu + hu_hx_hu*K_u_ux + kron(hu,K_x_x*hx_hu)*K_u_ux; + B(id_z6_xf_xf_xf , id_inov7_u_u_u ) = hu_hu_hu; - C = zeros(y_nbr,z_nbr); - C(1:y_nbr , id_z1_xf ) = ghx(indy,:) + 1/2*ghxss(indy,:); - C(1:y_nbr , id_z2_xs ) = ghx(indy,:); - C(1:y_nbr , id_z3_xf_xf ) = 1/2*ghxx(indy,:); - C(1:y_nbr , id_z4_xrd ) = ghx(indy,:); - C(1:y_nbr , id_z5_xf_xs ) = ghxx(indy,:); - C(1:y_nbr , id_z6_xf_xf_xf) = 1/6*ghxxx(indy,:); - - D = zeros(y_nbr,e_nbr); - D(1:y_nbr , id_e1_u ) = ghu(indy,:) + 1/2*ghuss(indy,:); - D(1:y_nbr , id_e2_u_u ) = 1/2*ghuu(indy,:); - D(1:y_nbr , id_e3_xf_u ) = ghxu(indy,:); - D(1:y_nbr , id_e5_xs_u ) = ghxu(indy,:); - D(1:y_nbr , id_e7_xf_xf_u) = 1/2*ghxxu(indy,:); - D(1:y_nbr , id_e10_xf_u_u) = 1/2*ghxuu(indy,:); - D(1:y_nbr , id_e13_u_u_u ) = 1/6*ghuuu(indy,:); - - c = zeros(z_nbr,1); - c(id_z2_xs , 1) = 1/2*ghs2(indx,:) + 1/2*ghuu(indx,:)*E_uu(:); + c = zeros(z_nbr, 1); + c(id_z2_xs , 1) = 1/2*hss + 1/2*huu*E_uu(:); c(id_z3_xf_xf , 1) = hu_hu*E_uu(:); - d = zeros(y_nbr,1); - d(1:y_nbr , 1) = 0.5*ghs2(indy,:) + 0.5*ghuu(indy,:)*E_uu(:); + C = zeros(y_nbr,z_nbr); + C(: , id_z1_xf ) = gx + 3/6*gxss; + C(: , id_z2_xs ) = gx; + C(: , id_z3_xf_xf ) = 1/2*gxx; + C(: , id_z4_xrd ) = gx; + C(: , id_z5_xf_xs ) = gxx; + C(: , id_z6_xf_xf_xf) = 1/6*gxxx; - Varinov = zeros(e_nbr,e_nbr); - Varinov(id_e1_u , id_e1_u ) = E_uu; - Varinov(id_e1_u , id_e5_xs_u ) = kron(E_xs',E_uu); - Varinov(id_e1_u , id_e6_u_xs ) = kron(E_uu,E_xs'); - Varinov(id_e1_u , id_e7_xf_xf_u) = kron(E_xfxf(:)',E_uu); - Varinov(id_e1_u , id_e8_xf_u_xf) = kron(E_uu,E_xfxf(:)')*kron(K_xu,speye(x_nbr)); - Varinov(id_e1_u , id_e9_u_xf_xf) = kron(E_uu,E_xfxf(:)'); - Varinov(id_e1_u , id_e13_u_u_u ) = reshape(E_u_u_u_u,u_nbr,u_nbr^3); - - Varinov(id_e2_u_u , id_e2_u_u ) = E_uu_uu-E_uu(:)*transpose(E_uu(:)); - - Varinov(id_e3_xf_u , id_e3_xf_u ) = E_xfxf_E_uu; - Varinov(id_e3_xf_u , id_e4_u_xf ) = E_xfxf_E_uu*K_ux'; - Varinov(id_e3_xf_u , id_e5_xs_u ) = E_xfxs_E_uu; - Varinov(id_e3_xf_u , id_e6_u_xs ) = E_xfxs_E_uu*K_ux'; - Varinov(id_e3_xf_u , id_e7_xf_xf_u) = kron(reshape(E_xf_xf_xf,x_nbr,x_nbr^2), E_uu); - Varinov(id_e3_xf_u , id_e8_xf_u_xf) = kron(reshape(E_xf_xf_xf,x_nbr,x_nbr^2), E_uu)*kron(speye(x_nbr),K_ux)'; - Varinov(id_e3_xf_u , id_e9_u_xf_xf) = K_xu*kron(E_uu, reshape(E_xf_xf_xf, x_nbr, x_nbr^2)); - - Varinov(id_e4_u_xf , id_e3_xf_u ) = K_ux*E_xfxf_E_uu; - Varinov(id_e4_u_xf , id_e4_u_xf ) = kron(E_uu,E_xfxf); - Varinov(id_e4_u_xf , id_e5_xs_u ) = K_ux*kron(E_xfxs,E_uu); - Varinov(id_e4_u_xf , id_e6_u_xs ) = kron(E_uu, E_xfxs); - Varinov(id_e4_u_xf , id_e7_xf_xf_u) = K_ux*kron(reshape(E_xf_xf_xf,x_nbr,x_nbr^2),E_uu); - Varinov(id_e4_u_xf , id_e8_xf_u_xf) = kron(E_uu,reshape(E_xf_xf_xf,x_nbr,x_nbr^2))*kron(K_xu,speye(x_nbr))'; - Varinov(id_e4_u_xf , id_e9_u_xf_xf) = kron(E_uu,reshape(E_xf_xf_xf,x_nbr,x_nbr^2)); - - Varinov(id_e5_xs_u , id_e1_u ) = kron(E_xs, E_uu); - Varinov(id_e5_xs_u , id_e3_xf_u ) = kron(E_xsxf, E_uu); - Varinov(id_e5_xs_u , id_e4_u_xf ) = kron(E_xsxf, E_uu)*K_ux'; - Varinov(id_e5_xs_u , id_e5_xs_u ) = kron(E_xsxs, E_uu); - Varinov(id_e5_xs_u , id_e6_u_xs ) = kron(E_xsxs, E_uu)*K_ux'; - Varinov(id_e5_xs_u , id_e7_xf_xf_u) = kron(reshape(E_xs_xf_xf,x_nbr,x_nbr^2),E_uu); - Varinov(id_e5_xs_u , id_e8_xf_u_xf) = kron(reshape(E_xs_xf_xf,x_nbr,x_nbr^2),E_uu)*kron(speye(x_nbr),K_ux)'; - Varinov(id_e5_xs_u , id_e9_u_xf_xf) = K_xu*kron(E_uu,reshape(E_xs_xf_xf,x_nbr,x_nbr^2)); - Varinov(id_e5_xs_u , id_e13_u_u_u ) = kron(E_xs,reshape(E_u_u_u_u,u_nbr,u_nbr^3)); - - Varinov(id_e6_u_xs , id_e1_u ) = kron(E_uu,E_xs); - Varinov(id_e6_u_xs , id_e3_xf_u ) = K_ux*kron(E_xsxf, E_uu); - Varinov(id_e6_u_xs , id_e4_u_xf ) = kron(E_uu, E_xsxf); - Varinov(id_e6_u_xs , id_e5_xs_u ) = K_ux*kron(E_xsxs,E_uu); - Varinov(id_e6_u_xs , id_e6_u_xs ) = kron(E_uu, E_xsxs); - Varinov(id_e6_u_xs , id_e7_xf_xf_u) = K_ux*kron(reshape(E_xs_xf_xf,x_nbr,x_nbr^2), E_uu); - Varinov(id_e6_u_xs , id_e8_xf_u_xf) = kron(E_uu, reshape(E_xs_xf_xf,x_nbr,x_nbr^2))*kron(K_xu,speye(x_nbr))'; - Varinov(id_e6_u_xs , id_e9_u_xf_xf) = kron(E_uu, reshape(E_xs_xf_xf,x_nbr,x_nbr^2)); - Varinov(id_e6_u_xs , id_e13_u_u_u ) = K_ux*kron(E_xs,reshape(E_u_u_u_u,u_nbr,u_nbr^3)); - - Varinov(id_e7_xf_xf_u , id_e1_u ) = kron(E_xfxf(:),E_uu); - Varinov(id_e7_xf_xf_u , id_e3_xf_u ) = kron(reshape(E_xf_xf_xf,x_nbr^2,x_nbr),E_uu); - Varinov(id_e7_xf_xf_u , id_e4_u_xf ) = kron(reshape(E_xf_xf_xf,x_nbr^2,x_nbr),E_uu)*K_ux'; - Varinov(id_e7_xf_xf_u , id_e5_xs_u ) = kron(reshape(E_xf_xf_xs,x_nbr^2,x_nbr),E_uu); - Varinov(id_e7_xf_xf_u , id_e6_u_xs ) = kron(reshape(E_xf_xf_xs,x_nbr^2,x_nbr),E_uu)*K_ux'; - Varinov(id_e7_xf_xf_u , id_e7_xf_xf_u) = kron(reshape(E_xf_xf_xf_xf,x_nbr^2,x_nbr^2),E_uu); - Varinov(id_e7_xf_xf_u , id_e8_xf_u_xf) = kron(reshape(E_xf_xf_xf_xf,x_nbr^2,x_nbr^2),E_uu)*kron(speye(x_nbr),K_ux)'; - Varinov(id_e7_xf_xf_u , id_e9_u_xf_xf) = kron(speye(x_nbr),K_ux)*kron(K_ux,speye(x_nbr))*kron(E_uu, reshape(E_xf_xf_xf_xf,x_nbr^2,x_nbr^2)); - Varinov(id_e7_xf_xf_u , id_e13_u_u_u ) = kron(E_xfxf(:),reshape(E_u_u_u_u,u_nbr,u_nbr^3)); - - Varinov(id_e8_xf_u_xf , id_e1_u ) = kron(K_xu,speye(x_nbr))*kron(E_uu,E_xfxf(:)); - Varinov(id_e8_xf_u_xf , id_e3_xf_u ) = kron(speye(x_nbr),K_xu)*kron(reshape(E_xf_xf_xf,x_nbr^2,x_nbr),E_uu); - Varinov(id_e8_xf_u_xf , id_e4_u_xf ) = kron(K_xu,speye(x_nbr))*kron(E_uu,reshape(E_xf_xf_xf,x_nbr^2,x_nbr)); - Varinov(id_e8_xf_u_xf , id_e5_xs_u ) = kron(speye(x_nbr),K_ux)*kron(reshape(E_xf_xf_xs,x_nbr^2,x_nbr),E_uu); - Varinov(id_e8_xf_u_xf , id_e6_u_xs ) = kron(K_xu,speye(x_nbr))*kron(E_uu,reshape(E_xf_xf_xs,x_nbr^2,x_nbr)); - Varinov(id_e8_xf_u_xf , id_e7_xf_xf_u) = kron(speye(x_nbr),K_ux)*kron(reshape(E_xf_xf_xf_xf,x_nbr^2,x_nbr^2),E_uu); - Varinov(id_e8_xf_u_xf , id_e8_xf_u_xf) = kron(K_xu,speye(x_nbr))*kron(E_uu,reshape(E_xf_xf_xf_xf,x_nbr^2,x_nbr^2))*kron(K_xu,speye(x_nbr))'; - Varinov(id_e8_xf_u_xf , id_e9_u_xf_xf) = kron(K_xu,speye(x_nbr))*kron(E_uu,reshape(E_xf_xf_xf_xf,x_nbr^2,x_nbr^2)); - Varinov(id_e8_xf_u_xf , id_e13_u_u_u ) = kron(K_xu,speye(x_nbr))*kron(reshape(E_u_u_u_u,u_nbr,u_nbr^3),E_xfxf(:)); - - Varinov(id_e9_u_xf_xf , id_e1_u ) = kron(E_uu, E_xfxf(:)); - Varinov(id_e9_u_xf_xf , id_e3_xf_u ) = kron(E_uu, reshape(E_xf_xf_xf,x_nbr^2,x_nbr))*K_xu'; - Varinov(id_e9_u_xf_xf , id_e4_u_xf ) = kron(E_uu, reshape(E_xf_xf_xf,x_nbr^2,x_nbr)); - Varinov(id_e9_u_xf_xf , id_e5_xs_u ) = kron(E_uu, reshape(E_xf_xf_xs,x_nbr^2,x_nbr))*K_xu'; - Varinov(id_e9_u_xf_xf , id_e6_u_xs ) = kron(E_uu, reshape(E_xf_xf_xs,x_nbr^2,x_nbr)); - Varinov(id_e9_u_xf_xf , id_e7_xf_xf_u) = kron(speye(x_nbr),K_ux)*kron(K_ux,speye(x_nbr))*kron(reshape(E_xf_xf_xf_xf,x_nbr^2,x_nbr^2),E_uu); - Varinov(id_e9_u_xf_xf , id_e8_xf_u_xf) = kron(E_uu,reshape(E_xf_xf_xf_xf,x_nbr^2,x_nbr^2))*kron(speye(x_nbr),K_ux)'; - Varinov(id_e9_u_xf_xf , id_e9_u_xf_xf) = kron(E_uu,reshape(E_xf_xf_xf_xf,x_nbr^2,x_nbr^2)); - Varinov(id_e9_u_xf_xf , id_e13_u_u_u ) = kron(reshape(E_u_u_u_u,u_nbr,u_nbr^3),E_xfxf(:)); - - Varinov(id_e10_xf_u_u , id_e10_xf_u_u) = kron(E_xfxf,reshape(E_u_u_u_u,u_nbr^2,u_nbr^2)); - Varinov(id_e10_xf_u_u , id_e11_u_xf_u) = kron(E_xfxf,reshape(E_u_u_u_u,u_nbr^2,u_nbr^2))*kron(K_ux,speye(u_nbr))'; - Varinov(id_e10_xf_u_u , id_e12_u_u_xf) = kron(E_xfxf,reshape(E_u_u_u_u,u_nbr^2,u_nbr^2))*kron(K_ux,speye(u_nbr))'*kron(speye(u_nbr),K_ux)'; - - Varinov(id_e11_u_xf_u , id_e10_xf_u_u) = kron(K_ux,speye(u_nbr))*kron(E_xfxf,reshape(E_u_u_u_u,u_nbr^2,u_nbr^2)); - Varinov(id_e11_u_xf_u , id_e11_u_xf_u) = kron(K_ux,speye(u_nbr))*kron(E_xfxf,reshape(E_u_u_u_u,u_nbr^2,u_nbr^2))*kron(K_ux,speye(u_nbr))'; - Varinov(id_e11_u_xf_u , id_e12_u_u_xf) = kron(speye(u_nbr),K_ux)*kron(reshape(E_u_u_u_u,u_nbr^2,u_nbr^2),E_xfxf); - - Varinov(id_e12_u_u_xf , id_e10_xf_u_u) = kron(K_ux,speye(u_nbr))*kron(speye(u_nbr),K_ux)*kron(E_xfxf,reshape(E_u_u_u_u,u_nbr^2,u_nbr^2)); - Varinov(id_e12_u_u_xf , id_e11_u_xf_u) = kron(reshape(E_u_u_u_u,u_nbr^2,u_nbr^2),E_xfxf)*kron(speye(u_nbr),K_xu)'; - Varinov(id_e12_u_u_xf , id_e12_u_u_xf) = kron(reshape(E_u_u_u_u,u_nbr^2,u_nbr^2),E_xfxf); - - Varinov(id_e13_u_u_u , id_e1_u ) = reshape(E_u_u_u_u,u_nbr^3,u_nbr); - Varinov(id_e13_u_u_u , id_e5_xs_u ) = kron(E_xs', reshape(E_u_u_u_u,u_nbr^3,u_nbr)); - Varinov(id_e13_u_u_u , id_e6_u_xs ) = kron(reshape(E_u_u_u_u,u_nbr^3,u_nbr),E_xs'); - Varinov(id_e13_u_u_u , id_e7_xf_xf_u ) = kron(E_xfxf(:)',reshape(E_u_u_u_u,u_nbr^3,u_nbr)); - Varinov(id_e13_u_u_u , id_e8_xf_u_xf ) = kron(E_xfxf(:)',reshape(E_u_u_u_u,u_nbr^3,u_nbr))*kron(speye(x_nbr),K_ux)'; - Varinov(id_e13_u_u_u , id_e9_u_xf_xf ) = kron(reshape(E_u_u_u_u,u_nbr^3,u_nbr), E_xfxf(:)'); - Varinov(id_e13_u_u_u , id_e13_u_u_u ) = reshape(Q6Pu*E_u_u_u_u_u_u,u_nbr^3,u_nbr^3); - - E_z = (speye(z_nbr)-A)\c; - + D = zeros(y_nbr,inov_nbr); + D(: , id_inov1_u ) = gu + 3/6*guss; + D(: , id_inov2_u_u ) = 1/2*guu; + D(: , id_inov3_xf_u ) = gxu; + D(: , id_inov4_xs_u ) = gxu; + D(: , id_inov5_xf_xf_u) = 3/6*gxxu; + D(: , id_inov6_xf_u_u) = 3/6*gxuu; + D(: , id_inov7_u_u_u ) = 1/6*guuu; + + d = 1/2*gss + 1/2*guu*E_uu(:); + + Varinov = zeros(inov_nbr,inov_nbr); + Varinov(id_inov1_u , id_inov1_u ) = E_uu; + %Varinov(id_inov1_u , id_inov2_u_u ) = zeros(u_nbr,u_nbr^2); + %Varinov(id_inov1_u , id_inov3_xf_u ) = zeros(u_nbr,x_nbr*u_nbr); + Varinov(id_inov1_u , id_inov4_xs_u ) = kron(E_xs',E_uu); + Varinov(id_inov1_u , id_inov5_xf_xf_u) = kron(E_xfxf(:)',E_uu); + %Varinov(id_inov1_u , id_inov6_xf_u_u ) = zeros(u_nbr,x_nbr*u_nbr^2); + Varinov(id_inov1_u , id_inov7_u_u_u ) = reshape(QPu*E_u_u_u_u,u_nbr,u_nbr^3); + + %Varinov(id_inov2_u_u , id_inov1_u ) = zeros(u_nbr^2,u_nbr); + Varinov(id_inov2_u_u , id_inov2_u_u ) = reshape(QPu*E_u_u_u_u,u_nbr^2,u_nbr^2)-E_uu(:)*E_uu(:)'; + %Varinov(id_inov2_u_u , id_inov3_xf_u ) = zeros(u_nbr^2,x_nbr*u_nbr); + %Varinov(id_inov2_u_u , id_inov4_xs_u ) = zeros(u_nbr^2,x_nbr*u_nbr); + %Varinov(id_inov2_u_u , id_inov5_xf_xf_u) = zeros(u_nbr^2,x_nbr^2,u_nbr); + %Varinov(id_inov2_u_u , id_inov6_xf_u_u ) = zeros(u_nbr^2,x_nbr*u_nbr^2); + %Varinov(id_inov2_u_u , id_inov7_u_u_u ) = zeros(u_nbr^2,u_nbr^3); + + %Varinov(id_inov3_xf_u , id_inov1_u ) = zeros(x_nbr*u_nbr,u_nbr); + %Varinov(id_inov3_xf_u , id_inov2_u_u ) = zeros(x_nbr*u_nbr,u_nbr^2); + Varinov(id_inov3_xf_u , id_inov3_xf_u ) = E_xfxf_uu; + %Varinov(id_inov3_xf_u , id_inov4_xs_u ) = zeros(x_nbr*u_nbr,x_nbr*u_nbr); + %Varinov(id_inov3_xf_u , id_inov5_xf_xf_u) = zeros(x_nbr*u_nbr,x_nbr^2*u_nbr); + %Varinov(id_inov3_xf_u , id_inov6_xf_u_u ) = zeros(x_nbr*u_nbr,x_nbr*u_nbr^2); + %Varinov(id_inov3_xf_u , id_inov7_u_u_u ) = zeros(x_nbr*u_nbr,u_nbr^3); + + Varinov(id_inov4_xs_u , id_inov1_u ) = kron(E_xs,E_uu); + %Varinov(id_inov4_xs_u , id_inov2_u_u ) = zeros(x_nbr*u_nbr,u_nbr^2); + %Varinov(id_inov4_xs_u , id_inov3_xf_u ) = zeros(x_nbr*u_nbr,x_nbr*u_nbr); + Varinov(id_inov4_xs_u , id_inov4_xs_u ) = kron(E_xsxs,E_uu); + Varinov(id_inov4_xs_u , id_inov5_xf_xf_u) = kron(E_xsxf_xf, E_uu); + %Varinov(id_inov4_xs_u , id_inov6_xf_u_u ) = zeros(x_nbr*u_nbr,x_nbr*u_nbr^2); + Varinov(id_inov4_xs_u , id_inov7_u_u_u ) = kron(E_xs,reshape(QPu*E_u_u_u_u,u_nbr,u_nbr^3)); + + Varinov(id_inov5_xf_xf_u , id_inov1_u ) = kron(E_xfxf(:),E_uu); + %Varinov(id_inov5_xf_xf_u , id_inov2_u_u ) = zeros(x_nbr^2*u_nbr,u_nbr^2); + %Varinov(id_inov5_xf_xf_u , id_inov3_xf_u ) = zeros(x_nbr^2*u_nbr,x_nbr*u_nbr); + Varinov(id_inov5_xf_xf_u , id_inov4_xs_u ) = kron(E_xf_xfxs,E_uu); + Varinov(id_inov5_xf_xf_u , id_inov5_xf_xf_u) = kron(E_xf_xfxf_xf,E_uu); + %Varinov(id_inov5_xf_xf_u , id_inov6_xf_u_u ) = zeros(x_nbr^2*u_nbr,x_nbr*u_nbr^2); + Varinov(id_inov5_xf_xf_u , id_inov7_u_u_u ) = kron(E_xfxf(:),reshape(QPu*E_u_u_u_u,u_nbr,u_nbr^3)); + + %Varinov(id_inov6_xf_u_u , id_inov1_u ) = zeros(x_nbr*u_nbr^2,u_nbr); + %Varinov(id_inov6_xf_u_u , id_inov2_u_u ) = zeros(x_nbr*u_nbr^2,u_nbr^2); + %Varinov(id_inov6_xf_u_u , id_inov3_xf_u ) = zeros(x_nbr*u_nbr^2,x_nbr*u_nbr); + %Varinov(id_inov6_xf_u_u , id_inov4_xs_u ) = zeros(x_nbr*u_nbr^2,x_nbr*u_nbr); + %Varinov(id_inov6_xf_u_u , id_inov5_xf_xf_u) = zeros(x_nbr*u_nbr^2,x_nbr^2*u_nbr); + Varinov(id_inov6_xf_u_u , id_inov6_xf_u_u ) = kron(E_xfxf,reshape(QPu*E_u_u_u_u,u_nbr^2,u_nbr^2)); + %Varinov(id_inov6_xf_u_u , id_inov7_u_u_u ) = zeros(x_nbr*u_nbr^2,u_nbr^3); + + Varinov(id_inov7_u_u_u , id_inov1_u ) = reshape(QPu*E_u_u_u_u,u_nbr^3,u_nbr); + %Varinov(id_inov7_u_u_u , id_inov2_u_u ) = zeros(u_nbr^3,u_nbr^2); + %Varinov(id_inov7_u_u_u , id_inov3_xf_u ) = zeros(u_nbr^3,x_nbr*u_nbr); + Varinov(id_inov7_u_u_u , id_inov4_xs_u ) = kron(E_xs',reshape(QPu*E_u_u_u_u,u_nbr^3,u_nbr)); + Varinov(id_inov7_u_u_u , id_inov5_xf_xf_u) = kron(transpose(E_xfxf(:)),reshape(QPu*E_u_u_u_u,u_nbr^3,u_nbr)); + %Varinov(id_inov7_u_u_u , id_inov6_xf_u_u ) = zeros(u_nbr^3,x_nbr*u_nbr^2); + Varinov(id_inov7_u_u_u , id_inov7_u_u_u ) = reshape(Q6Pu*E_u_u_u_u_u_u,u_nbr^3,u_nbr^3); + + E_xrd = zeros(x_nbr,1);%due to gaussianity + + E_inovzlag1 = zeros(inov_nbr,z_nbr); % Attention: E[inov*z(-1)'] is not equal to zero for a third-order approximation due to kron(kron(xf(-1),u),u) + E_inovzlag1(id_inov6_xf_u_u , id_z1_xf ) = kron(E_xfxf,E_uu(:)); + E_inovzlag1(id_inov6_xf_u_u , id_z4_xrd ) = kron(E_xrdxf',E_uu(:)); + E_inovzlag1(id_inov6_xf_u_u , id_z5_xf_xs ) = kron(reshape(commutation(x_nbr^2,x_nbr)*vec(E_xsxf_xf),x_nbr,x_nbr^2),vec(E_uu)) ; + E_inovzlag1(id_inov6_xf_u_u , id_z6_xf_xf_xf ) = kron(reshape(E_xf_xfxf_xf,x_nbr,x_nbr^3),E_uu(:)); + + Binovzlag1A= B*E_inovzlag1*transpose(A); + Om_z = B*Varinov*transpose(B) + Binovzlag1A + transpose(Binovzlag1A); + + lyapunov_symm_method = 1; %method=1 to initialize persistent variables + [Var_z, errorflag] = disclyap_fast(A,Om_z,options.lyapunov_doubling_tol); + if errorflag %use Schur-based method + fprintf('PRUNED_STATE_SPACE_SYSTEM: error flag in disclyap_fast at order=3, use lyapunov_symm\n'); + Var_z = lyapunov_symm(A,Om_z,... + options.lyapunov_fixed_point_tol,options.qz_criterium,options.lyapunov_complex_threshold,... + lyapunov_symm_method,... + options.debug); + lyapunov_symm_method = 2; %we can now make use of persistent variables from shur + end + %make sure some stuff is zero due to Gaussianity + Var_z(id_z1_xf , id_z2_xs) = zeros(x_nbr,x_nbr); + Var_z(id_z1_xf , id_z3_xf_xf) = zeros(x_nbr,x_nbr^2); + Var_z(id_z2_xs , id_z1_xf) = zeros(x_nbr,x_nbr); + Var_z(id_z2_xs , id_z4_xrd) = zeros(x_nbr,x_nbr); + Var_z(id_z2_xs , id_z5_xf_xs) = zeros(x_nbr,x_nbr^2); + Var_z(id_z2_xs , id_z6_xf_xf_xf) = zeros(x_nbr,x_nbr^3); + Var_z(id_z3_xf_xf , id_z1_xf) = zeros(x_nbr^2,x_nbr); + Var_z(id_z3_xf_xf , id_z4_xrd) = zeros(x_nbr^2,x_nbr); + Var_z(id_z3_xf_xf , id_z5_xf_xs) = zeros(x_nbr^2,x_nbr^2); + Var_z(id_z3_xf_xf , id_z6_xf_xf_xf) = zeros(x_nbr^2,x_nbr^3); + Var_z(id_z4_xrd , id_z2_xs) = zeros(x_nbr,x_nbr); + Var_z(id_z4_xrd , id_z3_xf_xf) = zeros(x_nbr,x_nbr^2); + Var_z(id_z5_xf_xs , id_z2_xs) = zeros(x_nbr^2,x_nbr); + Var_z(id_z5_xf_xs , id_z3_xf_xf) = zeros(x_nbr^2,x_nbr^2); + Var_z(id_z6_xf_xf_xf , id_z2_xs) = zeros(x_nbr^3,x_nbr); + Var_z(id_z6_xf_xf_xf , id_z3_xf_xf) = zeros(x_nbr^3,x_nbr^2); if compute_derivs - dA = zeros(z_nbr,z_nbr,totparam_nbr); - dB = zeros(z_nbr,e_nbr,totparam_nbr); - dC = zeros(y_nbr,z_nbr,totparam_nbr); - dD = zeros(y_nbr,e_nbr,totparam_nbr); - dc = zeros(z_nbr,totparam_nbr); - dd = zeros(y_nbr,totparam_nbr); - dVarinov = zeros(e_nbr,e_nbr,totparam_nbr); - dE_z =zeros(z_nbr,totparam_nbr); - - for jp=1:totparam_nbr - dA(id_z1_xf , id_z1_xf ,jp) = dghx(indx,:,jp); - dA(id_z2_xs , id_z2_xs ,jp) = dghx(indx,:,jp); - dA(id_z2_xs , id_z3_xf_xf ,jp) = 1/2*dghxx(indx,:,jp); - dA(id_z3_xf_xf , id_z3_xf_xf ,jp) = dhx_hx(:,:,jp); - dA(id_z4_xrd , id_z1_xf ,jp) = 3/6*dghxss(indx,:,jp); - dA(id_z4_xrd , id_z4_xrd ,jp) = dghx(indx,:,jp); - dA(id_z4_xrd , id_z5_xf_xs ,jp) = dghxx(indx,:,jp); - dA(id_z4_xrd , id_z6_xf_xf_xf ,jp) = 1/6*dghxxx(indx,:,jp); - dA(id_z5_xf_xs , id_z1_xf ,jp) = kron(dghx(indx,:,jp),1/2*ghs2(indx,:)) + kron(ghx(indx,:),1/2*dghs2(indx,jp)); - dA(id_z5_xf_xs , id_z5_xf_xs ,jp) = dhx_hx(:,:,jp); - dA(id_z5_xf_xs , id_z6_xf_xf_xf ,jp) = kron(dghx(indx,:,jp),1/2*ghxx(indx,:)) + kron(ghx(indx,:),1/2*dghxx(indx,:,jp)); - dA(id_z6_xf_xf_xf , id_z6_xf_xf_xf ,jp) = kron(dghx(indx,:,jp),hx_hx) + kron(ghx(indx,:),dhx_hx(:,:,jp)); + dA = zeros(z_nbr,z_nbr,totparam_nbr); + dB = zeros(z_nbr,inov_nbr,totparam_nbr); + dc = zeros(z_nbr,totparam_nbr); + dC = zeros(y_nbr,z_nbr,totparam_nbr); + dD = zeros(y_nbr,inov_nbr,totparam_nbr); + dd = zeros(y_nbr,totparam_nbr); + dVarinov = zeros(inov_nbr,inov_nbr,totparam_nbr); + dE_xrd = zeros(x_nbr,totparam_nbr); + dE_inovzlag1 = zeros(inov_nbr,z_nbr,totparam_nbr); + dVar_z = zeros(z_nbr,z_nbr,totparam_nbr); + + for jp3 = 1:totparam_nbr + if jp3 <= (stderrparam_nbr+corrparam_nbr) + dE_uu_jp3 = dE_uu(:,:,jp3); + dE_u_u_u_u_jp3 = QPu*dE_u_u_u_u(:,jp3); + dE_u_u_u_u_u_u_jp3 = Q6Pu*dE_u_u_u_u_u_u(:,jp3); + else + dE_uu_jp3 = zeros(u_nbr,u_nbr); + dE_u_u_u_u_jp3 = zeros(u_nbr^4,1); + dE_u_u_u_u_u_u_jp3 = zeros(u_nbr^6,1); + end + dhx_jp3 = dhx(:,:,jp3); + dhu_jp3 = dhu(:,:,jp3); + dhxx_jp3 = dhxx(:,:,jp3); + dhxu_jp3 = dhxu(:,:,jp3); + dhuu_jp3 = dhuu(:,:,jp3); + dhss_jp3 = dhss(:,jp3); + dhxxx_jp3 = dhxxx(:,:,jp3); + dhxxu_jp3 = dhxxu(:,:,jp3); + dhxuu_jp3 = dhxuu(:,:,jp3); + dhuuu_jp3 = dhuuu(:,:,jp3); + dhxss_jp3 = dhxss(:,:,jp3); + dhuss_jp3 = dhuss(:,:,jp3); + dgx_jp3 = dgx(:,:,jp3); + dgu_jp3 = dgu(:,:,jp3); + dgxx_jp3 = dgxx(:,:,jp3); + dgxu_jp3 = dgxu(:,:,jp3); + dguu_jp3 = dguu(:,:,jp3); + dgss_jp3 = dgss(:,jp3); + dgxxx_jp3 = dgxxx(:,:,jp3); + dgxxu_jp3 = dgxxu(:,:,jp3); + dgxuu_jp3 = dgxuu(:,:,jp3); + dguuu_jp3 = dguuu(:,:,jp3); + dgxss_jp3 = dgxss(:,:,jp3); + dguss_jp3 = dguss(:,:,jp3); - dB(id_z1_xf , id_e1_u , jp) = dghu(indx,:,jp); - dB(id_z2_xs , id_e2_u_u , jp) = 1/2*dghuu(indx,:,jp); - dB(id_z2_xs , id_e3_xf_u , jp) = dghxu(indx,:,jp); - dB(id_z3_xf_xf , id_e2_u_u , jp) = dhu_hu(:,:,jp); - dB(id_z3_xf_xf , id_e3_xf_u , jp) = dhx_hu(:,:,jp); - dB(id_z3_xf_xf , id_e4_u_xf , jp) = dhu_hx(:,:,jp); - dB(id_z4_xrd , id_e1_u , jp) = 3/6*dghuss(indx,:,jp); - dB(id_z4_xrd , id_e5_xs_u , jp) = dghxu(indx,:,jp); - dB(id_z4_xrd , id_e7_xf_xf_u , jp) = 3/6*dghxxu(indx,:,jp); - dB(id_z4_xrd , id_e10_xf_u_u , jp) = 3/6*dghxuu(indx,:,jp); - dB(id_z4_xrd , id_e13_u_u_u , jp) = 1/6*dghuuu(indx,:,jp); - dB(id_z5_xf_xs , id_e1_u , jp) = kron(dghu(indx,:,jp),1/2*ghs2(indx,:)) + kron(ghu(indx,:),1/2*dghs2(indx,jp)); - dB(id_z5_xf_xs , id_e6_u_xs , jp) = dhu_hx(:,:,jp); - dB(id_z5_xf_xs , id_e7_xf_xf_u , jp) = kron(dghx(indx,:,jp),ghxu(indx,:)) + kron(ghx(indx,:),dghxu(indx,:,jp)); - dB(id_z5_xf_xs , id_e9_u_xf_xf , jp) = kron(dghu(indx,:,jp),1/2*ghxx(indx,:)) + kron(ghu(indx,:),1/2*dghxx(indx,:,jp)); - dB(id_z5_xf_xs , id_e10_xf_u_u , jp) = kron(dghx(indx,:,jp),1/2*ghuu(indx,:)) + kron(ghx(indx,:),1/2*dghuu(indx,:,jp)); - dB(id_z5_xf_xs , id_e11_u_xf_u , jp) = kron(dghu(indx,:,jp),ghxu(indx,:)) + kron(ghu(indx,:),dghxu(indx,:,jp)); - dB(id_z5_xf_xs , id_e13_u_u_u , jp) = kron(dghu(indx,:,jp),1/2*ghuu(indx,:)) + kron(ghu(indx,:),1/2*dghuu(indx,:,jp)); - dB(id_z6_xf_xf_xf , id_e7_xf_xf_u , jp) = kron(dhx_hx(:,:,jp),ghu(indx,:)) + kron(hx_hx,dghu(indx,:,jp)); - dB(id_z6_xf_xf_xf , id_e8_xf_u_xf , jp) = kron(dghx(indx,:,jp),hu_hx) + kron(ghx(indx,:),dhu_hx(:,:,jp)); - dB(id_z6_xf_xf_xf , id_e9_u_xf_xf , jp) = kron(dghu(indx,:,jp),hx_hx) + kron(ghu(indx,:),dhx_hx(:,:,jp)); - dB(id_z6_xf_xf_xf , id_e10_xf_u_u , jp) = kron(dhx_hu(:,:,jp),ghu(indx,:)) + kron(hx_hu,dghu(indx,:,jp)); - dB(id_z6_xf_xf_xf , id_e11_u_xf_u , jp) = kron(dghu(indx,:,jp),hx_hu) + kron(ghu(indx,:),dhx_hu(:,:,jp)); - dB(id_z6_xf_xf_xf , id_e12_u_u_xf , jp) = kron(dhu_hu(:,:,jp),ghx(indx,:)) + kron(hu_hu,dghx(indx,:,jp)); - dB(id_z6_xf_xf_xf , id_e13_u_u_u , jp) = kron(dghu(indx,:,jp),hu_hu) + kron(ghu(indx,:),dhu_hu(:,:,jp)); - - dC(1:y_nbr , id_z1_xf , jp) = dghx(indy,:,jp) + 1/2*dghxss(indy,:,jp); - dC(1:y_nbr , id_z2_xs , jp) = dghx(indy,:,jp); - dC(1:y_nbr , id_z3_xf_xf , jp) = 1/2*dghxx(indy,:,jp); - dC(1:y_nbr , id_z4_xrd , jp) = dghx(indy,:,jp); - dC(1:y_nbr , id_z5_xf_xs , jp) = dghxx(indy,:,jp); - dC(1:y_nbr , id_z6_xf_xf_xf , jp) = 1/6*dghxxx(indy,:,jp); + dhx_hx_jp3 = kron(dhx_jp3,hx) + kron(hx,dhx_jp3); + dhx_hu_jp3 = kron(dhx_jp3,hu) + kron(hx,dhu_jp3); + dhu_hu_jp3 = kron(dhu_jp3,hu) + kron(hu,dhu_jp3); + dhx_hss2_jp3 = kron(dhx_jp3,1/2*hss) + kron(hx,1/2*dhss_jp3); + dhu_hss2_jp3 = kron(dhu_jp3,1/2*hss) + kron(hu,1/2*dhss_jp3); + dhx_hxx2_jp3 = kron(dhx_jp3,1/2*hxx) + kron(hx,1/2*dhxx_jp3); + dhxx2_hu_jp3 = kron(1/2*dhxx_jp3,hu) + kron(1/2*hxx,dhu_jp3); + dhx_hxu_jp3 = kron(dhx_jp3,hxu) + kron(hx,dhxu_jp3); + dhxu_hu_jp3 = kron(dhxu_jp3,hu) + kron(hxu,dhu_jp3); + dhx_huu2_jp3 = kron(dhx_jp3,1/2*huu) + kron(hx,1/2*dhuu_jp3); + dhu_huu2_jp3 = kron(dhu_jp3,1/2*huu) + kron(hu,1/2*dhuu_jp3); + dhx_hx_hx_jp3 = kron(dhx_jp3,hx_hx) + kron(hx,dhx_hx_jp3); + dhx_hx_hu_jp3 = kron(dhx_hx_jp3,hu) + kron(hx_hx,dhu_jp3); + dhu_hx_hx_jp3 = kron(dhu_jp3,hx_hx) + kron(hu,dhx_hx_jp3); + dhu_hu_hu_jp3 = kron(dhu_hu_jp3,hu) + kron(hu_hu,dhu_jp3); + dhx_hu_hu_jp3 = kron(dhx_jp3,hu_hu) + kron(hx,dhu_hu_jp3); + dhu_hx_hu_jp3 = kron(dhu_jp3,hx_hu) + kron(hu,dhx_hu_jp3); - dD(1:y_nbr , id_e1_u , jp) = dghu(indy,:,jp) + 1/2*dghuss(indy,:,jp); - dD(1:y_nbr , id_e2_u_u , jp) = 1/2*dghuu(indy,:,jp); - dD(1:y_nbr , id_e3_xf_u , jp) = dghxu(indy,:,jp); - dD(1:y_nbr , id_e5_xs_u , jp) = dghxu(indy,:,jp); - dD(1:y_nbr , id_e7_xf_xf_u , jp) = 1/2*dghxxu(indy,:,jp); - dD(1:y_nbr , id_e10_xf_u_u , jp) = 1/2*dghxuu(indy,:,jp); - dD(1:y_nbr , id_e13_u_u_u , jp) = 1/6*dghuuu(indy,:,jp); - - dc(id_z2_xs , jp) = 1/2*dghs2(indx,jp) + 1/2*dghuu(indx,:,jp)*E_uu(:) + 1/2*ghuu(indx,:)*vec(dE_uu(:,:,jp)); - dc(id_z3_xf_xf , jp) = dhu_hu(:,:,jp)*E_uu(:) + hu_hu*vec(dE_uu(:,:,jp)); - - dd(1:y_nbr , jp) = 0.5*dghs2(indy,jp) + 0.5*dghuu(indy,:,jp)*E_uu(:) + 0.5*ghuu(indy,:)*vec(dE_uu(:,:,jp)); - - dVarinov(id_e1_u , id_e1_u , jp) = dE_uu(:,:,jp); - dVarinov(id_e1_u , id_e5_xs_u , jp) = kron(dE_xs(:,jp)',E_uu) + kron(E_xs',dE_uu(:,:,jp)); - dVarinov(id_e1_u , id_e6_u_xs , jp) = kron(dE_uu(:,:,jp),E_xs') + kron(E_uu,dE_xs(:,jp)'); - dVarinov(id_e1_u , id_e7_xf_xf_u , jp) = kron(vec(dE_xfxf(:,:,jp))',E_uu) + kron(E_xfxf(:)',dE_uu(:,:,jp)); - dVarinov(id_e1_u , id_e8_xf_u_xf , jp) = (kron(dE_uu(:,:,jp),E_xfxf(:)') + kron(E_uu,vec(dE_xfxf(:,:,jp))'))*kron(K_xu,speye(x_nbr)); - dVarinov(id_e1_u , id_e9_u_xf_xf , jp) = kron(dE_uu(:,:,jp),E_xfxf(:)') + kron(E_uu,vec(dE_xfxf(:,:,jp))'); - if jp <= (stderrparam_nbr+corrparam_nbr) - dVarinov(id_e1_u , id_e13_u_u_u , jp) = reshape(QPu*dE_u_u_u_u(:,jp),u_nbr,u_nbr^3); - dVarinov(id_e2_u_u , id_e2_u_u , jp) = reshape(QPu*dE_u_u_u_u(:,jp),u_nbr^2,u_nbr^2) - vec(dE_uu(:,:,jp))*transpose(E_uu(:)) - E_uu(:)*transpose(vec(dE_uu(:,:,jp))); - end + dE_xs_jp3 = dE_xs(:,jp3); + dE_xfxf_jp3 = dE_xfxf(:,:,jp3); + dE_xsxs_jp3 = dE_xsxs(:,:,jp3); + dE_xsxf_xf_jp3 = dE_xsxf_xf(:,:,jp3); + dE_xfxf_uu_jp3 = kron(dE_xfxf_jp3,E_uu) + kron(E_xfxf,dE_uu_jp3); + dE_xf_xfxs_jp3 = dE_xf_xfxs(:,:,jp3); + dE_xf_xfxf_xf_jp3 = dE_xf_xfxf_xf(:,:,jp3); + dE_xrdxf_jp3 = dE_xrdxf(:,:,jp3); - dVarinov(id_e3_xf_u , id_e3_xf_u , jp) = dE_xfxf_E_uu(:,:,jp); - dVarinov(id_e3_xf_u , id_e4_u_xf , jp) = dE_xfxf_E_uu(:,:,jp)*K_ux'; - dVarinov(id_e3_xf_u , id_e5_xs_u , jp) = dE_xfxs_E_uu(:,:,jp); - dVarinov(id_e3_xf_u , id_e6_u_xs , jp) = dE_xfxs_E_uu(:,:,jp)*K_ux'; - dVarinov(id_e3_xf_u , id_e7_xf_xf_u , jp) = kron(reshape(dE_xf_xf_xf(:,jp),x_nbr,x_nbr^2), E_uu) + kron(reshape(E_xf_xf_xf,x_nbr,x_nbr^2), dE_uu(:,:,jp)); - dVarinov(id_e3_xf_u , id_e8_xf_u_xf , jp) = (kron(reshape(dE_xf_xf_xf(:,jp),x_nbr,x_nbr^2), E_uu) + kron(reshape(E_xf_xf_xf,x_nbr,x_nbr^2), dE_uu(:,:,jp)))*kron(speye(x_nbr),K_ux)'; - dVarinov(id_e3_xf_u , id_e9_u_xf_xf , jp) = K_xu*(kron(dE_uu(:,:,jp), reshape(E_xf_xf_xf, x_nbr, x_nbr^2)) + kron(E_uu, reshape(dE_xf_xf_xf(:,jp), x_nbr, x_nbr^2))); + dA(id_z1_xf , id_z1_xf , jp3) = dhx_jp3; + dA(id_z2_xs , id_z2_xs , jp3) = dhx_jp3; + dA(id_z2_xs , id_z3_xf_xf , jp3) = 1/2*dhxx_jp3; + dA(id_z3_xf_xf , id_z3_xf_xf , jp3) = dhx_hx_jp3; + dA(id_z4_xrd , id_z1_xf , jp3) = 3/6*dhxss_jp3; + dA(id_z4_xrd , id_z4_xrd , jp3) = dhx_jp3; + dA(id_z4_xrd , id_z5_xf_xs , jp3) = dhxx_jp3; + dA(id_z4_xrd , id_z6_xf_xf_xf , jp3) = 1/6*dhxxx_jp3; + dA(id_z5_xf_xs , id_z1_xf , jp3) = dhx_hss2_jp3; + dA(id_z5_xf_xs , id_z5_xf_xs , jp3) = dhx_hx_jp3; + dA(id_z5_xf_xs , id_z6_xf_xf_xf , jp3) = dhx_hxx2_jp3; + dA(id_z6_xf_xf_xf , id_z6_xf_xf_xf , jp3) = dhx_hx_hx_jp3; - dVarinov(id_e4_u_xf , id_e3_xf_u , jp) = K_ux*dE_xfxf_E_uu(:,:,jp); - dVarinov(id_e4_u_xf , id_e4_u_xf , jp) = kron(dE_uu(:,:,jp),E_xfxf) + kron(E_uu,dE_xfxf(:,:,jp)); - dVarinov(id_e4_u_xf , id_e5_xs_u , jp) = K_ux*(kron(dE_xfxs(:,:,jp),E_uu) + kron(E_xfxs,dE_uu(:,:,jp))); - dVarinov(id_e4_u_xf , id_e6_u_xs , jp) = kron(dE_uu(:,:,jp), E_xfxs) + kron(E_uu, dE_xfxs(:,:,jp)); - dVarinov(id_e4_u_xf , id_e7_xf_xf_u , jp) = K_ux*(kron(reshape(dE_xf_xf_xf(:,jp),x_nbr,x_nbr^2),E_uu) + kron(reshape(E_xf_xf_xf,x_nbr,x_nbr^2),dE_uu(:,:,jp))); - dVarinov(id_e4_u_xf , id_e8_xf_u_xf , jp) = (kron(dE_uu(:,:,jp),reshape(E_xf_xf_xf,x_nbr,x_nbr^2)) + kron(E_uu,reshape(dE_xf_xf_xf(:,jp),x_nbr,x_nbr^2)))*kron(K_xu,speye(x_nbr))'; - dVarinov(id_e4_u_xf , id_e9_u_xf_xf , jp) = kron(dE_uu(:,:,jp),reshape(E_xf_xf_xf,x_nbr,x_nbr^2)) + kron(E_uu,reshape(dE_xf_xf_xf(:,jp),x_nbr,x_nbr^2)); + dB(id_z1_xf , id_inov1_u , jp3) = dhu_jp3; + dB(id_z2_xs , id_inov2_u_u , jp3) = 1/2*dhuu_jp3; + dB(id_z2_xs , id_inov3_xf_u , jp3) = dhxu_jp3; + dB(id_z3_xf_xf , id_inov2_u_u , jp3) = dhu_hu_jp3; + dB(id_z3_xf_xf , id_inov3_xf_u , jp3) = (I_xx+K_x_x)*dhx_hu_jp3; + dB(id_z4_xrd , id_inov1_u , jp3) = 3/6*dhuss_jp3; + dB(id_z4_xrd , id_inov4_xs_u , jp3) = dhxu_jp3; + dB(id_z4_xrd , id_inov5_xf_xf_u , jp3) = 3/6*dhxxu_jp3; + dB(id_z4_xrd , id_inov6_xf_u_u , jp3) = 3/6*dhxuu_jp3; + dB(id_z4_xrd , id_inov7_u_u_u , jp3) = 1/6*dhuuu_jp3; + dB(id_z5_xf_xs , id_inov1_u , jp3) = dhu_hss2_jp3; + dB(id_z5_xf_xs , id_inov4_xs_u , jp3) = K_x_x*dhx_hu_jp3; + dB(id_z5_xf_xs , id_inov5_xf_xf_u , jp3) = dhx_hxu_jp3 + K_x_x*dhxx2_hu_jp3; + dB(id_z5_xf_xs , id_inov6_xf_u_u , jp3) = dhx_huu2_jp3 + K_x_x*dhxu_hu_jp3; + dB(id_z5_xf_xs , id_inov7_u_u_u , jp3) = dhu_huu2_jp3; + dB(id_z6_xf_xf_xf , id_inov5_xf_xf_u , jp3) = dhx_hx_hu_jp3 + kron(dhx_jp3,K_x_x*hx_hu) + kron(hx,K_x_x*dhx_hu_jp3) + dhu_hx_hx_jp3*K_u_xx; + dB(id_z6_xf_xf_xf , id_inov6_xf_u_u , jp3) = dhx_hu_hu_jp3 + dhu_hx_hu_jp3*K_u_ux + kron(dhu_jp3,K_x_x*hx_hu)*K_u_ux + kron(hu,K_x_x*dhx_hu_jp3)*K_u_ux; + dB(id_z6_xf_xf_xf , id_inov7_u_u_u , jp3) = dhu_hu_hu_jp3; - dVarinov(id_e5_xs_u , id_e1_u , jp) = kron(dE_xs(:,jp), E_uu) + kron(E_xs, dE_uu(:,:,jp)); - dVarinov(id_e5_xs_u , id_e3_xf_u , jp) = kron(dE_xsxf(:,:,jp), E_uu) + kron(E_xsxf, dE_uu(:,:,jp)); - dVarinov(id_e5_xs_u , id_e4_u_xf , jp) = (kron(dE_xsxf(:,:,jp), E_uu) + kron(E_xsxf, dE_uu(:,:,jp)))*K_ux'; - dVarinov(id_e5_xs_u , id_e5_xs_u , jp) = kron(dE_xsxs(:,:,jp), E_uu) + kron(E_xsxs, dE_uu(:,:,jp)); - dVarinov(id_e5_xs_u , id_e6_u_xs , jp) = (kron(dE_xsxs(:,:,jp), E_uu) + kron(E_xsxs, dE_uu(:,:,jp)))*K_ux'; - dVarinov(id_e5_xs_u , id_e7_xf_xf_u , jp) = kron(reshape(dE_xs_xf_xf(:,jp),x_nbr,x_nbr^2),E_uu) + kron(reshape(E_xs_xf_xf,x_nbr,x_nbr^2),dE_uu(:,:,jp)); - dVarinov(id_e5_xs_u , id_e8_xf_u_xf , jp) = (kron(reshape(dE_xs_xf_xf(:,jp),x_nbr,x_nbr^2),E_uu) + kron(reshape(E_xs_xf_xf,x_nbr,x_nbr^2),dE_uu(:,:,jp)))*kron(speye(x_nbr),K_ux)'; - dVarinov(id_e5_xs_u , id_e9_u_xf_xf , jp) = K_xu*(kron(dE_uu(:,:,jp),reshape(E_xs_xf_xf,x_nbr,x_nbr^2)) + kron(E_uu,reshape(dE_xs_xf_xf(:,jp),x_nbr,x_nbr^2))); - if jp <= (stderrparam_nbr+corrparam_nbr) - dVarinov(id_e5_xs_u , id_e13_u_u_u , jp) = kron(dE_xs(:,jp),reshape(E_u_u_u_u,u_nbr,u_nbr^3)) + kron(E_xs,reshape(QPu*dE_u_u_u_u(:,jp),u_nbr,u_nbr^3)); - else - dVarinov(id_e5_xs_u , id_e13_u_u_u , jp) = kron(dE_xs(:,jp),reshape(E_u_u_u_u,u_nbr,u_nbr^3)); - end + dc(id_z2_xs , jp3) = 1/2*dhss_jp3 + 1/2*dhuu_jp3*E_uu(:) + 1/2*huu*dE_uu_jp3(:); + dc(id_z3_xf_xf , jp3) = dhu_hu_jp3*E_uu(:) + hu_hu*dE_uu_jp3(:); - dVarinov(id_e6_u_xs , id_e1_u , jp) = kron(dE_uu(:,:,jp),E_xs) + kron(E_uu,dE_xs(:,jp)); - dVarinov(id_e6_u_xs , id_e3_xf_u , jp) = K_ux*(kron(dE_xsxf(:,:,jp), E_uu) + kron(E_xsxf, dE_uu(:,:,jp))); - dVarinov(id_e6_u_xs , id_e4_u_xf , jp) = kron(dE_uu(:,:,jp), E_xsxf) + kron(E_uu, dE_xsxf(:,:,jp)); - dVarinov(id_e6_u_xs , id_e5_xs_u , jp) = K_ux*(kron(dE_xsxs(:,:,jp),E_uu) + kron(E_xsxs,dE_uu(:,:,jp))); - dVarinov(id_e6_u_xs , id_e6_u_xs , jp) = kron(dE_uu(:,:,jp), E_xsxs) + kron(E_uu, dE_xsxs(:,:,jp)); - dVarinov(id_e6_u_xs , id_e7_xf_xf_u , jp) = K_ux*(kron(reshape(dE_xs_xf_xf(:,jp),x_nbr,x_nbr^2), E_uu) + kron(reshape(E_xs_xf_xf,x_nbr,x_nbr^2), dE_uu(:,:,jp))); - dVarinov(id_e6_u_xs , id_e8_xf_u_xf , jp) = (kron(dE_uu(:,:,jp), reshape(E_xs_xf_xf,x_nbr,x_nbr^2)) + kron(E_uu, reshape(dE_xs_xf_xf(:,jp),x_nbr,x_nbr^2)))*kron(K_xu,speye(x_nbr))'; - dVarinov(id_e6_u_xs , id_e9_u_xf_xf , jp) = kron(dE_uu(:,:,jp), reshape(E_xs_xf_xf,x_nbr,x_nbr^2)) + kron(E_uu, reshape(dE_xs_xf_xf(:,jp),x_nbr,x_nbr^2)); - if jp <= (stderrparam_nbr+corrparam_nbr) - dVarinov(id_e6_u_xs , id_e13_u_u_u , jp) = K_ux*(kron(dE_xs(:,jp),reshape(E_u_u_u_u,u_nbr,u_nbr^3)) + kron(E_xs,reshape(QPu*dE_u_u_u_u(:,jp),u_nbr,u_nbr^3))); - else - dVarinov(id_e6_u_xs , id_e13_u_u_u , jp) = K_ux*kron(dE_xs(:,jp),reshape(E_u_u_u_u,u_nbr,u_nbr^3)); - end + dC(: , id_z1_xf , jp3) = dgx_jp3 + 3/6*dgxss_jp3; + dC(: , id_z2_xs , jp3) = dgx_jp3; + dC(: , id_z3_xf_xf , jp3) = 1/2*dgxx_jp3; + dC(: , id_z4_xrd , jp3) = dgx_jp3; + dC(: , id_z5_xf_xs , jp3) = dgxx_jp3; + dC(: , id_z6_xf_xf_xf , jp3) = 1/6*dgxxx_jp3; - dVarinov(id_e7_xf_xf_u , id_e1_u , jp) = kron(vec(dE_xfxf(:,:,jp)),E_uu) + kron(E_xfxf(:),dE_uu(:,:,jp)); - dVarinov(id_e7_xf_xf_u , id_e3_xf_u , jp) = kron(reshape(dE_xf_xf_xf(:,jp),x_nbr^2,x_nbr),E_uu) + kron(reshape(E_xf_xf_xf,x_nbr^2,x_nbr),dE_uu(:,:,jp)); - dVarinov(id_e7_xf_xf_u , id_e4_u_xf , jp) = (kron(reshape(dE_xf_xf_xf(:,jp),x_nbr^2,x_nbr),E_uu) + kron(reshape(E_xf_xf_xf,x_nbr^2,x_nbr),dE_uu(:,:,jp)))*K_ux'; - dVarinov(id_e7_xf_xf_u , id_e5_xs_u , jp) = kron(reshape(dE_xf_xf_xs(:,jp),x_nbr^2,x_nbr),E_uu) + kron(reshape(E_xf_xf_xs,x_nbr^2,x_nbr),dE_uu(:,:,jp)); - dVarinov(id_e7_xf_xf_u , id_e6_u_xs , jp) = (kron(reshape(dE_xf_xf_xs(:,jp),x_nbr^2,x_nbr),E_uu) + kron(reshape(E_xf_xf_xs,x_nbr^2,x_nbr),dE_uu(:,:,jp)))*K_ux'; - dVarinov(id_e7_xf_xf_u , id_e7_xf_xf_u , jp) = kron(reshape(dE_xf_xf_xf_xf(:,jp),x_nbr^2,x_nbr^2),E_uu) + kron(reshape(E_xf_xf_xf_xf,x_nbr^2,x_nbr^2),dE_uu(:,:,jp)); - dVarinov(id_e7_xf_xf_u , id_e8_xf_u_xf , jp) = (kron(reshape(dE_xf_xf_xf_xf(:,jp),x_nbr^2,x_nbr^2),E_uu) + kron(reshape(E_xf_xf_xf_xf,x_nbr^2,x_nbr^2),dE_uu(:,:,jp)))*kron(speye(x_nbr),K_ux)'; - dVarinov(id_e7_xf_xf_u , id_e9_u_xf_xf , jp) = kron(speye(x_nbr),K_ux)*kron(K_ux,speye(x_nbr))*(kron(dE_uu(:,:,jp), reshape(E_xf_xf_xf_xf,x_nbr^2,x_nbr^2)) + kron(E_uu, reshape(dE_xf_xf_xf_xf(:,jp),x_nbr^2,x_nbr^2))); - if jp <= (stderrparam_nbr+corrparam_nbr) - dVarinov(id_e7_xf_xf_u , id_e13_u_u_u , jp) = kron(vec(dE_xfxf(:,:,jp)),reshape(E_u_u_u_u,u_nbr,u_nbr^3)) + kron(E_xfxf(:),reshape(QPu*dE_u_u_u_u(:,jp),u_nbr,u_nbr^3)); - else - dVarinov(id_e7_xf_xf_u , id_e13_u_u_u , jp) = kron(vec(dE_xfxf(:,:,jp)),reshape(E_u_u_u_u,u_nbr,u_nbr^3)); - end + dD(: , id_inov1_u , jp3) = dgu_jp3 + 3/6*dguss_jp3; + dD(: , id_inov2_u_u , jp3) = 1/2*dguu_jp3; + dD(: , id_inov3_xf_u , jp3) = dgxu_jp3; + dD(: , id_inov4_xs_u , jp3) = dgxu_jp3; + dD(: , id_inov5_xf_xf_u , jp3) = 3/6*dgxxu_jp3; + dD(: , id_inov6_xf_u_u , jp3) = 3/6*dgxuu_jp3; + dD(: , id_inov7_u_u_u , jp3) = 1/6*dguuu_jp3; - dVarinov(id_e8_xf_u_xf , id_e1_u , jp) = kron(K_xu,speye(x_nbr))*(kron(dE_uu(:,:,jp),E_xfxf(:)) + kron(E_uu,vec(dE_xfxf(:,:,jp)))); - dVarinov(id_e8_xf_u_xf , id_e3_xf_u , jp) = kron(speye(x_nbr),K_xu)*(kron(reshape(dE_xf_xf_xf(:,jp),x_nbr^2,x_nbr),E_uu) + kron(reshape(E_xf_xf_xf,x_nbr^2,x_nbr),dE_uu(:,:,jp))); - dVarinov(id_e8_xf_u_xf , id_e4_u_xf , jp) = kron(K_xu,speye(x_nbr))*(kron(dE_uu(:,:,jp),reshape(E_xf_xf_xf,x_nbr^2,x_nbr)) + kron(E_uu,reshape(dE_xf_xf_xf(:,jp),x_nbr^2,x_nbr))); - dVarinov(id_e8_xf_u_xf , id_e5_xs_u , jp) = kron(speye(x_nbr),K_ux)*(kron(reshape(dE_xf_xf_xs(:,jp),x_nbr^2,x_nbr),E_uu) + kron(reshape(E_xf_xf_xs,x_nbr^2,x_nbr),dE_uu(:,:,jp))); - dVarinov(id_e8_xf_u_xf , id_e6_u_xs , jp) = kron(K_xu,speye(x_nbr))*(kron(dE_uu(:,:,jp),reshape(E_xf_xf_xs,x_nbr^2,x_nbr)) + kron(E_uu,reshape(dE_xf_xf_xs(:,jp),x_nbr^2,x_nbr))); - dVarinov(id_e8_xf_u_xf , id_e7_xf_xf_u , jp) = kron(speye(x_nbr),K_ux)*(kron(reshape(dE_xf_xf_xf_xf(:,jp),x_nbr^2,x_nbr^2),E_uu) + kron(reshape(E_xf_xf_xf_xf,x_nbr^2,x_nbr^2),dE_uu(:,:,jp))); - dVarinov(id_e8_xf_u_xf , id_e8_xf_u_xf , jp) = kron(K_xu,speye(x_nbr))*(kron(dE_uu(:,:,jp),reshape(E_xf_xf_xf_xf,x_nbr^2,x_nbr^2)) + kron(E_uu,reshape(dE_xf_xf_xf_xf(:,jp),x_nbr^2,x_nbr^2)))*kron(K_xu,speye(x_nbr))'; - dVarinov(id_e8_xf_u_xf , id_e9_u_xf_xf , jp) = kron(K_xu,speye(x_nbr))*(kron(dE_uu(:,:,jp),reshape(E_xf_xf_xf_xf,x_nbr^2,x_nbr^2)) + kron(E_uu,reshape(dE_xf_xf_xf_xf(:,jp),x_nbr^2,x_nbr^2))); - if jp <= (stderrparam_nbr+corrparam_nbr) - dVarinov(id_e8_xf_u_xf , id_e13_u_u_u , jp) = kron(K_xu,speye(x_nbr))*(kron(reshape(QPu*dE_u_u_u_u(:,jp),u_nbr,u_nbr^3),E_xfxf(:)) + kron(reshape(E_u_u_u_u,u_nbr,u_nbr^3),vec(dE_xfxf(:,:,jp)))); - else - dVarinov(id_e8_xf_u_xf , id_e13_u_u_u , jp) = kron(K_xu,speye(x_nbr))*kron(reshape(E_u_u_u_u,u_nbr,u_nbr^3),vec(dE_xfxf(:,:,jp))); - end + dd(:,jp3) = 1/2*dgss_jp3 + 1/2*dguu_jp3*E_uu(:) + 1/2*guu*dE_uu_jp3(:); - dVarinov(id_e9_u_xf_xf , id_e1_u , jp) = kron(dE_uu(:,:,jp), E_xfxf(:)) + kron(E_uu, vec(dE_xfxf(:,:,jp))); - dVarinov(id_e9_u_xf_xf , id_e3_xf_u , jp) = (kron(dE_uu(:,:,jp), reshape(E_xf_xf_xf,x_nbr^2,x_nbr)) + kron(E_uu, reshape(dE_xf_xf_xf(:,jp),x_nbr^2,x_nbr)))*K_xu'; - dVarinov(id_e9_u_xf_xf , id_e4_u_xf , jp) = kron(dE_uu(:,:,jp), reshape(E_xf_xf_xf,x_nbr^2,x_nbr)) + kron(E_uu, reshape(dE_xf_xf_xf(:,jp),x_nbr^2,x_nbr)); - dVarinov(id_e9_u_xf_xf , id_e5_xs_u , jp) = (kron(dE_uu(:,:,jp), reshape(E_xf_xf_xs,x_nbr^2,x_nbr)) + kron(E_uu, reshape(dE_xf_xf_xs(:,jp),x_nbr^2,x_nbr)))*K_xu'; - dVarinov(id_e9_u_xf_xf , id_e6_u_xs , jp) = kron(dE_uu(:,:,jp), reshape(E_xf_xf_xs,x_nbr^2,x_nbr)) + kron(E_uu, reshape(dE_xf_xf_xs(:,jp),x_nbr^2,x_nbr)); - dVarinov(id_e9_u_xf_xf , id_e7_xf_xf_u , jp) = kron(speye(x_nbr),K_ux)*kron(K_ux,speye(x_nbr))*(kron(reshape(dE_xf_xf_xf_xf(:,jp),x_nbr^2,x_nbr^2),E_uu) + kron(reshape(E_xf_xf_xf_xf,x_nbr^2,x_nbr^2),dE_uu(:,:,jp))); - dVarinov(id_e9_u_xf_xf , id_e8_xf_u_xf , jp) = (kron(dE_uu(:,:,jp),reshape(E_xf_xf_xf_xf,x_nbr^2,x_nbr^2)) + kron(E_uu,reshape(dE_xf_xf_xf_xf(:,jp),x_nbr^2,x_nbr^2)))*kron(speye(x_nbr),K_ux)'; - dVarinov(id_e9_u_xf_xf , id_e9_u_xf_xf , jp) = kron(dE_uu(:,:,jp),reshape(E_xf_xf_xf_xf,x_nbr^2,x_nbr^2)) + kron(E_uu,reshape(dE_xf_xf_xf_xf(:,jp),x_nbr^2,x_nbr^2)); - if jp <= (stderrparam_nbr+corrparam_nbr) - dVarinov(id_e9_u_xf_xf , id_e13_u_u_u , jp) = kron(reshape(QPu*dE_u_u_u_u(:,jp),u_nbr,u_nbr^3),E_xfxf(:)) + kron(reshape(E_u_u_u_u,u_nbr,u_nbr^3),vec(dE_xfxf(:,:,jp))); - else - dVarinov(id_e9_u_xf_xf , id_e13_u_u_u , jp) = kron(reshape(E_u_u_u_u,u_nbr,u_nbr^3),vec(dE_xfxf(:,:,jp))); - end + dVarinov(id_inov1_u , id_inov1_u , jp3) = dE_uu_jp3; + dVarinov(id_inov1_u , id_inov4_xs_u , jp3) = kron(dE_xs_jp3',E_uu) + kron(E_xs',dE_uu_jp3); + dVarinov(id_inov1_u , id_inov5_xf_xf_u , jp3) = kron(dE_xfxf_jp3(:)',E_uu) + kron(E_xfxf(:)',dE_uu_jp3); + dVarinov(id_inov1_u , id_inov7_u_u_u , jp3) = reshape(dE_u_u_u_u_jp3,u_nbr,u_nbr^3); + dVarinov(id_inov2_u_u , id_inov2_u_u , jp3) = reshape(dE_u_u_u_u_jp3,u_nbr^2,u_nbr^2) - dE_uu_jp3(:)*E_uu(:)' - E_uu(:)*dE_uu_jp3(:)'; + dVarinov(id_inov3_xf_u , id_inov3_xf_u , jp3) = dE_xfxf_uu_jp3; + dVarinov(id_inov4_xs_u , id_inov1_u , jp3) = kron(dE_xs_jp3,E_uu) + kron(E_xs,dE_uu_jp3); + dVarinov(id_inov4_xs_u , id_inov4_xs_u , jp3) = kron(dE_xsxs_jp3,E_uu) + kron(E_xsxs,dE_uu_jp3); + dVarinov(id_inov4_xs_u , id_inov5_xf_xf_u , jp3) = kron(dE_xsxf_xf_jp3, E_uu) + kron(E_xsxf_xf, dE_uu_jp3); + dVarinov(id_inov4_xs_u , id_inov7_u_u_u , jp3) = kron(dE_xs_jp3,reshape(QPu*E_u_u_u_u,u_nbr,u_nbr^3)) + kron(E_xs,reshape(dE_u_u_u_u_jp3,u_nbr,u_nbr^3)); + dVarinov(id_inov5_xf_xf_u , id_inov1_u , jp3) = kron(dE_xfxf_jp3(:),E_uu) + kron(E_xfxf(:),dE_uu_jp3); + dVarinov(id_inov5_xf_xf_u , id_inov4_xs_u , jp3) = kron(dE_xf_xfxs_jp3,E_uu) + kron(E_xf_xfxs,dE_uu_jp3); + dVarinov(id_inov5_xf_xf_u , id_inov5_xf_xf_u , jp3) = kron(dE_xf_xfxf_xf_jp3,E_uu) + kron(E_xf_xfxf_xf,dE_uu_jp3); + dVarinov(id_inov5_xf_xf_u , id_inov7_u_u_u , jp3) = kron(dE_xfxf_jp3(:),reshape(QPu*E_u_u_u_u,u_nbr,u_nbr^3)) + kron(E_xfxf(:),reshape(dE_u_u_u_u_jp3,u_nbr,u_nbr^3)); + dVarinov(id_inov6_xf_u_u , id_inov6_xf_u_u , jp3) = kron(dE_xfxf_jp3,reshape(QPu*E_u_u_u_u,u_nbr^2,u_nbr^2)) + kron(E_xfxf,reshape(dE_u_u_u_u_jp3,u_nbr^2,u_nbr^2)); + dVarinov(id_inov7_u_u_u , id_inov1_u , jp3) = reshape(dE_u_u_u_u_jp3,u_nbr^3,u_nbr); + dVarinov(id_inov7_u_u_u , id_inov4_xs_u , jp3) = kron(dE_xs_jp3',reshape(QPu*E_u_u_u_u,u_nbr^3,u_nbr)) + kron(E_xs',reshape(dE_u_u_u_u_jp3,u_nbr^3,u_nbr)); + dVarinov(id_inov7_u_u_u , id_inov5_xf_xf_u , jp3) = kron(transpose(dE_xfxf_jp3(:)),reshape(QPu*E_u_u_u_u,u_nbr^3,u_nbr)) + kron(transpose(E_xfxf(:)),reshape(dE_u_u_u_u_jp3,u_nbr^3,u_nbr)); + dVarinov(id_inov7_u_u_u , id_inov7_u_u_u , jp3) = reshape(dE_u_u_u_u_u_u_jp3,u_nbr^3,u_nbr^3); + + dE_inovzlag1(id_inov6_xf_u_u , id_z1_xf , jp3) = kron(dE_xfxf_jp3,E_uu(:)) + kron(E_xfxf,dE_uu_jp3(:)); + dE_inovzlag1(id_inov6_xf_u_u , id_z4_xrd , jp3) = kron(dE_xrdxf_jp3',E_uu(:)) + kron(E_xrdxf',dE_uu_jp3(:)); + dE_inovzlag1(id_inov6_xf_u_u , id_z5_xf_xs , jp3) = kron(reshape(commutation(x_nbr^2,x_nbr)*vec(dE_xsxf_xf_jp3),x_nbr,x_nbr^2),vec(E_uu)) + kron(reshape(commutation(x_nbr^2,x_nbr)*vec(E_xsxf_xf),x_nbr,x_nbr^2),vec(dE_uu_jp3)) ; + dE_inovzlag1(id_inov6_xf_u_u , id_z6_xf_xf_xf , jp3) = kron(reshape(dE_xf_xfxf_xf_jp3,x_nbr,x_nbr^3),E_uu(:)) + kron(reshape(E_xf_xfxf_xf,x_nbr,x_nbr^3),dE_uu_jp3(:)); + + dBinovzlag1A_jp3 = dB(:,:,jp3)*E_inovzlag1*transpose(A) + B*dE_inovzlag1(:,:,jp3)*transpose(A) + B*E_inovzlag1*transpose(dA(:,:,jp3)); + dOm_z_jp3 = dB(:,:,jp3)*Varinov*transpose(B) + B*dVarinov(:,:,jp3)*transpose(B) + B*Varinov*transpose(dB(:,:,jp3)) + dBinovzlag1A_jp3 + transpose(dBinovzlag1A_jp3); - if jp <= (stderrparam_nbr+corrparam_nbr) - dVarinov(id_e10_xf_u_u , id_e10_xf_u_u , jp) = kron(dE_xfxf(:,:,jp),reshape(E_u_u_u_u,u_nbr^2,u_nbr^2)) + kron(E_xfxf,reshape(QPu*dE_u_u_u_u(:,jp),u_nbr^2,u_nbr^2)); - dVarinov(id_e10_xf_u_u , id_e11_u_xf_u , jp) = (kron(dE_xfxf(:,:,jp),reshape(E_u_u_u_u,u_nbr^2,u_nbr^2)) + kron(E_xfxf,reshape(QPu*dE_u_u_u_u(:,jp),u_nbr^2,u_nbr^2)))*kron(K_ux,speye(u_nbr))'; - dVarinov(id_e10_xf_u_u , id_e12_u_u_xf , jp) = (kron(dE_xfxf(:,:,jp),reshape(E_u_u_u_u,u_nbr^2,u_nbr^2)) + kron(E_xfxf,reshape(QPu*dE_u_u_u_u(:,jp),u_nbr^2,u_nbr^2)))*kron(K_ux,speye(u_nbr))'*kron(speye(u_nbr),K_ux)'; - else - dVarinov(id_e10_xf_u_u , id_e10_xf_u_u , jp) = kron(dE_xfxf(:,:,jp),reshape(E_u_u_u_u,u_nbr^2,u_nbr^2)); - dVarinov(id_e10_xf_u_u , id_e11_u_xf_u , jp) = kron(dE_xfxf(:,:,jp),reshape(E_u_u_u_u,u_nbr^2,u_nbr^2))*kron(K_ux,speye(u_nbr))'; - dVarinov(id_e10_xf_u_u , id_e12_u_u_xf , jp) = kron(dE_xfxf(:,:,jp),reshape(E_u_u_u_u,u_nbr^2,u_nbr^2))*kron(K_ux,speye(u_nbr))'*kron(speye(u_nbr),K_ux)'; + [dVar_z(:,:,jp3), errorflag] = disclyap_fast(A, dA(:,:,jp3)*Var_z*A' + A*Var_z*dA(:,:,jp3)' + dOm_z_jp3, options.lyapunov_doubling_tol); + if errorflag + dVar_z(:,:,jp3) = lyapunov_symm(A, dA(:,:,jp3)*Var_z*A' + A*Var_z*dA(:,:,jp3)' + dOm_z_jp3,... + options.lyapunov_fixed_point_tol,options.qz_criterium,options.lyapunov_complex_threshold,... + lyapunov_symm_method,... + options.debug); + if lyapunov_symm_method == 1 + lyapunov_symm_method = 2; %now we can reuse persistent schur + end end - - if jp <= (stderrparam_nbr+corrparam_nbr) - dVarinov(id_e11_u_xf_u , id_e10_xf_u_u , jp) = kron(K_ux,speye(u_nbr))*(kron(dE_xfxf(:,:,jp),reshape(E_u_u_u_u,u_nbr^2,u_nbr^2)) + kron(E_xfxf,reshape(QPu*dE_u_u_u_u(:,jp),u_nbr^2,u_nbr^2))); - dVarinov(id_e11_u_xf_u , id_e11_u_xf_u , jp) = kron(K_ux,speye(u_nbr))*(kron(dE_xfxf(:,:,jp),reshape(E_u_u_u_u,u_nbr^2,u_nbr^2)) + kron(E_xfxf,reshape(QPu*dE_u_u_u_u(:,jp),u_nbr^2,u_nbr^2)))*kron(K_ux,speye(u_nbr))'; - dVarinov(id_e11_u_xf_u , id_e12_u_u_xf , jp) = kron(speye(u_nbr),K_ux)*(kron(reshape(QPu*dE_u_u_u_u(:,jp),u_nbr^2,u_nbr^2),E_xfxf) + kron(reshape(E_u_u_u_u,u_nbr^2,u_nbr^2),dE_xfxf(:,:,jp))); - else - dVarinov(id_e11_u_xf_u , id_e10_xf_u_u , jp) = kron(K_ux,speye(u_nbr))*kron(dE_xfxf(:,:,jp),reshape(E_u_u_u_u,u_nbr^2,u_nbr^2)); - dVarinov(id_e11_u_xf_u , id_e11_u_xf_u , jp) = kron(K_ux,speye(u_nbr))*kron(dE_xfxf(:,:,jp),reshape(E_u_u_u_u,u_nbr^2,u_nbr^2))*kron(K_ux,speye(u_nbr))'; - dVarinov(id_e11_u_xf_u , id_e12_u_u_xf , jp) = kron(speye(u_nbr),K_ux)*kron(reshape(E_u_u_u_u,u_nbr^2,u_nbr^2),dE_xfxf(:,:,jp)); - end - - if jp <= (stderrparam_nbr+corrparam_nbr) - dVarinov(id_e12_u_u_xf , id_e10_xf_u_u , jp) = kron(K_ux,speye(u_nbr))*kron(speye(u_nbr),K_ux)*(kron(dE_xfxf(:,:,jp),reshape(E_u_u_u_u,u_nbr^2,u_nbr^2)) + kron(E_xfxf,reshape(QPu*dE_u_u_u_u(:,jp),u_nbr^2,u_nbr^2))); - dVarinov(id_e12_u_u_xf , id_e11_u_xf_u , jp) = (kron(reshape(QPu*dE_u_u_u_u(:,jp),u_nbr^2,u_nbr^2),E_xfxf) + kron(reshape(E_u_u_u_u,u_nbr^2,u_nbr^2),dE_xfxf(:,:,jp)))*kron(speye(u_nbr),K_xu)'; - dVarinov(id_e12_u_u_xf , id_e12_u_u_xf , jp) = kron(reshape(QPu*dE_u_u_u_u(:,jp),u_nbr^2,u_nbr^2),E_xfxf) + kron(reshape(E_u_u_u_u,u_nbr^2,u_nbr^2),dE_xfxf(:,:,jp)); - else - dVarinov(id_e12_u_u_xf , id_e10_xf_u_u , jp) = kron(K_ux,speye(u_nbr))*kron(speye(u_nbr),K_ux)*kron(dE_xfxf(:,:,jp),reshape(E_u_u_u_u,u_nbr^2,u_nbr^2)); - dVarinov(id_e12_u_u_xf , id_e11_u_xf_u , jp) = kron(reshape(E_u_u_u_u,u_nbr^2,u_nbr^2),dE_xfxf(:,:,jp))*kron(speye(u_nbr),K_xu)'; - dVarinov(id_e12_u_u_xf , id_e12_u_u_xf , jp) = kron(reshape(E_u_u_u_u,u_nbr^2,u_nbr^2),dE_xfxf(:,:,jp)); - end - - if jp <= (stderrparam_nbr+corrparam_nbr) - dVarinov(id_e13_u_u_u , id_e1_u , jp) = reshape(QPu*dE_u_u_u_u(:,jp),u_nbr^3,u_nbr); - dVarinov(id_e13_u_u_u , id_e5_xs_u , jp) = kron(dE_xs(:,jp)', reshape(E_u_u_u_u,u_nbr^3,u_nbr)) + kron(E_xs', reshape(QPu*dE_u_u_u_u(:,jp),u_nbr^3,u_nbr)); - dVarinov(id_e13_u_u_u , id_e6_u_xs , jp) = kron(reshape(QPu*dE_u_u_u_u(:,jp),u_nbr^3,u_nbr),E_xs') + kron(reshape(E_u_u_u_u,u_nbr^3,u_nbr),dE_xs(:,jp)'); - dVarinov(id_e13_u_u_u , id_e7_xf_xf_u , jp) = kron(vec(dE_xfxf(:,:,jp))',reshape(E_u_u_u_u,u_nbr^3,u_nbr)) + kron(E_xfxf(:)',reshape(QPu*dE_u_u_u_u(:,jp),u_nbr^3,u_nbr)); - dVarinov(id_e13_u_u_u , id_e8_xf_u_xf , jp) = (kron(vec(dE_xfxf(:,:,jp))',reshape(E_u_u_u_u,u_nbr^3,u_nbr)) + kron(E_xfxf(:)',reshape(QPu*dE_u_u_u_u(:,jp),u_nbr^3,u_nbr)))*kron(speye(x_nbr),K_ux)'; - dVarinov(id_e13_u_u_u , id_e9_u_xf_xf , jp) = kron(reshape(QPu*dE_u_u_u_u(:,jp),u_nbr^3,u_nbr), E_xfxf(:)') + kron(reshape(E_u_u_u_u,u_nbr^3,u_nbr), vec(dE_xfxf(:,:,jp))'); - else - dVarinov(id_e13_u_u_u , id_e5_xs_u , jp) = kron(dE_xs(:,jp)', reshape(E_u_u_u_u,u_nbr^3,u_nbr)); - dVarinov(id_e13_u_u_u , id_e6_u_xs , jp) = kron(reshape(E_u_u_u_u,u_nbr^3,u_nbr),dE_xs(:,jp)'); - dVarinov(id_e13_u_u_u , id_e7_xf_xf_u , jp) = kron(vec(dE_xfxf(:,:,jp))',reshape(E_u_u_u_u,u_nbr^3,u_nbr)); - dVarinov(id_e13_u_u_u , id_e8_xf_u_xf , jp) = kron(vec(dE_xfxf(:,:,jp))',reshape(E_u_u_u_u,u_nbr^3,u_nbr))*kron(speye(x_nbr),K_ux)'; - dVarinov(id_e13_u_u_u , id_e9_u_xf_xf , jp) = kron(reshape(E_u_u_u_u,u_nbr^3,u_nbr), vec(dE_xfxf(:,:,jp))'); - end - if jp <= (stderrparam_nbr+corrparam_nbr) - dVarinov(id_e13_u_u_u , id_e13_u_u_u , jp) = reshape(Q6Pu*dE_u_u_u_u_u_u(:,jp),u_nbr^3,u_nbr^3); - end - - dE_z(:,jp) = (speye(z_nbr)-A)\(dc(:,jp) + dA(:,:,jp)*E_z); + %make sure some stuff is zero due to Gaussianity + dVar_z(id_z1_xf , id_z2_xs , jp3) = zeros(x_nbr,x_nbr); + dVar_z(id_z1_xf , id_z3_xf_xf , jp3) = zeros(x_nbr,x_nbr^2); + dVar_z(id_z2_xs , id_z1_xf , jp3) = zeros(x_nbr,x_nbr); + dVar_z(id_z2_xs , id_z4_xrd , jp3) = zeros(x_nbr,x_nbr); + dVar_z(id_z2_xs , id_z5_xf_xs , jp3) = zeros(x_nbr,x_nbr^2); + dVar_z(id_z2_xs , id_z6_xf_xf_xf , jp3) = zeros(x_nbr,x_nbr^3); + dVar_z(id_z3_xf_xf , id_z1_xf , jp3) = zeros(x_nbr^2,x_nbr); + dVar_z(id_z3_xf_xf , id_z4_xrd , jp3) = zeros(x_nbr^2,x_nbr); + dVar_z(id_z3_xf_xf , id_z5_xf_xs , jp3) = zeros(x_nbr^2,x_nbr^2); + dVar_z(id_z3_xf_xf , id_z6_xf_xf_xf , jp3) = zeros(x_nbr^2,x_nbr^3); + dVar_z(id_z4_xrd , id_z2_xs , jp3) = zeros(x_nbr,x_nbr); + dVar_z(id_z4_xrd , id_z3_xf_xf , jp3) = zeros(x_nbr,x_nbr^2); + dVar_z(id_z5_xf_xs , id_z2_xs , jp3) = zeros(x_nbr^2,x_nbr); + dVar_z(id_z5_xf_xs , id_z3_xf_xf , jp3) = zeros(x_nbr^2,x_nbr^2); + dVar_z(id_z6_xf_xf_xf , id_z2_xs , jp3) = zeros(x_nbr^3,x_nbr); + dVar_z(id_z6_xf_xf_xf , id_z3_xf_xf , jp3) = zeros(x_nbr^3,x_nbr^2); end end end end +%% Covariance/Correlation of control variables +Var_y = NaN*ones(y_nbr,y_nbr); +if order < 3 + Var_y(stationary_vars,stationary_vars) = C(stationary_vars,:)*Var_z*C(stationary_vars,:)'... + + D(stationary_vars,:)*Varinov*D(stationary_vars,:)'; +else + Var_y(stationary_vars,stationary_vars) = C(stationary_vars,:)*Var_z*C(stationary_vars,:)'... + + D(stationary_vars,:)*E_inovzlag1*C(stationary_vars,:)'... + + C(stationary_vars,:)*transpose(E_inovzlag1)*D(stationary_vars,:)'... + + D(stationary_vars,:)*Varinov*D(stationary_vars,:)'; +end +indzeros = find(abs(Var_y) < 1e-12); %find values that are numerical zero +Var_y(indzeros) = 0; +if useautocorr + sy = sqrt(diag(Var_y)); %theoretical standard deviation + sy = sy(stationary_vars); + sy = sy*sy'; %cross products of standard deviations + Corr_y = NaN*ones(y_nbr,y_nbr); + Corr_y(stationary_vars,stationary_vars) = Var_y(stationary_vars,stationary_vars)./sy; + Corr_yi = NaN*ones(y_nbr,y_nbr,nlags); +end -E_y = Yss(indy,:) + C*E_z + d; -Om_z = B*Varinov*transpose(B); -Om_y = D*Varinov*transpose(D); +if compute_derivs + dVar_y = NaN*ones(y_nbr,y_nbr,totparam_nbr); + if useautocorr + dCorr_y = NaN*ones(y_nbr,y_nbr,totparam_nbr); + dCorr_yi = NaN*ones(y_nbr,y_nbr,nlags,totparam_nbr); + end + for jpV=1:totparam_nbr + if order < 3 + dVar_y(stationary_vars,stationary_vars,jpV) = dC(stationary_vars,:,jpV)*Var_z*C(stationary_vars,:)' + C(stationary_vars,:)*dVar_z(:,:,jpV)*C(stationary_vars,:)' + C(stationary_vars,:)*Var_z*dC(stationary_vars,:,jpV)'... + + dD(stationary_vars,:,jpV)*Varinov*D(stationary_vars,:)' + D(stationary_vars,:)*dVarinov(:,:,jpV)*D(stationary_vars,:)' + D(stationary_vars,:)*Varinov*dD(stationary_vars,:,jpV)'; + else + dVar_y(stationary_vars,stationary_vars,jpV) = dC(stationary_vars,:,jpV)*Var_z*C(stationary_vars,:)' + C(stationary_vars,:)*dVar_z(:,:,jpV)*C(stationary_vars,:)' + C(stationary_vars,:)*Var_z*dC(stationary_vars,:,jpV)'... + + dD(stationary_vars,:,jpV)*E_inovzlag1*C(stationary_vars,:)' + D(stationary_vars,:)*dE_inovzlag1(:,:,jpV)*C(stationary_vars,:)' + D(stationary_vars,:)*E_inovzlag1*dC(stationary_vars,:,jpV)'... + + dC(stationary_vars,:,jpV)*transpose(E_inovzlag1)*D(stationary_vars,:)' + C(stationary_vars,:)*transpose(dE_inovzlag1(:,:,jpV))*D(stationary_vars,:)' + C(stationary_vars,:)*transpose(E_inovzlag1)*dD(stationary_vars,:,jpV)'... + + dD(stationary_vars,:,jpV)*Varinov*D(stationary_vars,:)' + D(stationary_vars,:)*dVarinov(:,:,jpV)*D(stationary_vars,:)' + D(stationary_vars,:)*Varinov*dD(stationary_vars,:,jpV)'; + end + [indzerosrow,indzeroscol] = find(abs(dVar_y(:,:,jpV)) < 1e-12); %find values that are numerical zero + dVar_y(indzerosrow,indzeroscol,jpV) = 0; + if useautocorr + %is this correct?[@wmutschl] + dsy = 1/2./sy.*diag(dVar_y(:,:,jpV)); + dsy = dsy(stationary_vars); + dsy = dsy*sy'+sy*dsy'; + dCorr_y(stationary_vars,stationary_vars,jpV) = (dVar_y(stationary_vars,stationary_vars,jpV).*sy-dsy.*Var_y(stationary_vars,stationary_vars))./(sy.*sy); + dCorr_y(stationary_vars,stationary_vars,jpV) = dCorr_y(stationary_vars,stationary_vars,jpV)-diag(diag(dCorr_y(stationary_vars,stationary_vars,jpV)))+diag(diag(dVar_y(stationary_vars,stationary_vars,jpV))); + end + end +end + +%% Autocovariances/autocorrelations of lagged control variables +Var_yi = NaN*ones(y_nbr,y_nbr,nlags); +Ai = eye(z_nbr); %this is A^0 +hxi = eye(x_nbr); +E_inovzlagi = E_inovzlag1; +Var_zi = Var_z; +if order <= 2 + tmp = A*Var_z*C(stationary_vars,:)' + B*Varinov*D(stationary_vars,:)'; +else + tmp = A*E_inovzlag1'*D(stationary_vars,:)' + B*Varinov*D(stationary_vars,:)'; +end +for i = 1:nlags + if order <= 2 + Var_yi(stationary_vars,stationary_vars,i) = C(stationary_vars,:)*Ai*tmp; + else + Var_zi = A*Var_zi + B*E_inovzlagi; + hxi = hx*hxi; + E_inovzlagi = zeros(inov_nbr,z_nbr); + E_inovzlagi(id_inov6_xf_u_u , id_z1_xf ) = kron(hxi*E_xfxf,E_uu(:)); + E_inovzlagi(id_inov6_xf_u_u , id_z4_xrd ) = kron(hxi*E_xrdxf',E_uu(:)); + E_inovzlagi(id_inov6_xf_u_u , id_z5_xf_xs ) = kron(hxi*reshape(commutation(x_nbr^2,x_nbr)*vec(E_xsxf_xf),x_nbr,x_nbr^2),vec(E_uu)); + E_inovzlagi(id_inov6_xf_u_u , id_z6_xf_xf_xf ) = kron(hxi*reshape(E_xf_xfxf_xf,x_nbr,x_nbr^3),E_uu(:)); + Var_yi(stationary_vars,stationary_vars,i) = C(stationary_vars,:)*Var_zi*C(stationary_vars,:)' + C(stationary_vars,:)*Ai*tmp + D(stationary_vars,:)*E_inovzlagi*C(stationary_vars,:)'; + end + if useautocorr + Corr_yi(stationary_vars,stationary_vars,i) = Var_yi(stationary_vars,stationary_vars,i)./sy; + end + Ai = Ai*A; %note that this is A^(i-1) +end + +if compute_derivs + dVar_yi = NaN*ones(y_nbr,y_nbr,nlags,totparam_nbr); + for jpVi=1:totparam_nbr + Ai = eye(z_nbr); dAi_jpVi = zeros(z_nbr,z_nbr); + hxi = eye(x_nbr); dhxi_jpVi = zeros(x_nbr,x_nbr); + E_inovzlagi = E_inovzlag1; dE_inovzlagi_jpVi = dE_inovzlag1(:,:,jpVi); + Var_zi = Var_z; dVar_zi_jpVi = dVar_z(:,:,jpVi); + if order <= 2 + dtmp_jpVi = dA(:,:,jpVi)*Var_z*C(stationary_vars,:)' + A*dVar_z(:,:,jpVi)*C(stationary_vars,:)' + A*Var_z*dC(stationary_vars,:,jpVi)'... + + dB(:,:,jpVi)*Varinov*D(stationary_vars,:)' + B*dVarinov(:,:,jpVi)*D(stationary_vars,:)' + B*Varinov*dD(stationary_vars,:,jpVi)'; + else + dtmp_jpVi = dA(:,:,jpVi)*E_inovzlag1'*D(stationary_vars,:)' + A*dE_inovzlag1(:,:,jpVi)'*D(stationary_vars,:)' + A*E_inovzlag1'*dD(stationary_vars,:,jpVi)'... + + dB(:,:,jpVi)*Varinov*D(stationary_vars,:)' + B*dVarinov(:,:,jpVi)*D(stationary_vars,:)' + B*Varinov*dD(stationary_vars,:,jpVi)'; + end + + for i = 1:nlags + if order <= 2 + dVar_yi(stationary_vars,stationary_vars,i,jpVi) = dC(stationary_vars,:,jpVi)*Ai*tmp + C(stationary_vars,:)*dAi_jpVi*tmp + C(stationary_vars,:)*Ai*dtmp_jpVi; + else + Var_zi = A*Var_zi + B*E_inovzlagi; + dVar_zi_jpVi = dA(:,:,jpVi)*Var_zi + A*dVar_zi_jpVi + dB(:,:,jpVi)*E_inovzlagi + + B*dE_inovzlagi_jpVi; + dhxi_jpVi = dhx(:,:,jpVi)*hxi + hx*dhxi_jpVi; + hxi = hx*hxi; + E_inovzlagi = zeros(inov_nbr,z_nbr); + E_inovzlagi(id_inov6_xf_u_u , id_z1_xf ) = kron(hxi*E_xfxf,E_uu(:)); + E_inovzlagi(id_inov6_xf_u_u , id_z4_xrd ) = kron(hxi*E_xrdxf',E_uu(:)); + E_inovzlagi(id_inov6_xf_u_u , id_z5_xf_xs ) = kron(hxi*reshape(commutation(x_nbr^2,x_nbr)*vec(E_xsxf_xf),x_nbr,x_nbr^2),vec(E_uu)); + E_inovzlagi(id_inov6_xf_u_u , id_z6_xf_xf_xf ) = kron(hxi*reshape(E_xf_xfxf_xf,x_nbr,x_nbr^3),E_uu(:)); + dE_inovzlagi_jpVi = zeros(inov_nbr,z_nbr); + dE_inovzlagi_jpVi(id_inov6_xf_u_u , id_z1_xf ) = kron(dhxi_jpVi*E_xfxf,E_uu(:)) + kron(hxi*dE_xfxf(:,:,jpVi),E_uu(:)) + kron(hxi*E_xfxf,vec(dE_uu(:,:,jpVi))); + dE_inovzlagi_jpVi(id_inov6_xf_u_u , id_z4_xrd ) = kron(dhxi_jpVi*E_xrdxf',E_uu(:)) + kron(hxi*dE_xrdxf(:,:,jpVi)',E_uu(:)) + kron(hxi*E_xrdxf',vec(dE_uu(:,:,jpVi))); + dE_inovzlagi_jpVi(id_inov6_xf_u_u , id_z5_xf_xs ) = kron(dhxi_jpVi*reshape(commutation(x_nbr^2,x_nbr)*vec(E_xsxf_xf),x_nbr,x_nbr^2),vec(E_uu)) + kron(hxi*reshape(commutation(x_nbr^2,x_nbr)*vec(dE_xsxf_xf(:,:,jpVi)),x_nbr,x_nbr^2),vec(E_uu)) + kron(hxi*reshape(commutation(x_nbr^2,x_nbr)*vec(E_xsxf_xf),x_nbr,x_nbr^2),vec(dE_uu(:,:,jpVi))); + dE_inovzlagi_jpVi(id_inov6_xf_u_u , id_z6_xf_xf_xf ) = kron(dhxi_jpVi*reshape(E_xf_xfxf_xf,x_nbr,x_nbr^3),E_uu(:)) + kron(hxi*reshape(dE_xf_xfxf_xf(:,:,jpVi),x_nbr,x_nbr^3),E_uu(:)) + kron(hxi*reshape(E_xf_xfxf_xf,x_nbr,x_nbr^3),vec(dE_uu(:,:,jpVi))); + dVar_yi(stationary_vars,stationary_vars,i,jpVi) = dC(stationary_vars,:,jpVi)*Var_zi*C(stationary_vars,:)' + C(stationary_vars,:)*dVar_zi_jpVi*C(stationary_vars,:)' + C(stationary_vars,:)*Var_zi*dC(stationary_vars,:,jpVi)'... + + dC(stationary_vars,:,jpVi)*Ai*tmp + C(stationary_vars,:)*dAi_jpVi*tmp + C(stationary_vars,:)*Ai*dtmp_jpVi... + + dD(stationary_vars,:,jpVi)*E_inovzlagi*C(stationary_vars,:)' + D(stationary_vars,:)*dE_inovzlagi_jpVi*C(stationary_vars,:)' + D(stationary_vars,:)*E_inovzlagi*dC(stationary_vars,:,jpVi)'; + end + if useautocorr + dCorr_yi(stationary_vars,stationary_vars,i,jpVi) = (dVar_yi(stationary_vars,stationary_vars,i,jpVi).*sy-dsy.*Var_yi(stationary_vars,stationary_vars,i))./(sy.*sy); + end + dAi_jpVi = dAi_jpVi*A + Ai*dA(:,:,jpVi); + Ai = Ai*A; + end + end +end + + +%% Mean of control variables +E_z = E_xf; +if order > 1 + E_z = [E_xf;E_xs;E_xfxf(:)]; +end +if order > 2 + E_xf_xs = zeros(x_nbr^2,1); + E_xf_xf_xf = zeros(x_nbr^3,1); + E_z = [E_xf;E_xs;E_xfxf(:);E_xrd;E_xf_xs;E_xf_xf_xf]; +end +E_y = Yss(indy,:) + C*E_z + d; if compute_derivs dE_y = zeros(y_nbr,totparam_nbr); - dOm_z = zeros(z_nbr,z_nbr,totparam_nbr); - dOm_y = zeros(y_nbr,y_nbr,totparam_nbr); - for jp = 1:totparam_nbr - dE_y(:,jp) = dC(:,:,jp)*E_z + C*dE_z(:,jp) + dd(:,jp); - dOm_z(:,:,jp) = dB(:,:,jp)*Varinov*B' + B*dVarinov(:,:,jp)*B' + B*Varinov*dB(:,:,jp)'; - dOm_y(:,:,jp) = dD(:,:,jp)*Varinov*D' + D*dVarinov(:,:,jp)*D' + D*Varinov*dD(:,:,jp)'; - if jp > (stderrparam_nbr+corrparam_nbr) - dE_y(:,jp) = dE_y(:,jp) + dYss(indy,jp-stderrparam_nbr-corrparam_nbr); %add steady state + for jpE = 1:totparam_nbr + if order == 1 + dE_z_jpE = dE_xf(:,jpE); + elseif order == 2 + dE_z_jpE = [dE_xf(:,jpE);dE_xs(:,jpE);vec(dE_xfxf(:,:,jpE))]; + elseif order == 3 + dE_xf_xs_jpE = zeros(x_nbr^2,1); + dE_xf_xf_xf_jpE = zeros(x_nbr^3,1); + dE_z_jpE = [dE_xf(:,jpE);dE_xs(:,jpE);vec(dE_xfxf(:,:,jpE)); dE_xrd(:,jpE); dE_xf_xs_jpE; dE_xf_xf_xf_jpE]; + end + dE_y(:,jpE) = dC(:,:,jpE)*E_z + C*dE_z_jpE + dd(:,jpE); + if jpE > (stderrparam_nbr+corrparam_nbr) + dE_y(:,jpE) = dE_y(:,jpE) + dYss(indy,jpE-stderrparam_nbr-corrparam_nbr); %add steady state end end end - +non_stationary_vars = setdiff(1:y_nbr,stationary_vars); +E_y(non_stationary_vars) = NaN; +if compute_derivs + dE_y(non_stationary_vars,:) = NaN; +end %% Store into output structure -dr.pruned.indx = indx; -dr.pruned.indy = indy; -%dr.pruned.E_xfxf = E_xfxf; -dr.pruned.A = A; -dr.pruned.B = B; -dr.pruned.C = C; -dr.pruned.D = D; -dr.pruned.c = c; -dr.pruned.d = d; -dr.pruned.Om_z = Om_z; -dr.pruned.Om_y = Om_y; -dr.pruned.Varinov = Varinov; -dr.pruned.E_z = E_z; -dr.pruned.E_y = E_y; -if compute_derivs == 1 - %dr.pruned.dE_xfxf = dE_xfxf; - dr.pruned.dA = dA; - dr.pruned.dB = dB; - dr.pruned.dC = dC; - dr.pruned.dD = dD; - dr.pruned.dc = dc; - dr.pruned.dd = dd; - dr.pruned.dOm_z = dOm_z; - dr.pruned.dOm_y = dOm_y; - dr.pruned.dVarinov = dVarinov; - dr.pruned.dE_z = dE_z; - dr.pruned.dE_y = dE_y; +pruned_state_space.indx = indx; +pruned_state_space.indy = indy; +pruned_state_space.A = A; +pruned_state_space.B = B; +pruned_state_space.C = C; +pruned_state_space.D = D; +pruned_state_space.c = c; +pruned_state_space.d = d; +pruned_state_space.Varinov = Varinov; +pruned_state_space.Var_z = Var_z; %remove in future [@wmutschl] +pruned_state_space.Var_y = Var_y; +pruned_state_space.Var_yi = Var_yi; +if useautocorr + pruned_state_space.Corr_y = Corr_y; + pruned_state_space.Corr_yi = Corr_yi; +end +pruned_state_space.E_y = E_y; + +if compute_derivs == 1 + pruned_state_space.dA = dA; + pruned_state_space.dB = dB; + pruned_state_space.dC = dC; + pruned_state_space.dD = dD; + pruned_state_space.dc = dc; + pruned_state_space.dd = dd; + pruned_state_space.dVarinov = dVarinov; + pruned_state_space.dVar_y = dVar_y; + pruned_state_space.dVar_yi = dVar_yi; + if useautocorr + pruned_state_space.dCorr_y = dCorr_y; + pruned_state_space.dCorr_yi = dCorr_yi; + end + pruned_state_space.dE_y = dE_y; end diff --git a/matlab/quadruplication.m b/matlab/quadruplication.m index 462857ff0..2d47fe740 100644 --- a/matlab/quadruplication.m +++ b/matlab/quadruplication.m @@ -1,46 +1,54 @@ -% 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. 112–134 +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. +% +% 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 . +% ========================================================================= np = p*(p+1)*(p+2)*(p+3)/24; - -if sparseflag - QP = spalloc(p^4,p*(p+1)*(p+2)*(p+3)/24,p^4); -else - QP = zeros(p^4,p*(p+1)*(p+2)*(p+3)/24); -end +QP = spalloc(p^4,np,p^4); if nargout > 1 - if sparseflag - QPinv = spalloc(p*(p+1)*(p+2)*(p+3)/24,p*(p+1)*(p+2)*(p+3)/24,p^4); - else - QPinv = zeros(p*(p+1)*(p+2)*(p+3)/24,p*(p+1)*(p+2)*(p+3)/24); - end + QPinv = spalloc(np,np,p^4); end - +counti = 1; for l=1:p for k=l:p for j=k:p for i=j:p - if progress && (rem(counti,100)== 0) - msg = sprintf(' Quadruplication Matrix Processed %d/%d', counti, np); fprintf([reverseStr, msg]); reverseStr = repmat(sprintf('\b'), 1, length(msg)); - elseif progress && (counti==np) - msg = sprintf(' Quadruplication Matrix Processed %d/%d\n', counti, np); fprintf([reverseStr, msg]); reverseStr = repmat(sprintf('\b'), 1, length(msg)); - end idx = uperm([i j k l]); for r = 1:size(idx,1) ii = idx(r,1); jj= idx(r,2); kk=idx(r,3); ll=idx(r,4); @@ -75,6 +83,7 @@ end %QPinv = (transpose(QP)*QP)\transpose(QP); function m = mue(p,i,j,k,l) + % Auxiliary expression, see page 118 of Meijer (2005) 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; m = round(m); end diff --git a/matlab/uperm.m b/matlab/uperm.m index 0f4161b61..3322af96b 100644 --- a/matlab/uperm.m +++ b/matlab/uperm.m @@ -1,4 +1,5 @@ function p = uperm(a) +% Return all unique permutations of possibly-repeating array elements % ========================================================================= % Copyright (C) 2014 Bruno Luong % Copyright (C) 2020 Dynare Team diff --git a/tests/.gitignore b/tests/.gitignore index e7925b2d4..35ff3c844 100644 --- a/tests/.gitignore +++ b/tests/.gitignore @@ -126,6 +126,7 @@ wsOct !/practicing/dataHST.mat !/practicing/data_consRicardoypg.mat !/practicing/datasaver.m +!/pruning/Andreasen_et_al_2018_Dynare44Pruning_v2.mat !/utils/printMakeCheckMatlabErrMsg.m !/utils/printMakeCheckOctaveErrMsg.m !/prior_posterior_function/posterior_function_demo.m diff --git a/tests/Makefile.am b/tests/Makefile.am index b957098e9..7648c3818 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -21,6 +21,7 @@ MODFILES = \ analytic_derivatives/fs2000_analytic_derivation.mod \ analytic_derivatives/BrockMirman_PertParamsDerivs.mod \ analytic_derivatives/burnside_3_order_PertParamsDerivs.mod \ + pruning/AnSchorfheide_pruned_state_space.mod \ measurement_errors/fs2000_corr_me_ml_mcmc/fs2000_corr_ME.mod \ TeX/fs2000_corr_ME.mod \ estimation/MH_recover/fs2000_recover_tarb.mod \ @@ -716,6 +717,10 @@ analytic_derivatives: m/analytic_derivatives o/analytic_derivatives m/analytic_derivatives: $(patsubst %.mod, %.m.trs, $(filter analytic_derivatives/%.mod, $(MODFILES))) o/analytic_derivatives: $(patsubst %.mod, %.o.trs, $(filter analytic_derivatives/%.mod, $(MODFILES))) +pruning: m/pruning o/pruning +m/pruning: $(patsubst %.mod, %.m.trs, $(filter pruning/%.mod, $(MODFILES))) +o/pruning: $(patsubst %.mod, %.o.trs, $(filter pruning/%.mod, $(MODFILES))) + fs2000: m/fs2000 o/fs2000 m/fs2000: $(patsubst %.mod, %.m.trs, $(filter fs2000/%.mod, $(MODFILES))) o/fs2000: $(patsubst %.mod, %.o.trs, $(filter fs2000/%.mod, $(MODFILES))) @@ -837,6 +842,7 @@ EXTRA_DIST = \ gsa/data_ca1.mat \ analytic_derivatives/fsdat_simul.m \ analytic_derivatives/nBrockMirmanSYM.mat \ + pruning/Andreasen_et_al_2018_Dynare44Pruning_v2.mat \ fs2000/fsdat_simul.m \ fs2000/fsdat_simul_dseries.m \ fs2000/fsdat_simul_missing_obs.m \ diff --git a/tests/pruning/AnSchorfheide_pruned_state_space.mod b/tests/pruning/AnSchorfheide_pruned_state_space.mod new file mode 100644 index 000000000..9afd133b8 --- /dev/null +++ b/tests/pruning/AnSchorfheide_pruned_state_space.mod @@ -0,0 +1,193 @@ +% This mod file compares the functionality of Dynare's pruned_state_space.m with the +% external Dynare pruning toolbox of Andreasen, Fernández-Villaverde and Rubio-Ramírez (2018): +% "The Pruned State-Space System for Non-Linear DSGE Models: Theory and Empirical Applications", +% Review of Economic Studies, Volume 85, Issue 1, Pages 1–49. +% The model under study is taken from An and Schorfheide (2007): "Bayesian Analysis of DSGE Models", +% Econometric Reviews, Volume 26, Issue 2-4, Pages 113-172. +% Note that we use version 2 of the toolbox, i.e. the one which is called +% "Third-order GMM estimate package for DSGE models (version 2)" and can be +% downloaded from https://sites.google.com/site/mandreasendk/home-1 +% +% Created by @wmutschl (Willi Mutschler, willi@mutschler.eu) +% +% ========================================================================= +% Copyright (C) 2020 Dynare Team +% +% This file is part of Dynare. +% +% 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 . +% ========================================================================= + +% set this to 1 if you want to recompute using the Andreasen et al toolbox +% otherwise the results are loaded from Andreasen_et_al_2018_Dynare44Pruning_v2.mat +@#define Andreasen_et_al_toolbox = 0 + +var YGR INFL INT + c p R g y z; %if ordering of var is changed comparison code below needs to be adapted +varexo e_r e_g e_z; +parameters tau nu kap cyst psi1 psi2 rhor rhog rhoz rrst pist gamst; + +tau = 2; +nu = 0.1; +kap = 0.33; +cyst = 0.85; +psi1 = 1.5; +psi2 = 0.125; +rhor = 0.75; +rhog = 0.95; +rhoz = 0.9; +rrst = 1; +pist = 3.2; +gamst = 0.55; +sig_r = .2; +sig_g = .6; +sig_z = .3; + +model; +#pist2 = exp(pist/400); +#rrst2 = exp(rrst/400); +#bet = 1/rrst2; +#phi = tau*(1-nu)/nu/kap/pist2^2; +#gst = 1/cyst; +#cst = (1-nu)^(1/tau); +#yst = cst*gst; +#dy = y-y(-1); +1 = exp(-tau*c(+1)+tau*c+R-z(+1)-p(+1)); +(1-nu)/nu/phi/(pist2^2)*(exp(tau*c)-1) = (exp(p)-1)*((1-1/2/nu)*exp(p)+1/2/nu) - bet*(exp(p(+1))-1)*exp(-tau*c(+1)+tau*c+y(+1)-y+p(+1)); +exp(c-y) = exp(-g) - phi*pist2^2*gst/2*(exp(p)-1)^2; +R = rhor*R(-1) + (1-rhor)*psi1*p + (1-rhor)*psi2*(y-g) + e_r; +g = rhog*g(-1) + e_g; +z = rhoz*z(-1) + e_z; +YGR = gamst+100*(dy+z); +INFL = pist+400*p; +INT = pist+rrst+4*gamst+400*R; +end; + +shocks; +var e_r = sig_r^2; +var e_g = sig_g^2; +var e_z = sig_z^2; +end; + +steady_state_model; +y = 0; +R = 0; +g = 0; +z = 0; +c = 0; +p = 0; +YGR = gamst; +INFL = pist; +INT = pist + rrst + 4*gamst; +end; + +steady; check; model_diagnostics; + +@#for orderApp in [1, 2, 3] + stoch_simul(order=@{orderApp},pruning,irf=0,periods=0); + pruned_state_space.order_@{orderApp} = pruned_state_space_system(M_, options_, oo_.dr, [], options_.ar, 1, 0); + @#if Andreasen_et_al_toolbox + addpath('Dynare44Pruning_v2/simAndMoments3order'); %provide path to toolbox + optPruning.orderApp = @{orderApp}; + outAndreasenetal.order_@{orderApp} = RunDynarePruning(optPruning,oo_,M_,[oo_.dr.ghx oo_.dr.ghu]); + rmpath('Dynare44Pruning_v2/simAndMoments3order'); + close all; + @#endif +@#endfor + +@#if Andreasen_et_al_toolbox + delete Andreasen_et_al_2018_Dynare44Pruning_v2.mat; + pause(3); + save('Andreasen_et_al_2018_Dynare44Pruning_v2.mat', 'outAndreasenetal') + pause(3); +@#endif + +load('Andreasen_et_al_2018_Dynare44Pruning_v2.mat') + +% Make comparisons only at orders 1 and 2 +for iorder = 1:3 + fprintf('ORDER %d:\n',iorder); + pruned = pruned_state_space.(sprintf('order_%d',iorder)); + outAndreasen = outAndreasenetal.(sprintf('order_%d',iorder)); + %make sure variable ordering is correct + if ~isequal(M_.endo_names,[outAndreasen.label_y; outAndreasen.label_v(1:M_.nspred)]) + error('variable ordering is not the same, change declaration order'); + end + norm_E_yx = norm(pruned.E_y(oo_.dr.inv_order_var) - [outAndreasen.Mean_y; outAndreasen.Mean_v(1:M_.nspred)] , Inf); + fprintf('max(sum(abs(E[y;x]''))): %d\n',norm_E_yx); + norm_Var_y = norm(pruned.Var_y(oo_.dr.inv_order_var(1:(M_.endo_nbr-M_.nspred)),oo_.dr.inv_order_var(1:(M_.endo_nbr-M_.nspred))) - outAndreasen.Var_y , Inf); + fprintf('max(sum(abs(Var[y]''))):: %d\n',norm_Var_y); + norm_Var_x = norm(pruned.Var_y(M_.nstatic+(1:M_.nspred),M_.nstatic+(1:M_.nspred)) - outAndreasen.Var_v(1:M_.nspred,1:M_.nspred) , Inf); + fprintf('max(sum(abs(Var[x]''))): %d\n',norm_Var_x); + norm_Corr_yi1 = norm(pruned.Corr_yi(oo_.dr.inv_order_var(1:(M_.endo_nbr-M_.nspred)),oo_.dr.inv_order_var(1:(M_.endo_nbr-M_.nspred)),1) - outAndreasen.Corr_y(:,:,1) , Inf); + fprintf('max(sum(abs(Corr[y,y(-1)]''))): %d\n',norm_Corr_yi1); + norm_Corr_yi2 = norm(pruned.Corr_yi(oo_.dr.inv_order_var(1:(M_.endo_nbr-M_.nspred)),oo_.dr.inv_order_var(1:(M_.endo_nbr-M_.nspred)),2) - outAndreasen.Corr_y(:,:,2) , Inf); + fprintf('max(sum(abs(Corr[y,y(-2)]''))): %d\n',norm_Corr_yi2); + norm_Corr_xi1 = norm(pruned.Corr_yi(M_.nstatic+(1:M_.nspred),M_.nstatic+(1:M_.nspred),1) - outAndreasen.Corr_v(1:M_.nspred,1:M_.nspred,1) , Inf); + fprintf('max(sum(abs(Corr[x,x(-1)]''))): %d\n',norm_Corr_xi1); + norm_Corr_xi2 = norm(pruned.Corr_yi(M_.nstatic+(1:M_.nspred),M_.nstatic+(1:M_.nspred),2) - outAndreasen.Corr_v(1:M_.nspred,1:M_.nspred,2) , Inf); + fprintf('max(sum(abs(Corr[x,x(-2)]''))): %d\n',norm_Corr_xi2); + + if iorder < 3 && any([norm_E_yx norm_Var_y norm_Var_x norm_Corr_yi1 norm_Corr_yi2 norm_Corr_xi1 norm_Corr_xi2] > 1e-5) + error('Something wrong with pruned_state_space.m compared to Andreasen et al 2018 Toolbox v2 at order %d.',iorder); + end +end +skipline(); +fprintf('Note that at third order, there is an error in the computation of Var_z in Andreasen et al (2018)''s toolbox, @wmutschl is in contact to clarify this.\n'); +fprintf('EXAMPLE:\n') +fprintf(' Consider Var[kron(kron(xf,xf),xf)] = E[kron(kron(kron(kron(kron(xf,xf),xf),xf),xf),xf)] - E[kron(kron(xf,xf),xf)]*E[kron(kron(xf,xf),xf)].''\n'); +fprintf(' Now note that xf=hx*xf(-1)+hu*u is Gaussian, that is E[kron(kron(xf,xf),xf)]=0, and Var[kron(kron(xf,xf),xf)] are the sixth-order product moments\n'); +fprintf(' which can be computed using the prodmom.m function by providing E[xf*xf''] as covariance matrix.\n'); +fprintf(' In order to replicate this you have to change UnconditionalMoments_3rd_Lyap.m to also output Var_z.\n') + +dynare_nx = M_.nspred; +dynare_E_xf2 = pruned_state_space.order_3.Var_z(1:dynare_nx,1:dynare_nx); +dynare_E_xf6 = pruned_state_space.order_3.Var_z((end-dynare_nx^3+1):end,(end-dynare_nx^3+1):end); +dynare_E_xf6 = dynare_E_xf6(:); + +Andreasen_nx = M_.nspred+M_.exo_nbr; +Andreasen_E_xf2 = outAndreasenetal.order_3.Var_z(1:Andreasen_nx,1:Andreasen_nx); +Andreasen_E_xf6 = outAndreasenetal.order_3.Var_z((end-Andreasen_nx^3+1):end,(end-Andreasen_nx^3+1):end); +Andreasen_E_xf6 = Andreasen_E_xf6(:); + +fprintf('Second-order product moments of xf and u are the same:\n') +norm_E_xf2 = norm(dynare_E_xf2-Andreasen_E_xf2(1:M_.nspred,1:M_.nspred),Inf) +norm_E_uu = norm(M_.Sigma_e-Andreasen_E_xf2(M_.nspred+(1:M_.exo_nbr),M_.nspred+(1:M_.exo_nbr)),Inf) + +% Compute unique sixth-order product moments of xf, i.e. unique(E[kron(kron(kron(kron(kron(xf,xf),xf),xf),xf),xf)],'stable') +dynare_nx6 = dynare_nx*(dynare_nx+1)/2*(dynare_nx+2)/3*(dynare_nx+3)/4*(dynare_nx+4)/5*(dynare_nx+5)/6; +dynare_Q6Px = Q6_plication(dynare_nx); +dynare_COMBOS6 = flipud(allVL1(dynare_nx, 6)); %all possible (unique) combinations of powers that sum up to six +dynare_true_E_xf6 = zeros(dynare_nx6,1); %only unique entries +for j6 = 1:size(dynare_COMBOS6,1) + dynare_true_E_xf6(j6) = prodmom(dynare_E_xf2, 1:dynare_nx, dynare_COMBOS6(j6,:)); +end +dynare_true_E_xf6 = dynare_Q6Px*dynare_true_E_xf6; %add duplicate entries +norm_dynare_E_xf6 = norm(dynare_true_E_xf6 - dynare_E_xf6, Inf); + +Andreasen_nx6 = Andreasen_nx*(Andreasen_nx+1)/2*(Andreasen_nx+2)/3*(Andreasen_nx+3)/4*(Andreasen_nx+4)/5*(Andreasen_nx+5)/6; +Andreasen_Q6Px = Q6_plication(Andreasen_nx); +Andreasen_COMBOS6 = flipud(allVL1(Andreasen_nx, 6)); %all possible (unique) combinations of powers that sum up to six +Andreasen_true_E_xf6 = zeros(Andreasen_nx6,1); %only unique entries +for j6 = 1:size(Andreasen_COMBOS6,1) + Andreasen_true_E_xf6(j6) = prodmom(Andreasen_E_xf2, 1:Andreasen_nx, Andreasen_COMBOS6(j6,:)); +end +Andreasen_true_E_xf6 = Andreasen_Q6Px*Andreasen_true_E_xf6; %add duplicate entries +norm_Andreasen_E_xf6 = norm(Andreasen_true_E_xf6 - Andreasen_E_xf6, Inf); + +fprintf('Sixth-order product moments of xf and u are not the same!\n'); +fprintf(' Dynare maximum absolute deviations of sixth-order product moments of xf: %d\n',norm_dynare_E_xf6) +fprintf(' Andreasen et al maximum absolute deviations of sixth-order product moments of xf: %d\n',norm_Andreasen_E_xf6) +skipline(); +fprintf('Note that the standard deviations are set quite high to make the numerical differences more apparent.\n'); diff --git a/tests/pruning/Andreasen_et_al_2018_Dynare44Pruning_v2.mat b/tests/pruning/Andreasen_et_al_2018_Dynare44Pruning_v2.mat new file mode 100644 index 000000000..7bfbab582 Binary files /dev/null and b/tests/pruning/Andreasen_et_al_2018_Dynare44Pruning_v2.mat differ