diff --git a/matlab/Q6_plication.m b/matlab/Q6_plication.m new file mode 100644 index 000000000..c08e861f3 --- /dev/null +++ b/matlab/Q6_plication.m @@ -0,0 +1,64 @@ +% 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 +% +% Inputs: +% p: size of vector +% Outputs: +% QP: quadruplication matrix +% QPinv: Moore-Penrose inverse of QP +% + +function [DP6,DP6inv] = Q6_plication(p,progress) +if nargin <2 + progress =0; +end +reverseStr = ''; counti=1; +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); + +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); + n = ii1 + (ii2-1)*p + (ii3-1)*p^2 + (ii4-1)*p^3 + (ii5-1)*p^4 + (ii6-1)*p^5; + m = mue(p,i6,i5,i4,i3,i2,i1); + DP6(n,m)=1; + end + counti = counti+1; + end + end + end + end + end +end +DP6inv = (transpose(DP6)*DP6)\transpose(DP6); + +function m = mue(p,i1,i2,i3,i4,i5,i6) + 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) + t = q; r =p+q-i; + if t==0 + N=1; + else + N=1; + for h = 0:(t-1) + N = N*(r-h); + end + N=N/factorial(t); + end +end +end \ No newline at end of file diff --git a/matlab/allVL1.m b/matlab/allVL1.m new file mode 100644 index 000000000..3219ba1b7 --- /dev/null +++ b/matlab/allVL1.m @@ -0,0 +1,180 @@ +function v = allVL1(n, L1, L1ops, MaxNbSol) +% All integer permutations with sum criteria +% +% function v=allVL1(n, L1); OR +% v=allVL1(n, L1, L1opt); +% v=allVL1(n, L1, L1opt, MaxNbSol); +% +% INPUT +% n: length of the vector +% L1: target L1 norm +% L1ops: optional string ('==' or '<=' or '<') +% default value is '==' +% MaxNbSol: integer, returns at most MaxNbSol permutations. +% When MaxNbSol is NaN, allVL1 returns the total number of all possible +% permutations, which is useful to check the feasibility before getting +% the permutations. +% OUTPUT: +% v: (m x n) array such as: sum(v,2) == L1, +% (or <= or < depending on L1ops) +% all elements of v is naturel numbers {0,1,...} +% v contains all (=m) possible combinations +% v is sorted by sum (L1 norm), then by dictionnary sorting criteria +% class(v) is same as class(L1) +% Algorithm: +% Recursive +% Remark: +% allVL1(n,L1-n)+1 for natural numbers defined as {1,2,...} +% Example: +% This function can be used to generate all orders of all +% multivariable polynomials of degree p in R^n: +% Order = allVL1(n, p) +% Author: Bruno Luong +% Original, 30/nov/2007 +% Version 1.1, 30/apr/2008: Add H1 line as suggested by John D'Errico +% 1.2, 17/may/2009: Possibility to get the number of permutations +% alone (set fourth parameter MaxNbSol to NaN) +% 1.3, 16/Sep/2009: Correct bug for number of solution +% 1.4, 18/Dec/2010: + non-recursive engine + +global MaxCounter; + +if nargin<3 || isempty(L1ops) + L1ops = '=='; +end + +n = floor(n); % make sure n is integer + +if n<1 + v = []; + return +end + +if nargin<4 || isempty(MaxNbSol) + MaxCounter = Inf; +else + MaxCounter = MaxNbSol; +end +Counter(0); + +switch L1ops + case {'==' '='}, + if isnan(MaxCounter) + % return the number of solutions + v = nchoosek(n+L1-1,L1); % nchoosek(n+L1-1,n-1) + else + v = allVL1eq(n, L1); + end + case '<=', % call allVL1eq for various sum targets + if isnan(MaxCounter) + % return the number of solutions + %v = nchoosek(n+L1,L1)*factorial(n-L1); BUG <- 16/Sep/2009: + v = 0; + for j=0:L1 + v = v + nchoosek(n+j-1,j); + end + % See pascal's 11th identity, the sum doesn't seem to + % simplify to a fix formula + else + v = cell2mat(arrayfun(@(j) allVL1eq(n, j), (0:L1)', ... + 'UniformOutput', false)); + end + case '<', + v = allVL1(n, L1-1, '<=', MaxCounter); + otherwise + error('allVL1: unknown L1ops') +end + +end % allVL1 + +%% +function v = allVL1eq(n, L1) + +global MaxCounter; + +n = feval(class(L1),n); +s = n+L1; +sd = double(n)+double(L1); +notoverflowed = double(s)==sd; +if isinf(MaxCounter) && notoverflowed + v = allVL1nonrecurs(n, L1); +else + v = allVL1recurs(n, L1); +end + +end % allVL1eq + +%% Recursive engine +function v = allVL1recurs(n, L1, head) +% function v=allVL1eq(n, L1); +% INPUT +% n: length of the vector +% L1: desired L1 norm +% head: optional parameter to by concatenate in the first column +% of the output +% OUTPUT: +% if head is not defined +% v: (m x n) array such as sum(v,2)==L1 +% all elements of v is naturel numbers {0,1,...} +% v contains all (=m) possible combinations +% v is (dictionnary) sorted +% Algorithm: +% Recursive + +global MaxCounter; + +if n==1 + if Counter < MaxCounter + v = L1; + else + v = zeros(0,1,class(L1)); + end +else % recursive call + v = cell2mat(arrayfun(@(j) allVL1recurs(n-1, L1-j, j), (0:L1)', ... + 'UniformOutput', false)); +end + +if nargin>=3 % add a head column + v = [head+zeros(size(v,1),1,class(head)) v]; +end + +end % allVL1recurs + +%% +function res=Counter(newval) + persistent counter; + if nargin>=1 + counter = newval; + res = counter; + else + res = counter; + counter = counter+1; + end +end % Counter + +%% Non-recursive engine +function v = allVL1nonrecurs(n, L1) +% function v=allVL1eq(n, L1); +% INPUT +% n: length of the vector +% L1: desired L1 norm +% OUTPUT: +% if head is not defined +% v: (m x n) array such as sum(v,2)==L1 +% all elements of v is naturel numbers {0,1,...} +% v contains all (=m) possible combinations +% v is (dictionnary) sorted +% Algorithm: +% NonRecursive + +% Chose (n-1) the splitting points of the array [0:(n+L1)] +s = nchoosek(1:n+L1-1,n-1); +m = size(s,1); + +s1 = zeros(m,1,class(L1)); +s2 = (n+L1)+s1; + +v = diff([s1 s s2],1,2); % m x n +v = v-1; + +end % allVL1nonrecurs diff --git a/matlab/bivmom.m b/matlab/bivmom.m new file mode 100644 index 000000000..aa22eabdf --- /dev/null +++ b/matlab/bivmom.m @@ -0,0 +1,53 @@ +% +% 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. +% 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) +% +function [y,dy] = bivmom(p,rho) +s1 = p(1); +s2 = p(2); +rho2 = rho^2; +if nargout > 1 + drho2 = 2*rho; +end +if rem(s1+s2,2)==1 + y = 0; + return +end +r = fix(s1/2); +s = fix(s2/2); +y = 1; +c = 1; +if nargout > 1 + dy = 0; + dc = 0; +end +odd = 2*rem(s1,2); +for j=1:min(r,s) + if nargout > 1 + dc = 2*dc*(r+1-j)*(s+1-j)*rho2/(j*(2*j-1+odd)) + 2*c*(r+1-j)*(s+1-j)*drho2/(j*(2*j-1+odd)); + end + c = 2*c*(r+1-j)*(s+1-j)*rho2/(j*(2*j-1+odd)); + y = y+c; + if nargout > 1 + dy = dy + dc; + end +end +if odd + if nargout > 1 + dy = y + dy*rho; + end + y = y*rho; +end +y = prod([1:2:s1])*prod([1:2:s2])*y; +if nargout > 1 + dy = prod([1:2:s1])*prod([1:2:s2])*dy; +end + diff --git a/matlab/commutation.m b/matlab/commutation.m index 2dd3bc209..133ba97c9 100644 --- a/matlab/commutation.m +++ b/matlab/commutation.m @@ -13,7 +13,7 @@ function k = commutation(n, m, sparseflag) % k: [n by m] commutation matrix % ------------------------------------------------------------------------- % This function is called by -% * get_first_order_solution_params_deriv.m (previously getH.m) +% * get_perturbation_params_derivs.m (previously getH.m) % * get_identification_jacobians.m (previously getJJ.m) % ------------------------------------------------------------------------- % This function calls diff --git a/matlab/disp_identification.m b/matlab/disp_identification.m index bad140934..e7cfdc3f4 100644 --- a/matlab/disp_identification.m +++ b/matlab/disp_identification.m @@ -23,9 +23,6 @@ function disp_identification(pdraws, ide_reducedform, ide_moments, ide_spectrum, % ------------------------------------------------------------------------- % This function is called by % * dynare_identification.m -% ------------------------------------------------------------------------- -% This function calls -% * dynare_identification.m % ========================================================================= % Copyright (C) 2010-2019 Dynare Team % @@ -84,45 +81,56 @@ for jide = 1:4 no_warning_message_display = 1; %% Set output strings depending on test if jide == 1 - strTest = 'REDUCED-FORM'; strJacobian = 'Tau'; strMeaning = 'reduced-form solution'; + strTest = 'REDUCED-FORM'; strJacobian = 'Tau'; strMeaning = 'Jacobian of steady state and reduced-form solution matrices'; if ~no_identification_reducedform noidentification = 0; ide = ide_reducedform; if SampleSize == 1 - Jacob = ide.dTAU; + Jacob = ide.dREDUCEDFORM; end else %skip test noidentification = 1; no_warning_message_display = 0; end elseif jide == 2 - strTest = 'Iskrev (2010)'; strJacobian = 'J'; strMeaning = 'moments'; - if ~no_identification_moments - noidentification = 0; ide = ide_moments; + 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'; + elseif options_ident.order == 3 + strMeaning = 'Jacobian of second-order accurate mean and first-order minimal system'; + end + if ~no_identification_minimal + noidentification = 0; ide = ide_minimal; if SampleSize == 1 - Jacob = ide.si_J; + Jacob = ide.dMINIMAL; end else %skip test noidentification = 1; no_warning_message_display = 0; - end + end elseif jide == 3 - strTest = 'Komunjer and NG (2011)'; strJacobian = 'D'; strMeaning = 'minimal system'; - if ~no_identification_minimal - noidentification = 0; ide = ide_minimal; + strTest = 'SPECTRUM (Qu and Tkachenko, 2012)'; strJacobian = 'Gbar'; strMeaning = 'Jacobian of mean and spectrum'; + if options_ident.order > 1 + strTest = 'SPECTRUM (Mutschler, 2015)'; + end + if ~no_identification_spectrum + noidentification = 0; ide = ide_spectrum; if SampleSize == 1 - Jacob = ide.D; + Jacob = ide.dSPECTRUM; end else %skip test noidentification = 1; no_warning_message_display = 0; end elseif jide == 4 - strTest = 'Qu and Tkachenko (2012)'; strJacobian = 'G'; strMeaning = 'spectrum'; - if ~no_identification_spectrum - noidentification = 0; ide = ide_spectrum; + strTest = 'MOMENTS (Iskrev, 2010)'; strJacobian = 'J'; strMeaning = 'Jacobian of first two moments'; + if options_ident.order > 1 + strTest = 'MOMENTS (Mutschler, 2015)'; strJacobian = 'Mbar'; + end + if ~no_identification_moments + noidentification = 0; ide = ide_moments; if SampleSize == 1 - Jacob = ide.G; + Jacob = ide.si_dMOMENTS; end else %skip test noidentification = 1; no_warning_message_display = 0; - end + end end if ~noidentification @@ -176,7 +184,7 @@ for jide = 1:4 end end - %% display problematic parameters computed by identification_checks_via_subsets (only for debugging) + %% display problematic parameters computed by identification_checks_via_subsets elseif checks_via_subsets if ide.rank < size(Jacob,2) no_warning_message_display = 0; @@ -236,7 +244,7 @@ end %% Advanced identificaton patterns if SampleSize==1 && options_ident.advanced skipline() - for j=1:size(ide_moments.cosnJ,2) + for j=1:size(ide_moments.cosndMOMENTS,2) pax=NaN(totparam_nbr,totparam_nbr); fprintf('\n') disp(['Collinearity patterns with ', int2str(j) ,' parameter(s)']) @@ -249,10 +257,10 @@ if SampleSize==1 && options_ident.advanced namx=[namx ' ' sprintf('%-15s','--')]; else namx=[namx ' ' sprintf('%-15s',name{dumpindx})]; - pax(i,dumpindx)=ide_moments.cosnJ(i,j); + pax(i,dumpindx)=ide_moments.cosndMOMENTS(i,j); end end - fprintf('%-15s [%s] %14.7f\n',name{i},namx,ide_moments.cosnJ(i,j)) + fprintf('%-15s [%s] %14.7f\n',name{i},namx,ide_moments.cosndMOMENTS(i,j)) end end end diff --git a/matlab/dynare_estimation_init.m b/matlab/dynare_estimation_init.m index ab678c33b..729f5e43c 100644 --- a/matlab/dynare_estimation_init.m +++ b/matlab/dynare_estimation_init.m @@ -104,6 +104,10 @@ end if options_.order>2 && options_.particle.pruning error('Higher order nonlinear filters are not compatible with pruning option.') end +% Check the perturbation order (nonlinear filters with third order perturbation, or higher order, are not yet implemented). +if options_.order>2 && ~isfield(options_,'options_ident') %For identification at order=3 we skip this check. + error(['I cannot estimate a model with a ' int2str(options_.order) ' order approximation of the model!']) +end % analytical derivation is not yet available for kalman_filter_fast if options_.analytic_derivation && options_.fast_kalman_filter diff --git a/matlab/dynare_identification.m b/matlab/dynare_identification.m index 6a09e8705..62164b213 100644 --- a/matlab/dynare_identification.m +++ b/matlab/dynare_identification.m @@ -1,32 +1,32 @@ -function [pdraws, STO_TAU, STO_MOMENTS, STO_LRE, STO_si_dLRE, STO_si_dTAU, STO_si_J, STO_G, STO_D] = dynare_identification(options_ident, pdraws0) -%function [pdraws, STO_TAU, STO_MOMENTS, STO_LRE, STO_dLRE, STO_dTAU, STO_J, STO_G, STO_D] = dynare_identification(options_ident, pdraws0) +function [pdraws, STO_REDUCEDFORM, STO_MOMENTS, STO_DYNAMIC, STO_si_dDYNAMIC, STO_si_dREDUCEDFORM, STO_si_dMOMENTS, STO_dSPECTRUM, STO_dMINIMAL] = dynare_identification(options_ident, pdraws0) +%function [pdraws, STO_REDUCEDFORM, STO_MOMENTS, STO_DYNAMIC, STO_si_dDYNAMIC, STO_si_dREDUCEDFORM, STO_si_dMOMENTS, STO_dSPECTRUM, STO_dMINIMAL] = dynare_identification(options_ident, pdraws0) % ------------------------------------------------------------------------- % This function is called, when the user specifies identification(...); in % the mod file. It prepares all identification analysis, i.e. -% (1) sets options, local/persistent/global variables for a new identification -% analysis either for a single point or MC Sample and displays and plots the results -% or -% (2) loads, displays and plots a previously saved identification analysis +% (1) sets options, local and persistent variables for a new identification +% analysis either for a single point or a MC Sample. It also displays +% and plots the results. +% or this function can be used to +% (2) load, display and plot a previously saved identification analysis +% Note: This function does not output the arguments to the workspace if only called by +% "identification" in the mod file, but saves results to the folder identification. +% If you want to use this function directly in the mod file and workspace, you still have +% to put identification once in the mod file, otherwise the preprocessor won't compute all necessary objects % ========================================================================= % INPUTS % * options_ident [structure] identification options % * pdraws0 [SampleSize by totparam_nbr] optional: matrix of MC sample of model parameters % ------------------------------------------------------------------------- % OUTPUTS -% Note: This function does not output the arguments to the workspace if only called by -% "identification" in the mod file, but saves results to the folder identification. -% One can, however, just use -% [pdraws, STO_TAU, STO_MOMENTS, STO_LRE, STO_dLRE, STO_dTAU, STO_J, STO_G, STO_D] = dynare_identification(options_ident, pdraws0) -% in the mod file to get the results directly in the workspace -% * pdraws [matrix] MC sample of model params used -% * STO_TAU, [matrix] MC sample of entries in the model solution (stacked vertically) -% * STO_MOMENTS, [matrix] MC sample of entries in the moments (stacked vertically) -% * STO_LRE, [matrix] MC sample of entries in LRE model (stacked vertically) -% * STO_dLRE, [matrix] MC sample of derivatives of the Jacobian (dLRE) -% * STO_dTAU, [matrix] MC sample of derivatives of the model solution and steady state (dTAU) -% * STO_J [matrix] MC sample of Iskrev (2010)'s J matrix -% * STO_G [matrix] MC sample of Qu and Tkachenko (2012)'s G matrix -% * STO_D [matrix] MC sample of Komunjer and Ng (2011)'s D matrix +% * pdraws [matrix] MC sample of model params used +% * STO_REDUCEDFORM, [matrix] MC sample of entries in steady state and reduced form model solution (stacked vertically) +% * STO_MOMENTS, [matrix] MC sample of entries in theoretical first two moments (stacked vertically) +% * STO_DYNAMIC, [matrix] MC sample of entries in steady state and dynamic model derivatives (stacked vertically) +% * STO_si_dDYNAMIC, [matrix] MC sample of derivatives of steady state and dynamic derivatives +% * STO_si_dREDUCEDFORM, [matrix] MC sample of derivatives of steady state and reduced form model solution +% * STO_si_dMOMENTS [matrix] MC sample of Iskrev (2010)'s J matrix +% * STO_dSPECTRUM [matrix] MC sample of Qu and Tkachenko (2012)'s \bar{G} matrix +% * STO_dMINIMAL [matrix] MC sample of Komunjer and Ng (2011)'s \bar{\Delta} matrix % ------------------------------------------------------------------------- % This function is called by % * driver.m @@ -107,10 +107,10 @@ options_ident = set_default_option(options_ident,'parameter_set','prior_mean'); options_ident = set_default_option(options_ident,'load_ident_files',0); % 1: load previously computed analysis from identification/fname_identif.mat options_ident = set_default_option(options_ident,'useautocorr',0); - % 1: use autocorrelations in Iskrev (2010)'s J criteria - % 0: use autocovariances in Iskrev (2010)'s J criteria + % 1: use autocorrelations in Iskrev (2010)'s criteria + % 0: use autocovariances in Iskrev (2010)'s criteria options_ident = set_default_option(options_ident,'ar',1); - % number of lags to consider for autocovariances/autocorrelations in Iskrev (2010)'s J criteria + % number of lags to consider for autocovariances/autocorrelations in Iskrev (2010)'s criteria options_ident = set_default_option(options_ident,'prior_mc',1); % size of Monte-Carlo sample of parameter draws options_ident = set_default_option(options_ident,'prior_range',0); @@ -121,16 +121,17 @@ options_ident = set_default_option(options_ident,'periods',300); options_ident = set_default_option(options_ident,'replic',100); % number of replicas to compute simulated moment uncertainty, when analytic Hessian is not available options_ident = set_default_option(options_ident,'advanced',0); - % 1: show a more detailed analysis based on reduced-form solution and Jacobian of dynamic model (LRE). Further, performs a brute force - % search of the groups of parameters best reproducing the behavior of each single parameter of Iskrev (2010)'s J. + % 1: show a more detailed analysis based on reduced-form model solution and dynamic model derivatives. Further, performs a brute force + % search of the groups of parameters best reproducing the behavior of each single parameter. options_ident = set_default_option(options_ident,'normalize_jacobians',1); - % 1: normalize Jacobians by rescaling each row by its largest element in absolute value + % 1: normalize Jacobians by either rescaling each row by its largest element in absolute value or for Gram (or Hessian-type) matrices by transforming into correlation-type matrices options_ident = set_default_option(options_ident,'grid_nbr',5000); - % number of grid points in [-pi;pi] to approximate the integral to compute Qu and Tkachenko (2012)'s G criteria + % number of grid points in [-pi;pi] to approximate the integral to compute Qu and Tkachenko (2012)'s criteria % note that grid_nbr needs to be even and actually we use (grid_nbr+1) points, as we add the 0 frequency and use symmetry, i.e. grid_nbr/2 % negative as well as grid_nbr/2 positive values to speed up the compuations if mod(options_ident.grid_nbr,2) ~= 0 options_ident.grid_nbr = options_ident.grid_nbr+1; + fprintf('IDENTIFICATION: ''grid_nbr'' needs to be even, hence it is reset to %d\n',options_ident.grid_nbr) if mod(options_ident.grid_nbr,2) ~= 0 error('IDENTIFICATION: You need to set an even value for ''grid_nbr'''); end @@ -148,25 +149,25 @@ if ~isfield(options_ident,'no_identification_strength') else options_ident.no_identification_strength = 1; end -%check whether to compute and display identification criteria based on steady state and reduced-form solution +%check whether to compute and display identification criteria based on steady state and reduced-form model solution if ~isfield(options_ident,'no_identification_reducedform') options_ident.no_identification_reducedform = 0; else options_ident.no_identification_reducedform = 1; end -%check whether to compute and display identification criteria based on Iskrev (2010)'s J, i.e. derivative of first two moments +%check whether to compute and display identification criteria based on Iskrev (2010), i.e. derivative of first two moments if ~isfield(options_ident,'no_identification_moments') options_ident.no_identification_moments = 0; else options_ident.no_identification_moments = 1; end -%check whether to compute and display identification criteria based on Komunjer and Ng (2011)'s D, i.e. derivative of first moment, minimal state space system and observational equivalent spectral density transformation +%check whether to compute and display identification criteria based on Komunjer and Ng (2011), i.e. derivative of first moment, minimal state space system and observational equivalent spectral density transformation if ~isfield(options_ident,'no_identification_minimal') options_ident.no_identification_minimal = 0; else - options_ident.no_identification_minimal = 1; + options_ident.no_identification_minimal = 1; end -%Check whether to compute and display identification criteria based on Qu and Tkachenko (2012)'s G, i.e. Gram matrix of derivatives of first moment plus outer product of derivatives of spectral density +%Check whether to compute and display identification criteria based on Qu and Tkachenko (2012), i.e. Gram matrix of derivatives of first moment plus outer product of derivatives of spectral density if ~isfield(options_ident,'no_identification_spectrum') options_ident.no_identification_spectrum = 0; else @@ -211,6 +212,10 @@ options_ident = set_default_option(options_ident,'lik_init',1); % ii) option 1 for the stationary elements options_ident = set_default_option(options_ident,'analytic_derivation',1); % 1: analytic derivation of gradient and hessian of likelihood in dsge_likelihood.m, only works for stationary models, i.e. kalman_algo<3 +options_ident = set_default_option(options_ident,'order',1); + % 1: first-order perturbation approximation, identification is based on linear state space system + % 2: second-order perturbation approximation, identification is based on second-order pruned state space system + % 3: third-order perturbation approximation, identification is based on third-order pruned state space system %overwrite values in options_, note that options_ is restored at the end of the function if isfield(options_ident,'TeX') @@ -273,7 +278,17 @@ else end % overwrite settings in options_ and prepare to call dynare_estimation_init -options_.order = 1; % Identification does not support order>1 (yet) +options_.order = options_ident.order; +if options_ident.order > 1 + %order>1 is not compatible with analytic derivation in dsge_likelihood.m + options_ident.analytic_derivation=0; + options_.analytic_derivation=0; + %order>1 is based on pruned state space system + options_.pruning = true; +end +if options_ident.order == 3 + options_.k_order_solver = 1; +end options_.ar = options_ident.ar; options_.prior_mc = options_ident.prior_mc; options_.Schur_vec_tol = 1.e-8; @@ -284,15 +299,15 @@ options_ = set_default_option(options_,'datafile',''); options_.mode_compute = 0; options_.plot_priors = 0; options_.smoother = 1; +options_.options_ident = []; [dataset_, dataset_info, xparam1, hh, M_, options_, oo_, estim_params_, bayestopt_, bounds] = dynare_estimation_init(M_.endo_names, fname, 1, M_, options_, oo_, estim_params_, bayestopt_); %outputting dataset_ is needed for Octave - % set method to compute identification Jacobians (kronflag). Default:0 options_ident = set_default_option(options_ident,'analytic_derivation_mode', options_.analytic_derivation_mode); % if not set by user, inherit default global one - % 0: efficient sylvester equation method to compute analytical derivatives as in Ratto & Iskrev (2011) - % 1: kronecker products method to compute analytical derivatives as in Iskrev (2010) - % -1: numerical two-sided finite difference method to compute numerical derivatives of all Jacobians using function identification_numerical_objective.m (previously thet2tau.m) - % -2: numerical two-sided finite difference method to compute numerically dYss, dg1, d2Yss and d2g1, the Jacobians are then computed analytically as in options 0 + % 0: efficient sylvester equation method to compute analytical derivatives as in Ratto & Iskrev (2012) + % 1: kronecker products method to compute analytical derivatives as in Iskrev (2010) (only for order=1) + % -1: numerical two-sided finite difference method to compute numerical derivatives of all identification Jacobians using function identification_numerical_objective.m (previously thet2tau.m) + % -2: numerical two-sided finite difference method to compute numerically dYss, dg1, dg2, dg3, d2Yss and d2g1, the identification Jacobians are then computed analytically as if option is set to 0 options_.analytic_derivation_mode = options_ident.analytic_derivation_mode; %overwrite setting % initialize persistent variables in prior_draw @@ -365,7 +380,10 @@ end options_ident.name_tex = name_tex; skipline() -disp('==== Identification analysis ====') +fprintf('======== Identification Analysis ========\n') +if options_ident.order > 1 + fprintf('Using Pruned State Space System (order=%d)\n',options_ident.order); +end skipline() if totparam_nbr < 2 options_ident.advanced = 0; @@ -378,20 +396,19 @@ options_ident = set_default_option(options_ident,'max_dim_cova_group',min([2,tot options_ident.max_dim_cova_group = min([options_ident.max_dim_cova_group,totparam_nbr-1]); % In brute force search (ident_bruteforce.m) when advanced=1 this option sets the maximum dimension of groups of parameters that best reproduce the behavior of each single model parameter -options_ident = set_default_option(options_ident,'checks_via_subsets',0); %[ONLY FOR DEBUGGING] +options_ident = set_default_option(options_ident,'checks_via_subsets',0); % 1: uses identification_checks_via_subsets.m to compute problematic parameter combinations % 0: uses identification_checks.m to compute problematic parameter combinations [default] -options_ident = set_default_option(options_ident,'max_dim_subsets_groups',min([4,totparam_nbr-1])); %[ONLY FOR DEBUGGING] - % In identification_checks_via_subsets.m, when checks_via_subsets=1, this - % option sets the maximum dimension of groups of parameters for which - % the corresponding rank criteria is checked +options_ident = set_default_option(options_ident,'max_dim_subsets_groups',min([4,totparam_nbr-1])); + % In identification_checks_via_subsets.m, when checks_via_subsets=1, this option sets the maximum dimension of groups of parameters for which the corresponding rank criteria is checked options_.options_ident = options_ident; %store identification options into global options -MaxNumberOfBytes = options_.MaxNumberOfBytes; %threshold when to save from memory to files store_options_ident = options_ident; +MaxNumberOfBytes = options_.MaxNumberOfBytes; %threshold when to save from memory to files iload = options_ident.load_ident_files; SampleSize = options_ident.prior_mc; + if iload <=0 %% Perform new identification analysis, i.e. do not load previous analysis if prior_exist @@ -442,14 +459,14 @@ if iload <=0 end else % no estimated_params block is available, all stderr and model parameters, but no corr parameters are chosen - params = [sqrt(diag(M_.Sigma_e))', M_.params']; % use current values + params = [sqrt(diag(M_.Sigma_e))', M_.params']; parameters = 'Current_params'; parameters_TeX = 'Current parameter values'; disp('Testing all current stderr and model parameter values') end options_ident.tittxt = parameters; %title text for graphs and figures % perform identification analysis for single point - [ide_moments_point, ide_spectrum_point, ide_minimal_point, ide_hess_point, ide_reducedform_point, ide_lre_point, derivatives_info_point, info, options_ident] = ... + [ide_moments_point, ide_spectrum_point, ide_minimal_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, derivatives_info_point, info, options_ident] = ... identification_analysis(params, indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1); %the 1 at the end implies initialization of persistent variables if info(1)~=0 % there are errors in the solution algorithm @@ -496,7 +513,7 @@ if iload <=0 params = prior_draw(); options_ident.tittxt = 'Random_prior_params'; %title text for graphs and figures % perform identification analysis - [ide_moments_point, ide_spectrum_point, ide_minimal_point, ide_hess_point, ide_reducedform_point, ide_lre_point, derivatives_info, info, options_ident] = ... + [ide_moments_point, ide_spectrum_point, ide_minimal_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, derivatives_info_point, info, options_ident] = ... identification_analysis(params,indpmodel,indpstderr,indpcorr,options_ident,dataset_info, prior_exist, 1); end end @@ -521,13 +538,13 @@ if iload <=0 end ide_hess_point.params = params; % save all output into identification folder - save([IdentifDirectoryName '/' fname '_identif.mat'], 'ide_moments_point', 'ide_spectrum_point', 'ide_minimal_point', 'ide_hess_point', 'ide_reducedform_point', 'ide_lre_point','store_options_ident'); - save([IdentifDirectoryName '/' fname '_' parameters '_identif.mat'], 'ide_moments_point', 'ide_spectrum_point', 'ide_minimal_point', 'ide_hess_point', 'ide_reducedform_point', 'ide_lre_point','store_options_ident'); + save([IdentifDirectoryName '/' fname '_identif.mat'], 'ide_moments_point', 'ide_spectrum_point', 'ide_minimal_point', 'ide_hess_point', 'ide_reducedform_point', 'ide_dynamic_point','store_options_ident'); + save([IdentifDirectoryName '/' fname '_' parameters '_identif.mat'], 'ide_moments_point', 'ide_spectrum_point', 'ide_minimal_point', 'ide_hess_point', 'ide_reducedform_point', 'ide_dynamic_point','store_options_ident'); % display results of identification analysis disp_identification(params, ide_reducedform_point, ide_moments_point, ide_spectrum_point, ide_minimal_point, name, options_ident); if ~options_ident.no_identification_strength && ~options_.nograph % plot (i) identification strength and sensitivity measure based on the sample information matrix and (ii) advanced analysis graphs - plot_identification(params, ide_moments_point, ide_hess_point, ide_reducedform_point, ide_lre_point, options_ident.advanced, parameters, name, IdentifDirectoryName, parameters_TeX, name_tex); + plot_identification(params, ide_moments_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, options_ident.advanced, parameters, name, IdentifDirectoryName, parameters_TeX, name_tex); end if SampleSize > 1 @@ -553,30 +570,30 @@ if iload <=0 end options_ident.tittxt = []; % clear title text for graphs and figures % run identification analysis - [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide_lre, ide_derivatives_info, info, options_MC] = ... + [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide_dynamic, ide_derivatives_info, info, options_MC] = ... identification_analysis(params, indpmodel, indpstderr, indpcorr, options_MC, dataset_info, prior_exist, 0); % the 0 implies that we do not initialize persistent variables anymore - + if iteration==0 && info(1)==0 % preallocate storage in the first admissable run delete([IdentifDirectoryName '/' fname '_identif_*.mat']) % delete previously saved results - MAX_RUNS_BEFORE_SAVE_TO_FILE = min(SampleSize,ceil(MaxNumberOfBytes/(size(ide_reducedform.si_dTAU,1)*totparam_nbr)/8)); % set how many runs can be stored before we save to files + MAX_RUNS_BEFORE_SAVE_TO_FILE = min(SampleSize,ceil(MaxNumberOfBytes/(size(ide_reducedform.si_dREDUCEDFORM,1)*totparam_nbr)/8)); % set how many runs can be stored before we save to files pdraws = zeros(SampleSize,totparam_nbr); % preallocate storage for draws in each row - % preallocate storage for linear rational expectations model - STO_si_dLRE = zeros([size(ide_lre.si_dLRE,1),modparam_nbr,MAX_RUNS_BEFORE_SAVE_TO_FILE]); - STO_LRE = zeros(size(ide_lre.LRE,1),SampleSize); - IDE_LRE.ind_dLRE = ide_lre.ind_dLRE; - IDE_LRE.in0 = zeros(SampleSize,modparam_nbr); - IDE_LRE.ind0 = zeros(SampleSize,modparam_nbr); - IDE_LRE.jweak = zeros(SampleSize,modparam_nbr); - IDE_LRE.jweak_pair = zeros(SampleSize,modparam_nbr*(modparam_nbr+1)/2); - IDE_LRE.cond = zeros(SampleSize,1); - IDE_LRE.Mco = zeros(SampleSize,modparam_nbr); + % preallocate storage for steady state and dynamic model derivatives + STO_si_dDYNAMIC = zeros([size(ide_dynamic.si_dDYNAMIC,1),modparam_nbr,MAX_RUNS_BEFORE_SAVE_TO_FILE]); + STO_DYNAMIC = zeros(size(ide_dynamic.DYNAMIC,1),SampleSize); + IDE_DYNAMIC.ind_dDYNAMIC = ide_dynamic.ind_dDYNAMIC; + IDE_DYNAMIC.in0 = zeros(SampleSize,modparam_nbr); + IDE_DYNAMIC.ind0 = zeros(SampleSize,modparam_nbr); + IDE_DYNAMIC.jweak = zeros(SampleSize,modparam_nbr); + IDE_DYNAMIC.jweak_pair = zeros(SampleSize,modparam_nbr*(modparam_nbr+1)/2); + IDE_DYNAMIC.cond = zeros(SampleSize,1); + IDE_DYNAMIC.Mco = zeros(SampleSize,modparam_nbr); % preallocate storage for reduced form if ~options_MC.no_identification_reducedform - STO_si_dTAU = zeros([size(ide_reducedform.si_dTAU,1),totparam_nbr,MAX_RUNS_BEFORE_SAVE_TO_FILE]); - STO_TAU = zeros(size(ide_reducedform.TAU,1),SampleSize); - IDE_REDUCEDFORM.ind_dTAU = ide_reducedform.ind_dTAU; + STO_si_dREDUCEDFORM = zeros([size(ide_reducedform.si_dREDUCEDFORM,1),totparam_nbr,MAX_RUNS_BEFORE_SAVE_TO_FILE]); + STO_REDUCEDFORM = zeros(size(ide_reducedform.REDUCEDFORM,1),SampleSize); + IDE_REDUCEDFORM.ind_dREDUCEDFORM = ide_reducedform.ind_dREDUCEDFORM; IDE_REDUCEDFORM.in0 = zeros(SampleSize,1); IDE_REDUCEDFORM.ind0 = zeros(SampleSize,totparam_nbr); IDE_REDUCEDFORM.jweak = zeros(SampleSize,totparam_nbr); @@ -589,45 +606,45 @@ if iload <=0 % preallocate storage for moments if ~options_MC.no_identification_moments - STO_si_J = zeros([size(ide_moments.si_J,1),totparam_nbr,MAX_RUNS_BEFORE_SAVE_TO_FILE]); - STO_MOMENTS = zeros(size(ide_moments.MOMENTS,1),SampleSize); - IDE_MOMENTS.ind_J = ide_moments.ind_J; - IDE_MOMENTS.in0 = zeros(SampleSize,1); - IDE_MOMENTS.ind0 = zeros(SampleSize,totparam_nbr); - IDE_MOMENTS.jweak = zeros(SampleSize,totparam_nbr); - IDE_MOMENTS.jweak_pair = zeros(SampleSize,totparam_nbr*(totparam_nbr+1)/2); - IDE_MOMENTS.cond = zeros(SampleSize,1); - IDE_MOMENTS.Mco = zeros(SampleSize,totparam_nbr); - IDE_MOMENTS.S = zeros(SampleSize,min(8,totparam_nbr)); - IDE_MOMENTS.V = zeros(SampleSize,totparam_nbr,min(8,totparam_nbr)); + STO_si_dMOMENTS = zeros([size(ide_moments.si_dMOMENTS,1),totparam_nbr,MAX_RUNS_BEFORE_SAVE_TO_FILE]); + STO_MOMENTS = zeros(size(ide_moments.MOMENTS,1),SampleSize); + IDE_MOMENTS.ind_dMOMENTS = ide_moments.ind_dMOMENTS; + IDE_MOMENTS.in0 = zeros(SampleSize,1); + IDE_MOMENTS.ind0 = zeros(SampleSize,totparam_nbr); + IDE_MOMENTS.jweak = zeros(SampleSize,totparam_nbr); + IDE_MOMENTS.jweak_pair = zeros(SampleSize,totparam_nbr*(totparam_nbr+1)/2); + IDE_MOMENTS.cond = zeros(SampleSize,1); + IDE_MOMENTS.Mco = zeros(SampleSize,totparam_nbr); + IDE_MOMENTS.S = zeros(SampleSize,min(8,totparam_nbr)); + IDE_MOMENTS.V = zeros(SampleSize,totparam_nbr,min(8,totparam_nbr)); else IDE_MOMENTS = {}; end % preallocate storage for spectrum if ~options_MC.no_identification_spectrum - STO_G = zeros([size(ide_spectrum.G,1),size(ide_spectrum.G,2), MAX_RUNS_BEFORE_SAVE_TO_FILE]); - IDE_SPECTRUM.ind_G = ide_spectrum.ind_G; - IDE_SPECTRUM.in0 = zeros(SampleSize,1); - IDE_SPECTRUM.ind0 = zeros(SampleSize,totparam_nbr); - IDE_SPECTRUM.jweak = zeros(SampleSize,totparam_nbr); - IDE_SPECTRUM.jweak_pair = zeros(SampleSize,totparam_nbr*(totparam_nbr+1)/2); - IDE_SPECTRUM.cond = zeros(SampleSize,1); - IDE_SPECTRUM.Mco = zeros(SampleSize,totparam_nbr); + STO_dSPECTRUM = zeros([size(ide_spectrum.dSPECTRUM,1),size(ide_spectrum.dSPECTRUM,2), MAX_RUNS_BEFORE_SAVE_TO_FILE]); + IDE_SPECTRUM.ind_dSPECTRUM = ide_spectrum.ind_dSPECTRUM; + IDE_SPECTRUM.in0 = zeros(SampleSize,1); + IDE_SPECTRUM.ind0 = zeros(SampleSize,totparam_nbr); + IDE_SPECTRUM.jweak = zeros(SampleSize,totparam_nbr); + IDE_SPECTRUM.jweak_pair = zeros(SampleSize,totparam_nbr*(totparam_nbr+1)/2); + IDE_SPECTRUM.cond = zeros(SampleSize,1); + IDE_SPECTRUM.Mco = zeros(SampleSize,totparam_nbr); else IDE_SPECTRUM = {}; end % preallocate storage for minimal system if ~options_MC.no_identification_minimal - STO_D = zeros([size(ide_minimal.D,1),size(ide_minimal.D,2), MAX_RUNS_BEFORE_SAVE_TO_FILE]); - IDE_MINIMAL.ind_D = ide_minimal.ind_D; - IDE_MINIMAL.in0 = zeros(SampleSize,1); - IDE_MINIMAL.ind0 = zeros(SampleSize,totparam_nbr); - IDE_MINIMAL.jweak = zeros(SampleSize,totparam_nbr); - IDE_MINIMAL.jweak_pair = zeros(SampleSize,totparam_nbr*(totparam_nbr+1)/2); - IDE_MINIMAL.cond = zeros(SampleSize,1); - IDE_MINIMAL.Mco = zeros(SampleSize,totparam_nbr); + STO_dMINIMAL = zeros([size(ide_minimal.dMINIMAL,1),size(ide_minimal.dMINIMAL,2), MAX_RUNS_BEFORE_SAVE_TO_FILE]); + IDE_MINIMAL.ind_dMINIMAL = ide_minimal.ind_dMINIMAL; + IDE_MINIMAL.in0 = zeros(SampleSize,1); + IDE_MINIMAL.ind0 = zeros(SampleSize,totparam_nbr); + IDE_MINIMAL.jweak = zeros(SampleSize,totparam_nbr); + IDE_MINIMAL.jweak_pair = zeros(SampleSize,totparam_nbr*(totparam_nbr+1)/2); + IDE_MINIMAL.cond = zeros(SampleSize,1); + IDE_MINIMAL.Mco = zeros(SampleSize,totparam_nbr); else IDE_MINIMAL = {}; end @@ -638,20 +655,20 @@ if iload <=0 run_index = run_index + 1; %increase index of admissable draws after saving to files pdraws(iteration,:) = params; % store draw - % store results for linear rational expectations model - STO_LRE(:,iteration) = ide_lre.LRE; - STO_si_dLRE(:,:,run_index) = ide_lre.si_dLRE; - IDE_LRE.cond(iteration,1) = ide_lre.cond; - IDE_LRE.ino(iteration,1) = ide_lre.ino; - IDE_LRE.ind0(iteration,:) = ide_lre.ind0; - IDE_LRE.jweak(iteration,:) = ide_lre.jweak; - IDE_LRE.jweak_pair(iteration,:) = ide_lre.jweak_pair; - IDE_LRE.Mco(iteration,:) = ide_lre.Mco; + % store results for steady state and dynamic model derivatives + STO_DYNAMIC(:,iteration) = ide_dynamic.DYNAMIC; + STO_si_dDYNAMIC(:,:,run_index) = ide_dynamic.si_dDYNAMIC; + IDE_DYNAMIC.cond(iteration,1) = ide_dynamic.cond; + IDE_DYNAMIC.ino(iteration,1) = ide_dynamic.ino; + IDE_DYNAMIC.ind0(iteration,:) = ide_dynamic.ind0; + IDE_DYNAMIC.jweak(iteration,:) = ide_dynamic.jweak; + IDE_DYNAMIC.jweak_pair(iteration,:) = ide_dynamic.jweak_pair; + IDE_DYNAMIC.Mco(iteration,:) = ide_dynamic.Mco; - % store results for reduced form solution + % store results for reduced form model solution if ~options_MC.no_identification_reducedform - STO_TAU(:,iteration) = ide_reducedform.TAU; - STO_si_dTAU(:,:,run_index) = ide_reducedform.si_dTAU; + STO_REDUCEDFORM(:,iteration) = ide_reducedform.REDUCEDFORM; + STO_si_dREDUCEDFORM(:,:,run_index) = ide_reducedform.si_dREDUCEDFORM; IDE_REDUCEDFORM.cond(iteration,1) = ide_reducedform.cond; IDE_REDUCEDFORM.ino(iteration,1) = ide_reducedform.ino; IDE_REDUCEDFORM.ind0(iteration,:) = ide_reducedform.ind0; @@ -663,7 +680,7 @@ if iload <=0 % store results for moments if ~options_MC.no_identification_moments STO_MOMENTS(:,iteration) = ide_moments.MOMENTS; - STO_si_J(:,:,run_index) = ide_moments.si_J; + STO_si_dMOMENTS(:,:,run_index) = ide_moments.si_dMOMENTS; IDE_MOMENTS.cond(iteration,1) = ide_moments.cond; IDE_MOMENTS.ino(iteration,1) = ide_moments.ino; IDE_MOMENTS.ind0(iteration,:) = ide_moments.ind0; @@ -676,7 +693,7 @@ if iload <=0 % store results for spectrum if ~options_MC.no_identification_spectrum - STO_G(:,:,run_index) = ide_spectrum.G; + STO_dSPECTRUM(:,:,run_index) = ide_spectrum.dSPECTRUM; IDE_SPECTRUM.cond(iteration,1) = ide_spectrum.cond; IDE_SPECTRUM.ino(iteration,1) = ide_spectrum.ino; IDE_SPECTRUM.ind0(iteration,:) = ide_spectrum.ind0; @@ -687,7 +704,7 @@ if iload <=0 % store results for minimal system if ~options_MC.no_identification_minimal - STO_D(:,:,run_index) = ide_minimal.D; + STO_dMINIMAL(:,:,run_index) = ide_minimal.dMINIMAL; IDE_MINIMAL.cond(iteration,1) = ide_minimal.cond; IDE_MINIMAL.ino(iteration,1) = ide_minimal.ino; IDE_MINIMAL.ind0(iteration,:) = ide_minimal.ind0; @@ -701,37 +718,37 @@ if iload <=0 file_index = file_index + 1; if run_index 1 dyn_waitbar_close(h); - normalize_STO_LRE = std(STO_LRE,0,2); + normalize_STO_DYNAMIC = std(STO_DYNAMIC,0,2); if ~options_MC.no_identification_reducedform - normalize_STO_TAU = std(STO_TAU,0,2); + normalize_STO_REDUCEDFORM = std(STO_REDUCEDFORM,0,2); end if ~options_MC.no_identification_moments normalize_STO_MOMENTS = std(STO_MOMENTS,0,2); @@ -759,65 +776,65 @@ if iload <=0 normaliz1 = std(pdraws); iter = 0; for ifile_index = 1:file_index - load([IdentifDirectoryName '/' fname '_identif_' int2str(ifile_index) '.mat'], 'STO_si_dLRE'); - maxrun_dLRE = size(STO_si_dLRE,3); + load([IdentifDirectoryName '/' fname '_identif_' int2str(ifile_index) '.mat'], 'STO_si_dDYNAMIC'); + maxrun_dDYNAMIC = size(STO_si_dDYNAMIC,3); if ~options_MC.no_identification_reducedform - load([IdentifDirectoryName '/' fname '_identif_' int2str(ifile_index) '.mat'], 'STO_si_dTAU'); - maxrun_dTAU = size(STO_si_dTAU,3); + load([IdentifDirectoryName '/' fname '_identif_' int2str(ifile_index) '.mat'], 'STO_si_dREDUCEDFORM'); + maxrun_dREDUCEDFORM = size(STO_si_dREDUCEDFORM,3); else - maxrun_dTAU = 0; + maxrun_dREDUCEDFORM = 0; end if ~options_MC.no_identification_moments - load([IdentifDirectoryName '/' fname '_identif_' int2str(ifile_index) '.mat'], 'STO_si_J'); - maxrun_J = size(STO_si_J,3); + load([IdentifDirectoryName '/' fname '_identif_' int2str(ifile_index) '.mat'], 'STO_si_dMOMENTS'); + maxrun_dMOMENTS = size(STO_si_dMOMENTS,3); else - maxrun_J = 0; + maxrun_dMOMENTS = 0; end if ~options_MC.no_identification_spectrum - load([IdentifDirectoryName '/' fname '_identif_' int2str(ifile_index) '.mat'], 'STO_G'); - maxrun_G = size(STO_G,3); + load([IdentifDirectoryName '/' fname '_identif_' int2str(ifile_index) '.mat'], 'STO_dSPECTRUM'); + maxrun_dSPECTRUM = size(STO_dSPECTRUM,3); else - maxrun_G = 0; + maxrun_dSPECTRUM = 0; end if ~options_MC.no_identification_minimal - load([IdentifDirectoryName '/' fname '_identif_' int2str(ifile_index) '.mat'], 'STO_D'); - maxrun_D = size(STO_D,3); + load([IdentifDirectoryName '/' fname '_identif_' int2str(ifile_index) '.mat'], 'STO_dMINIMAL'); + maxrun_dMINIMAL = size(STO_dMINIMAL,3); else - maxrun_D = 0; + maxrun_dMINIMAL = 0; end - for irun=1:max([maxrun_dLRE, maxrun_dTAU, maxrun_J, maxrun_G, maxrun_D]) + for irun=1:max([maxrun_dDYNAMIC, maxrun_dREDUCEDFORM, maxrun_dMOMENTS, maxrun_dSPECTRUM, maxrun_dMINIMAL]) iter=iter+1; - % note that this is not the same si_dLREnorm as computed in identification_analysis + % note that this is not the same si_dDYNAMICnorm as computed in identification_analysis % given that we have the MC sample of the Jacobians, we also normalize by the std of the sample of Jacobian entries, to get a fully standardized sensitivity measure - si_dLREnorm(iter,:) = vnorm(STO_si_dLRE(:,:,irun)./repmat(normalize_STO_LRE,1,totparam_nbr-(stderrparam_nbr+corrparam_nbr))).*normaliz1((stderrparam_nbr+corrparam_nbr)+1:end); + si_dDYNAMICnorm(iter,:) = vnorm(STO_si_dDYNAMIC(:,:,irun)./repmat(normalize_STO_DYNAMIC,1,totparam_nbr-(stderrparam_nbr+corrparam_nbr))).*normaliz1((stderrparam_nbr+corrparam_nbr)+1:end); if ~options_MC.no_identification_reducedform - % note that this is not the same si_dTAUnorm as computed in identification_analysis + % note that this is not the same si_dREDUCEDFORMnorm as computed in identification_analysis % given that we have the MC sample of the Jacobians, we also normalize by the std of the sample of Jacobian entries, to get a fully standardized sensitivity measure - si_dTAUnorm(iter,:) = vnorm(STO_si_dTAU(:,:,irun)./repmat(normalize_STO_TAU,1,totparam_nbr)).*normaliz1; + si_dREDUCEDFORMnorm(iter,:) = vnorm(STO_si_dREDUCEDFORM(:,:,irun)./repmat(normalize_STO_REDUCEDFORM,1,totparam_nbr)).*normaliz1; end if ~options_MC.no_identification_moments - % note that this is not the same si_Jnorm as computed in identification_analysis + % note that this is not the same si_dMOMENTSnorm as computed in identification_analysis % given that we have the MC sample of the Jacobians, we also normalize by the std of the sample of Jacobian entries, to get a fully standardized sensitivity measure - si_Jnorm(iter,:) = vnorm(STO_si_J(:,:,irun)./repmat(normalize_STO_MOMENTS,1,totparam_nbr)).*normaliz1; + si_dMOMENTSnorm(iter,:) = vnorm(STO_si_dMOMENTS(:,:,irun)./repmat(normalize_STO_MOMENTS,1,totparam_nbr)).*normaliz1; end if ~options_MC.no_identification_spectrum - % note that this is not the same Gnorm as computed in identification_analysis - Gnorm(iter,:) = vnorm(STO_G(:,:,irun)); %not yet used + % note that this is not the same dSPECTRUMnorm as computed in identification_analysis + dSPECTRUMnorm(iter,:) = vnorm(STO_dSPECTRUM(:,:,irun)); %not yet used end if ~options_MC.no_identification_minimal - % note that this is not the same Dnorm as computed in identification_analysis - Dnorm(iter,:) = vnorm(STO_D(:,:,irun)); %not yet used + % note that this is not the same dMINIMALnorm as computed in identification_analysis + dMINIMALnorm(iter,:) = vnorm(STO_dMINIMAL(:,:,irun)); %not yet used end end end - IDE_LRE.si_dLREnorm = si_dLREnorm; - save([IdentifDirectoryName '/' fname '_identif.mat'], 'pdraws', 'IDE_LRE','STO_LRE','-append'); + IDE_DYNAMIC.si_dDYNAMICnorm = si_dDYNAMICnorm; + save([IdentifDirectoryName '/' fname '_identif.mat'], 'pdraws', 'IDE_DYNAMIC','STO_DYNAMIC','-append'); if ~options_MC.no_identification_reducedform - IDE_REDUCEDFORM.si_dTAUnorm = si_dTAUnorm; - save([IdentifDirectoryName '/' fname '_identif.mat'], 'IDE_REDUCEDFORM', 'STO_TAU','-append'); + IDE_REDUCEDFORM.si_dREDUCEDFORMnorm = si_dREDUCEDFORMnorm; + save([IdentifDirectoryName '/' fname '_identif.mat'], 'IDE_REDUCEDFORM', 'STO_REDUCEDFORM','-append'); end if ~options_MC.no_identification_moments - IDE_MOMENTS.si_Jnorm = si_Jnorm; + IDE_MOMENTS.si_dMOMENTSnorm = si_dMOMENTSnorm; save([IdentifDirectoryName '/' fname '_identif.mat'], 'IDE_MOMENTS', 'STO_MOMENTS','-append'); end @@ -836,25 +853,25 @@ end %% if dynare_identification is called as it own function (not through identification command) and if we load files if nargout>3 && iload filnam = dir([IdentifDirectoryName '/' fname '_identif_*.mat']); - STO_si_dLRE = []; - STO_si_dTAU=[]; - STO_si_J = []; - STO_G = []; - STO_D = []; + STO_si_dDYNAMIC = []; + STO_si_dREDUCEDFORM=[]; + STO_si_dMOMENTS = []; + STO_dSPECTRUM = []; + STO_dMINIMAL = []; for j=1:length(filnam) load([IdentifDirectoryName '/' fname '_identif_',int2str(j),'.mat']); - STO_si_dLRE = cat(3,STO_si_dLRE, STO_si_dLRE(:,abs(iload),:)); + STO_si_dDYNAMIC = cat(3,STO_si_dDYNAMIC, STO_si_dDYNAMIC(:,abs(iload),:)); if ~options_ident.no_identification_reducedform - STO_si_dTAU = cat(3,STO_si_dTAU, STO_si_dTAU(:,abs(iload),:)); + STO_si_dREDUCEDFORM = cat(3,STO_si_dREDUCEDFORM, STO_si_dREDUCEDFORM(:,abs(iload),:)); end if ~options_ident.no_identification_moments - STO_si_J = cat(3,STO_si_J, STO_si_J(:,abs(iload),:)); + STO_si_dMOMENTS = cat(3,STO_si_dMOMENTS, STO_si_dMOMENTS(:,abs(iload),:)); end if ~options_ident.no_identification_spectrum - STO_G = cat(3,STO_G, STO_G(:,abs(iload),:)); + STO_dSPECTRUM = cat(3,STO_dSPECTRUM, STO_dSPECTRUM(:,abs(iload),:)); end if ~options_ident.no_identification_minimal - STO_D = cat(3,STO_D, STO_D(:,abs(iload),:)); + STO_dMINIMAL = cat(3,STO_dMINIMAL, STO_dMINIMAL(:,abs(iload),:)); end end end @@ -865,7 +882,7 @@ if iload disp_identification(ide_hess_point.params, ide_reducedform_point, ide_moments_point, ide_spectrum_point, ide_minimal_point, name, options_ident); if ~options_.nograph % plot (i) identification strength and sensitivity measure based on the sample information matrix and (ii) advanced analysis graphs - plot_identification(ide_hess_point.params, ide_moments_point, ide_hess_point, ide_reducedform_point, ide_lre_point, options_ident.advanced, parameters, name, IdentifDirectoryName, [],name_tex); + plot_identification(ide_hess_point.params, ide_moments_point, ide_hess_point, ide_reducedform_point, ide_dynamic_point, options_ident.advanced, parameters, name, IdentifDirectoryName, [],name_tex); end end @@ -880,7 +897,7 @@ if SampleSize > 1 options_ident.advanced = advanced0; % reset advanced setting if ~options_.nograph % plot (i) identification strength and sensitivity measure based on the sample information matrix and (ii) advanced analysis graphs - plot_identification(pdraws, IDE_MOMENTS, ide_hess_point, IDE_REDUCEDFORM, IDE_LRE, options_ident.advanced, 'MC sample ', name, IdentifDirectoryName, [], name_tex); + plot_identification(pdraws, IDE_MOMENTS, ide_hess_point, IDE_REDUCEDFORM, IDE_DYNAMIC, options_ident.advanced, 'MC sample ', name, IdentifDirectoryName, [], name_tex); end %advanced display and plots for MC Sample, i.e. look at draws with highest/lowest condition number if options_ident.advanced @@ -898,16 +915,16 @@ if SampleSize > 1 disp(['Testing ',tittxt, '.']), if ~iload options_ident.tittxt = tittxt; %title text for graphs and figures - [ide_moments_max, ide_spectrum_max, ide_minimal_max, ide_hess_max, ide_reducedform_max, ide_lre_max, derivatives_info_max, info_max, options_ident] = ... + [ide_moments_max, ide_spectrum_max, ide_minimal_max, ide_hess_max, ide_reducedform_max, ide_dynamic_max, derivatives_info_max, info_max, options_ident] = ... identification_analysis(pdraws(jmax,:), indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1); %the 1 at the end initializes some persistent variables - save([IdentifDirectoryName '/' fname '_identif.mat'], 'ide_hess_max', 'ide_moments_max', 'ide_spectrum_max', 'ide_minimal_max','ide_reducedform_max', 'ide_lre_max', 'jmax', '-append'); + save([IdentifDirectoryName '/' fname '_identif.mat'], 'ide_hess_max', 'ide_moments_max', 'ide_spectrum_max', 'ide_minimal_max','ide_reducedform_max', 'ide_dynamic_max', 'jmax', '-append'); end advanced0 = options_ident.advanced; options_ident.advanced = 1; % make sure advanced setting is on disp_identification(pdraws(jmax,:), ide_reducedform_max, ide_moments_max, ide_spectrum_max, ide_minimal_max, name, options_ident); options_ident.advanced = advanced0; %reset advanced setting if ~options_.nograph % plot (i) identification strength and sensitivity measure based on the sample information matrix and (ii) advanced analysis graphs - plot_identification(pdraws(jmax,:), ide_moments_max, ide_hess_max, ide_reducedform_max, ide_lre_max, 1, tittxt, name, IdentifDirectoryName, tittxt, name_tex); + plot_identification(pdraws(jmax,:), ide_moments_max, ide_hess_max, ide_reducedform_max, ide_dynamic_max, 1, tittxt, name, IdentifDirectoryName, tittxt, name_tex); end % SMALLEST condition number @@ -918,16 +935,16 @@ if SampleSize > 1 disp(['Testing ',tittxt, '.']), if ~iload options_ident.tittxt = tittxt; %title text for graphs and figures - [ide_moments_min, ide_spectrum_min, ide_minimal_min, ide_hess_min, ide_reducedform_min, ide_lre_min, derivatives_info_min, info_min, options_ident] = ... + [ide_moments_min, ide_spectrum_min, ide_minimal_min, ide_hess_min, ide_reducedform_min, ide_dynamic_min, derivatives_info_min, info_min, options_ident] = ... identification_analysis(pdraws(jmin,:), indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1); %the 1 at the end initializes persistent variables - save([IdentifDirectoryName '/' fname '_identif.mat'], 'ide_hess_min', 'ide_moments_min','ide_spectrum_min','ide_minimal_min','ide_reducedform_min', 'ide_lre_min', 'jmin', '-append'); + save([IdentifDirectoryName '/' fname '_identif.mat'], 'ide_hess_min', 'ide_moments_min','ide_spectrum_min','ide_minimal_min','ide_reducedform_min', 'ide_dynamic_min', 'jmin', '-append'); end advanced0 = options_ident.advanced; options_ident.advanced = 1; % make sure advanced setting is on disp_identification(pdraws(jmin,:), ide_reducedform_min, ide_moments_min, ide_spectrum_min, ide_minimal_min, name, options_ident); options_ident.advanced = advanced0; %reset advanced setting if ~options_.nograph % plot (i) identification strength and sensitivity measure based on the sample information matrix and (ii) advanced analysis graphs - plot_identification(pdraws(jmin,:),ide_moments_min,ide_hess_min,ide_reducedform_min,ide_lre_min,1,tittxt,name,IdentifDirectoryName,tittxt,name_tex); + plot_identification(pdraws(jmin,:),ide_moments_min,ide_hess_min,ide_reducedform_min,ide_dynamic_min,1,tittxt,name,IdentifDirectoryName,tittxt,name_tex); end % reset nodisplay option options_.nodisplay = store_nodisplay; @@ -941,7 +958,7 @@ if SampleSize > 1 disp(['Testing ',tittxt, '.']), if ~iload options_ident.tittxt = tittxt; %title text for graphs and figures - [ide_moments_(j), ide_spectrum_(j), ide_minimal_(j), ide_hess_(j), ide_reducedform_(j), ide_lre_(j), derivatives_info_(j), info_resolve, options_ident] = ... + [ide_moments_(j), ide_spectrum_(j), ide_minimal_(j), ide_hess_(j), ide_reducedform_(j), ide_dynamic_(j), derivatives_info_(j), info_resolve, options_ident] = ... identification_analysis(pdraws(jcrit(j),:), indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, 1); end advanced0 = options_ident.advanced; options_ident.advanced = 1; %make sure advanced setting is on @@ -949,11 +966,11 @@ if SampleSize > 1 options_ident.advanced = advanced0; % reset advanced if ~options_.nograph % plot (i) identification strength and sensitivity measure based on the sample information matrix and (ii) advanced analysis graphs - plot_identification(pdraws(jcrit(j),:), ide_moments_(j), ide_hess_(j), ide_reducedform_(j), ide_lre_(j), 1, tittxt, name, IdentifDirectoryName, tittxt, name_tex); + plot_identification(pdraws(jcrit(j),:), ide_moments_(j), ide_hess_(j), ide_reducedform_(j), ide_dynamic_(j), 1, tittxt, name, IdentifDirectoryName, tittxt, name_tex); end end if ~iload - save([IdentifDirectoryName '/' fname '_identif.mat'], 'ide_hess_', 'ide_moments_', 'ide_reducedform_', 'ide_lre_', 'ide_spectrum_', 'ide_minimal_', 'jcrit', '-append'); + save([IdentifDirectoryName '/' fname '_identif.mat'], 'ide_hess_', 'ide_moments_', 'ide_reducedform_', 'ide_dynamic_', 'ide_spectrum_', 'ide_minimal_', 'jcrit', '-append'); end % reset nodisplay option options_.nodisplay = store_nodisplay; diff --git a/matlab/get_identification_jacobians.m b/matlab/get_identification_jacobians.m index 20e919a78..8d66d438f 100644 --- a/matlab/get_identification_jacobians.m +++ b/matlab/get_identification_jacobians.m @@ -1,16 +1,9 @@ -function [J, G, D, dTAU, dLRE, dA, dOm, dYss, MOMENTS] = get_identification_jacobians(A, B, estim_params, M, oo, options, options_ident, indpmodel, indpstderr, indpcorr, indvobs) -% [J, G, D, dTAU, dLRE, dA, dOm, dYss, MOMENTS] = get_identification_jacobians(A, B, estim_params, M, oo, options, options_ident, indpmodel, indpstderr, indpcorr, indvobs) +function [E_y, dE_y, 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, dMINIMAL, derivatives_info] = get_identification_jacobians(estim_params, M, oo, options, options_ident, indpmodel, indpstderr, indpcorr, indvobs) % previously getJJ.m -% ------------------------------------------------------------------------- -% Sets up the Jacobians needed for identification analysis based on the -% Iskrev's J, Qu and Tkachenko's G and Komunjer and Ng's D matrices as well -% as on the reduced-form model (dTAU) and the dynamic model (dLRE) +% % Sets up the Jacobians needed for identification analysis % ========================================================================= % INPUTS -% A: [endo_nbr by endo_nbr] in DR order -% Transition matrix from Kalman filter for all endogenous declared variables, -% B: [endo_nbr by exo_nbr] in DR order -% Transition matrix from Kalman filter mapping shocks today to endogenous variables today % estim_params: [structure] storing the estimation information % M: [structure] storing the model information % oo: [structure] storing the reduced-form solution results @@ -31,37 +24,50 @@ function [J, G, D, dTAU, dLRE, dA, dOm, dYss, MOMENTS] = get_identification_jaco % indvobs: [obs_nbr by 1] index of observed (VAROBS) variables % ------------------------------------------------------------------------- % OUTPUTS -% J: [obs_nbr+obs_nbr*(obs_nbr+1)/2+nlags*obs_nbr^2 by totparam_nbr] in DR Order -% Jacobian of 1st and 2nd order moments of observables wrt, -% all parameters, i.e. dgam (see below for definition of gam). -% Corresponds to Iskrev (2010)'s J matrix. -% G: [totparam_nbr by totparam_nbr] in DR Order -% Sum of (1) Gram Matrix of Jacobian of spectral density -% wrt to all parameters plus (2) Gram Matrix of Jacobian -% of steady state wrt to all parameters. -% Corresponds to Qu and Tkachenko (2012)'s G matrix. -% D: [obs_nbr+minstate_nbr^2+minstate_nbr*exo_nbr+obs_nbr*minstate_nbr+obs_nbr*exo_nbr+exo_nbr*(exo_nbr+1)/2 by totparam_nbr+minstate_nbr^2+exo_nbr^2] in DR order -% Jacobian of minimal System Matrices and unique Transformation -% of steady state and spectral density matrix. -% Corresponds to Komunjer and Ng (2011)'s Delta matrix. -% dTAU: [(endo_nbr+endo_nbr^2+endo_nbr*(endo_nbr+1)/2) by totparam_nbr] in DR order -% Jacobian of linearized reduced form state space model, given Yss [steady state], -% A [transition matrix], B [matrix of shocks], Sigma [covariance of shocks] -% tau = [ys; vec(A); dyn_vech(B*Sigma*B')] with respect to all parameters. -% dLRE: [endo_nbr+endo_nbr*(dynamicvar_nbr+exo_nbr) by modparam_nbr] in DR order -% Jacobian of steady state and linear rational expectation matrices -% (i.e. Jacobian of dynamic model) with respect to estimated model parameters only (indpmodel) -% dA: [endo_nbr by endo_nbr by totparam_nbr] in DR order -% Jacobian (wrt to all parameters) of transition matrix A -% dOm: [endo_nbr by endo_nbr by totparam_nbr] in DR order -% Jacobian (wrt to all paramters) of Om = (B*Sigma_e*B') -% dYss [endo_nbr by modparam_nbr] in DR order -% Jacobian (wrt model parameters only) of steady state -% MOMENTS: [obs_nbr+obs_nbr*(obs_nbr+1)/2+nlags*obs_nbr^2 by 1] -% vector of theoretical moments of observed (VAROBS) -% variables. Note that J is the Jacobian of MOMENTS. -% MOMENTS = [ys(indvobs); dyn_vech(GAM{1}); vec(GAM{j+1})]; for j=1:ar and where -% GAM is the first output of th_autocovariances +% +% MEAN [endo_nbr by 1], in DR order. Expectation of all model variables +% * order==1: corresponds to steady state +% * order==2: computed from pruned state space system (as in e.g. Andreasen, Fernandez-Villaverde, Rubio-Ramirez, 2018) +% dMEAN [endo_nbr by totparam_nbr], in DR Order, Jacobian (wrt all params) of MEAN +% +% REDUCEDFORM [rowredform_nbr by 1] in DR order. Steady state and reduced-form model solution matrices for all model variables +% * order==1: [Yss' vec(dghx)' vech(dghu*Sigma_e*dghu')']', +% where rowredform_nbr = endo_nbr+endo_nbr*nspred+endo_nbr*(endo_nbr+1)/2 +% * order==2: [Yss' vec(dghx)' vech(dghu*Sigma_e*dghu')' vec(dghxx)' vec(dghxu)' vec(dghuu)' vec(dghs2)']', +% where rowredform_nbr = endo_nbr+endo_nbr*nspred+endo_nbr*(endo_nbr+1)/2+endo_nbr*nspred^2*endo_nbr*nspred*exo_nr+endo_nbr*exo_nbr^2+endo_nbr +% dREDUCEDFORM: [rowrf_nbr by totparam_nbr] in DR order, Jacobian (wrt all params) of REDUCEDFORM +% * order==1: corresponds to Iskrev (2010)'s J_2 matrix +% * order==2: corresponds to Mutschler (2015)'s J matrix +% +% DYNAMIC [rowdyn_nbr by 1] in declaration order. Steady state and dynamic model derivatives for all model variables +% * order==1: [ys' vec(g1)']', rowdyn_nbr=endo_nbr+length(g1) +% * order==2: [ys' vec(g1)' vec(g2)']', rowdyn_nbr=endo_nbr+length(g1)+length(g2) +% dDYNAMIC [rowdyn_nbr by modparam_nbr] in declaration order. Jacobian (wrt model parameters) of DYNAMIC +% * order==1: corresponds to Ratto and Iskrev (2011)'s J_\Gamma matrix (or LRE) +% * order==2: additionally includes derivatives of dynamic hessian +% +% MOMENTS: [obs_nbr+obs_nbr*(obs_nbr+1)/2+nlags*obs_nbr^2 by 1] in DR order. First two theoretical moments for VAROBS variables, i.e. +% [E[yobs]' vech(E[yobs*yobs'])' vec(E[yobs*yobs(-1)'])' ... vec(E[yobs*yobs(-nlag)'])'] +% dMOMENTS: [obs_nbr+obs_nbr*(obs_nbr+1)/2+nlags*obs_nbr^2 by totparam_nbr] in DR order. Jacobian (wrt all params) of MOMENTS +% * 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 +% spectral density at frequency w: f(w) = (2*pi)^(-1)*H(exp(-i*w))*E[xi*xi']*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 xi and H are computed from linear state space system +% * order==2: corresponds to Mutschler (2015)'s G_2 matrix, where xi and H are computed from pruned state space system +% +% dMINIMAL: [obs_nbr+minx^2+minx*exo_nbr+obs_nbr*minx+obs_nbr*exo_nbr+exo_nbr*(exo_nbr+1)/2 by totparam_nbr+minx^2+exo_nbr^2] +% Jacobian (wrt all params, and similarity_transformation_matrices (T and U)) of observational equivalent minimal ABCD system, +% corresponds to Komunjer and Ng (2011)'s Deltabar matrix, where +% MINIMAL = [vec(E[yobs]' vec(minA)' vec(minB)' vec(minC)' vec(minD)' vech(Sigma_e)']' +% * order==1: E[yobs] is equal to steady state +% * order==2: E[yobs] is computed from the pruned state space system (second-order accurate), as noted in section 5 of Komunjer and Ng (2011) +% +% derivatives_info [structure] for use in dsge_likelihood to compute Hessian analytically. +% Contains dA, dB, and d(B*Sigma_e*B'), where A and B are from Kalman filter. Only used at order==1. +% % ------------------------------------------------------------------------- % This function is called by % * identification_analysis.m @@ -71,11 +77,10 @@ function [J, G, D, dTAU, dLRE, dA, dOm, dYss, MOMENTS] = get_identification_jaco % * get_minimal_state_representation % * duplication % * dyn_vech -% * get_perturbation_params_deriv (previously getH) +% * get_perturbation_params_derivs (previously getH) % * get_all_parameters % * fjaco % * lyapunov_symm -% * th_autocovariances % * identification_numerical_objective (previously thet2tau) % * vec % ========================================================================= @@ -98,13 +103,17 @@ function [J, G, D, dTAU, dLRE, dA, dOm, dYss, MOMENTS] = get_identification_jaco % ========================================================================= %get options +no_identification_moments = options_ident.no_identification_moments; +no_identification_minimal = options_ident.no_identification_minimal; +no_identification_spectrum = options_ident.no_identification_spectrum; + +order = options_ident.order; nlags = options_ident.ar; useautocorr = options_ident.useautocorr; grid_nbr = options_ident.grid_nbr; kronflag = options_ident.analytic_derivation_mode; -no_identification_moments = options_ident.no_identification_moments; -no_identification_minimal = options_ident.no_identification_minimal; -no_identification_spectrum = options_ident.no_identification_spectrum; + +% set values params0 = M.params; %values at which to evaluate dynamic, static and param_derivs files Sigma_e0 = M.Sigma_e; %covariance matrix of exogenous shocks Corr_e0 = M.Correlation_matrix; %correlation matrix of exogenous shocks @@ -114,6 +123,7 @@ if isempty(para0) %if there is no estimated_params block, consider all stderr and all model parameters, but no corr parameters para0 = [stderr_e0', params0']; end + %get numbers/lengths of vectors modparam_nbr = length(indpmodel); stderrparam_nbr = length(indpstderr); @@ -122,108 +132,189 @@ totparam_nbr = stderrparam_nbr + corrparam_nbr + modparam_nbr; obs_nbr = length(indvobs); exo_nbr = M.exo_nbr; endo_nbr = M.endo_nbr; +nspred = M.nspred; +nstatic = M.nstatic; +indvall = (1:endo_nbr)'; %index for all endogenous variables +d2flag = 0; % do not compute second parameter derivatives -%% Construct dTAU, dLRE, dA, dOm, dYss -DERIVS = get_perturbation_params_derivs(M, options, estim_params, oo, indpmodel, indpstderr, indpcorr, 0);%d2flag=0 -dA = zeros(M.endo_nbr, M.endo_nbr, size(DERIVS.dghx,3)); -dA(:,M.nstatic+(1:M.nspred),:) = DERIVS.dghx; -dB = DERIVS.dghu; -dSig = DERIVS.dSigma_e; -dOm = DERIVS.dOm; -dYss = DERIVS.dYss; +% Get Jacobians (wrt selected params) of steady state, dynamic model derivatives and perturbation solution matrices for all endogenous variables +oo.dr.derivs = get_perturbation_params_derivs(M, options, estim_params, oo, indpmodel, indpstderr, indpcorr, d2flag); -% Collect terms for derivative of tau=[Yss; vec(A); vech(Om)] -dTAU = zeros(endo_nbr*endo_nbr+endo_nbr*(endo_nbr+1)/2, totparam_nbr); -for j=1:totparam_nbr - dTAU(:,j) = [vec(dA(:,:,j)); dyn_vech(dOm(:,:,j))]; +[I,~] = find(M.lead_lag_incidence'); %I is used to select nonzero columns of the Jacobian of endogenous variables in dynamic model files +yy0 = oo.dr.ys(I); %steady state of dynamic (endogenous and auxiliary variables) in lead_lag_incidence order +Yss = oo.dr.ys(oo.dr.order_var); % steady state in DR order +if order == 1 + [~, g1 ] = feval([M.fname,'.dynamic'], yy0, oo.exo_steady_state', params0, oo.dr.ys, 1); + %g1 is [endo_nbr by yy0ex0_nbr first derivative (wrt all dynamic variables) of dynamic model equations, i.e. df/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order + DYNAMIC = [Yss; + vec(g1(oo.dr.order_var,:))]; %add steady state and put rows of g1 in DR order + dDYNAMIC = [oo.dr.derivs.dYss; + reshape(oo.dr.derivs.dg1(oo.dr.order_var,:,:),size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2),size(oo.dr.derivs.dg1,3)) ]; %reshape dg1 in DR order and add steady state + REDUCEDFORM = [Yss; + vec(oo.dr.ghx); + dyn_vech(oo.dr.ghu*Sigma_e0*transpose(oo.dr.ghu))]; %in DR order + dREDUCEDFORM = zeros(endo_nbr*nspred+endo_nbr*(endo_nbr+1)/2, totparam_nbr); + for j=1:totparam_nbr + dREDUCEDFORM(:,j) = [vec(oo.dr.derivs.dghx(:,:,j)); + dyn_vech(oo.dr.derivs.dOm(:,:,j))]; + end + dREDUCEDFORM = [ [zeros(endo_nbr, stderrparam_nbr+corrparam_nbr) oo.dr.derivs.dYss]; dREDUCEDFORM ]; % add steady state + +elseif order == 2 + [~, g1, g2 ] = feval([M.fname,'.dynamic'], yy0, oo.exo_steady_state', params0, oo.dr.ys, 1); + %g1 is [endo_nbr by yy0ex0_nbr first derivative (wrt all dynamic variables) of dynamic model equations, i.e. df/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order + %g2 is [endo_nbr by yy0ex0_nbr^2] second derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order + DYNAMIC = [Yss; + vec(g1(oo.dr.order_var,:)); + vec(g2(oo.dr.order_var,:))]; %add steady state and put rows of g1 and g2 in DR order + dDYNAMIC = [oo.dr.derivs.dYss; + reshape(oo.dr.derivs.dg1(oo.dr.order_var,:,:),size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2),size(oo.dr.derivs.dg1,3)); %reshape dg1 in DR order + reshape(oo.dr.derivs.dg2(oo.dr.order_var,:),size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2)^2,size(oo.dr.derivs.dg1,3))]; %reshape dg2 in DR order + REDUCEDFORM = [Yss; + vec(oo.dr.ghx); + dyn_vech(oo.dr.ghu*Sigma_e0*transpose(oo.dr.ghu)); + vec(oo.dr.ghxx); + vec(oo.dr.ghxu); + vec(oo.dr.ghuu); + vec(oo.dr.ghs2)]; %in DR order + dREDUCEDFORM = zeros(endo_nbr*nspred+endo_nbr*(endo_nbr+1)/2+endo_nbr*nspred^2+endo_nbr*nspred*exo_nbr+endo_nbr*exo_nbr^2+endo_nbr, totparam_nbr); + for j=1:totparam_nbr + dREDUCEDFORM(:,j) = [vec(oo.dr.derivs.dghx(:,:,j)); + dyn_vech(oo.dr.derivs.dOm(:,:,j)); + vec(oo.dr.derivs.dghxx(:,:,j)); + vec(oo.dr.derivs.dghxu(:,:,j)); + vec(oo.dr.derivs.dghuu(:,:,j)); + vec(oo.dr.derivs.dghs2(:,j))]; + end + dREDUCEDFORM = [ [zeros(endo_nbr, stderrparam_nbr+corrparam_nbr) oo.dr.derivs.dYss]; dREDUCEDFORM ]; % add steady state +elseif order == 3 + [~, g1, g2, g3 ] = feval([M.fname,'.dynamic'], yy0, oo.exo_steady_state', params0, oo.dr.ys, 1); + %g1 is [endo_nbr by yy0ex0_nbr first derivative (wrt all dynamic variables) of dynamic model equations, i.e. df/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order + %g2 is [endo_nbr by yy0ex0_nbr^2] second derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order + DYNAMIC = [Yss; + vec(g1(oo.dr.order_var,:)); + vec(g2(oo.dr.order_var,:)); + vec(g3(oo.dr.order_var,:))]; %add steady state and put rows of g1 and g2 in DR order + dDYNAMIC = [oo.dr.derivs.dYss; + reshape(oo.dr.derivs.dg1(oo.dr.order_var,:,:),size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2),size(oo.dr.derivs.dg1,3)); %reshape dg1 in DR order + reshape(oo.dr.derivs.dg2(oo.dr.order_var,:),size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2)^2,size(oo.dr.derivs.dg1,3)); + reshape(oo.dr.derivs.dg2(oo.dr.order_var,:),size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2)^2,size(oo.dr.derivs.dg1,3))]; %reshape dg3 in DR order + REDUCEDFORM = [Yss; + vec(oo.dr.ghx); + dyn_vech(oo.dr.ghu*Sigma_e0*transpose(oo.dr.ghu)); + vec(oo.dr.ghxx); vec(oo.dr.ghxu); vec(oo.dr.ghuu); vec(oo.dr.ghs2); + vec(oo.dr.ghxxx); vec(oo.dr.ghxxu); vec(oo.dr.ghxuu); vec(oo.dr.ghuuu); vec(oo.dr.ghxss); vec(oo.dr.ghuss)]; %in DR order + dREDUCEDFORM = zeros(size(REDUCEDFORM,1)-endo_nbr, totparam_nbr); + for j=1:totparam_nbr + dREDUCEDFORM(:,j) = [vec(oo.dr.derivs.dghx(:,:,j)); + dyn_vech(oo.dr.derivs.dOm(:,:,j)); + vec(oo.dr.derivs.dghxx(:,:,j)); vec(oo.dr.derivs.dghxu(:,:,j)); vec(oo.dr.derivs.dghuu(:,:,j)); vec(oo.dr.derivs.dghs2(:,j)) + vec(oo.dr.derivs.dghxxx(:,:,j)); vec(oo.dr.derivs.dghxxu(:,:,j)); vec(oo.dr.derivs.dghxuu(:,:,j)); vec(oo.dr.derivs.dghuuu(:,:,j)); vec(oo.dr.derivs.dghxss(:,:,j)); vec(oo.dr.derivs.dghuss(:,:,j))]; + end + dREDUCEDFORM = [ [zeros(endo_nbr, stderrparam_nbr+corrparam_nbr) oo.dr.derivs.dYss]; dREDUCEDFORM ]; % add steady state end -dTAU = [ [zeros(endo_nbr, stderrparam_nbr+corrparam_nbr) dYss]; dTAU ]; % add steady state -dLRE = [dYss; reshape(DERIVS.dg1,size(DERIVS.dg1,1)*size(DERIVS.dg1,2),size(DERIVS.dg1,3)) ]; %reshape dg1 and add steady state - -%State Space Matrices for VAROBS variables -C = A(indvobs,:); -dC = dA(indvobs,:,:); -D = B(indvobs,:); -dD = dB(indvobs,:,:); - -%% Iskrev (2010) + +% 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); +E_y = oo.dr.pruned.E_y; dE_y = 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; + +%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 + dT = zeros(endo_nbr,endo_nbr,totparam_nbr); + dT(:,(nstatic+1):(nstatic+nspred),:) = oo.dr.derivs.dghx; + derivatives_info.DT = dT; + derivatives_info.DOm = oo.dr.derivs.dOm; + derivatives_info.DYss = oo.dr.derivs.dYss; +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(para0, 1, estim_params, M, oo, options, indpmodel, indpstderr, indpcorr, indvobs, useautocorr, nlags, grid_nbr); %[outputflag=1] + MOMENTS = [E_y; MOMENTS]; if kronflag == -1 - %numerical derivative of autocovariogram [outputflag=1] - J = fjaco(str2func('identification_numerical_objective'), para0, 1, estim_params, M, oo, options, indpmodel, indpstderr, indpcorr, indvobs, useautocorr, nlags, grid_nbr); + %numerical derivative of autocovariogram + dMOMENTS = fjaco(str2func('identification_numerical_objective'), para0, 1, estim_params, M, oo, options, indpmodel, indpstderr, indpcorr, indvobs, useautocorr, nlags, grid_nbr); %[outputflag=1] M.params = params0; %make sure values are set back M.Sigma_e = Sigma_e0; %make sure values are set back M.Correlation_matrix = Corr_e0 ; %make sure values are set back - J = [[zeros(obs_nbr,stderrparam_nbr+corrparam_nbr) dYss(indvobs,:)]; J]; %add Jacobian of steady state of VAROBS variables + dMOMENTS = [dE_y; dMOMENTS]; %add Jacobian of steady state of VAROBS variables else - J = zeros(obs_nbr + obs_nbr*(obs_nbr+1)/2 + nlags*obs_nbr^2 , totparam_nbr); - J(1:obs_nbr,stderrparam_nbr+corrparam_nbr+1 : totparam_nbr) = dYss(indvobs,:); %add Jacobian of steady state of VAROBS variables - % Denote Ezz0 = E_t(z_t * z_t'), then the following Lyapunov equation defines the autocovariagram: Ezz0 -A*Ezz*A' = B*Sig_e*B' = Om - Gamma_y = lyapunov_symm(A, B*Sigma_e0*B', options.lyapunov_fixed_point_tol, options.qz_criterium, options.lyapunov_complex_threshold, 1, options.debug); + dMOMENTS = zeros(obs_nbr + obs_nbr*(obs_nbr+1)/2 + nlags*obs_nbr^2 , totparam_nbr); + dMOMENTS(1:obs_nbr,:) = dE_y; %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 = lyapunov_symm(A, B*Varinov*B', options.lyapunov_fixed_point_tol, options.qz_criterium, options.lyapunov_complex_threshold, 1, options.debug); + Eyy0 = C*Ezz0*C' + D*Varinov*D'; %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(Gamma_y) < 1e-12); %find values that are numerical zero - Gamma_y(indzeros) = 0; - % if useautocorr, - sdy = sqrt(diag(Gamma_y)); %theoretical standard deviation - sy = sdy*sdy'; %cross products of standard deviations - % end - for j = 1:(stderrparam_nbr+corrparam_nbr) - %Jacobian of Ezz0 wrt exogenous paramters: dEzz0(:,:,j)-A*dEzz0(:,:,j)*A'=dOm(:,:,j), because dA is zero by construction for stderr and corr parameters - dEzz0 = lyapunov_symm(A,dOm(:,:,j),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 - indzeros = find(abs(dEzz0) < 1e-12); - dEzz0(indzeros) = 0; - if useautocorr - dsy = 1/2./sdy.*diag(dEzz0); - dsy = dsy*sdy'+sdy*dsy'; - dEzz0corr = (dEzz0.*sy-dsy.*Gamma_y)./(sy.*sy); - dEzz0corr = dEzz0corr-diag(diag(dEzz0corr))+diag(diag(dEzz0)); - J(obs_nbr+1 : obs_nbr+obs_nbr*(obs_nbr+1)/2 , j) = dyn_vech(dEzz0corr(indvobs,indvobs)); %focus only on VAROBS variables - else - J(obs_nbr+1 : obs_nbr+obs_nbr*(obs_nbr+1)/2 , j) = dyn_vech(dEzz0(indvobs,indvobs)); %focus only on VAROBS variables - end - %Jacobian of Ezzi = E_t(z_t * z_{t-i}'): dEzzi(:,:,j) = A^i*dEzz0(:,:,j) wrt stderr and corr parameters, because dA is zero by construction for stderr and corr parameters - for i = 1:nlags - dEzzi = A^i*dEzz0; - if useautocorr - dEzzi = (dEzzi.*sy-dsy.*(A^i*Gamma_y))./(sy.*sy); - end - J(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, j) = vec(dEzzi(indvobs,indvobs)); %focus only on VAROBS variables - end + indzeros = find(abs(Eyy0) < 1e-12); %find values that are numerical zero + Eyy0(indzeros) = 0; + if useautocorr + sdy = sqrt(diag(Eyy0)); %theoretical standard deviation + sy = sdy*sdy'; %cross products of standard deviations end - for j=1:modparam_nbr - %Jacobian of Ezz0 wrt model parameters: dEzz0(:,:,j) - A*dEzz0(:,:,j)*A' = dOm(:,:,j) + dA(:,:,j)*Ezz*A'+ A*Ezz*dA(:,:,j)' - dEzz0 = lyapunov_symm(A,dA(:,:,j+stderrparam_nbr+corrparam_nbr)*Gamma_y*A'+A*Gamma_y*dA(:,:,j+stderrparam_nbr+corrparam_nbr)'+dOm(:,:,j+stderrparam_nbr+corrparam_nbr),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 - indzeros = find(abs(dEzz0) < 1e-12); - dEzz0(indzeros) = 0; + + 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(dEzz0); + dsy = 1/2./sdy.*diag(dEyy0); dsy = dsy*sdy'+sdy*dsy'; - dEzz0corr = (dEzz0.*sy-dsy.*Gamma_y)./(sy.*sy); - dEzz0corr = dEzz0corr-diag(diag(dEzz0corr))+diag(diag(dEzz0)); - J(obs_nbr+1 : obs_nbr+obs_nbr*(obs_nbr+1)/2 , stderrparam_nbr+corrparam_nbr+j) = dyn_vech(dEzz0corr(indvobs,indvobs)); %focus only on VAROBS variables + 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 - J(obs_nbr+1 : obs_nbr+obs_nbr*(obs_nbr+1)/2 , stderrparam_nbr+corrparam_nbr+j) = dyn_vech(dEzz0(indvobs,indvobs)); %focus only on VAROBS variables + dMOMENTS(obs_nbr+1 : obs_nbr+obs_nbr*(obs_nbr+1)/2 , jp) = dyn_vech(dEyy0); %focus only on VAROBS variables end - %Jacobian of Ezzi = E_t(z_t * z_{t-i}'): dEzzi(:,:,j) = A^i*dEzz0(:,:,j) + d(A^i)*dEzz0(:,:,j) wrt model parameters - for i = 1:nlags - dEzzi = A^i*dEzz0; - for ii=1:i - dEzzi = dEzzi + A^(ii-1)*dA(:,:,j+stderrparam_nbr+corrparam_nbr)*A^(i-ii)*Gamma_y; + tmpEyyi = A*Ezz0*C' + B*Varinov*D'; + %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) + for i = 1:nlags + Eyyi = C*Ai*tmpEyyi; + dEyyi = dC(:,:,jp)*Ai*tmpEyyi + C*dAi*tmpEyyi + C*Ai*dtmpEyyi; + if useautocorr + dEyyi = (dEyyi.*sy-dsy.*Eyyi)./(sy.*sy); end - if useautocorr - dEzzi = (dEzzi.*sy-dsy.*(A^i*Gamma_y))./(sy.*sy); - end - J(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, j+stderrparam_nbr+corrparam_nbr) = vec(dEzzi(indvobs,indvobs)); %focus only on VAROBS variables + 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) end end end else - J = []; + MOMENTS = []; + dMOMENTS = []; end - -%% Qu and Tkachenko (2012) + +%% Compute dSPECTRUM +%Note that our state space system for computing the spectrum 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_spectrum %Some info on the spectral density: Dynare's state space system is given for states (z_{t} = A*z_{t-1} + B*u_{t}) and observables (y_{t} = C*z_{t-1} + D*u_{t}) %The spectral density for y_{t} can be computed using different methods, which are numerically equivalent @@ -259,8 +350,8 @@ if ~no_identification_spectrum tneg = exp(-sqrt(-1)*freqs); %negative Fourier frequencies IA = eye(size(A,1)); if kronflag == -1 - %numerical derivative of spectral density [outputflag=2] - dOmega_tmp = fjaco(str2func('identification_numerical_objective'), para0, 2, estim_params, M, oo, options, indpmodel, indpstderr, indpcorr, indvobs, useautocorr, nlags, grid_nbr); + %numerical derivative of spectral density + dOmega_tmp = fjaco(str2func('identification_numerical_objective'), para0, 2, estim_params, M, oo, options, indpmodel, indpstderr, indpcorr, indvobs, useautocorr, nlags, grid_nbr); %[outputflag=2] M.params = params0; %make sure values are set back M.Sigma_e = Sigma_e0; %make sure values are set back M.Correlation_matrix = Corr_e0 ; %make sure values are set back @@ -269,72 +360,72 @@ 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 - G = dOmega'*dOmega; + dSPECTRUM = dOmega'*dOmega; else % due to symmetry to negative frequencies we can add positive frequencies twice - G = G + 2*(dOmega'*dOmega); + dSPECTRUM = dSPECTRUM + 2*(dOmega'*dOmega); end end - elseif kronflag == 1 + 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)); - dSig = reshape(dSig,size(dSig,1)*size(dSig,2),size(dSig,3)); - K_obs_exo = commutation(obs_nbr,exo_nbr); + dVarinov = reshape(dVarinov,size(dVarinov,1)*size(dVarinov,2),size(dVarinov,3)); + K_obs_exo = commutation(obs_nbr,size(Varinov,1)); for ig=1:length(freqs) z = tneg(ig); zIminusA = (z*IA - A); zIminusAinv = zIminusA\IA; Transferfct = D + C*zIminusAinv*B; % Transfer function dzIminusA = -dA; - dzIminusAinv = kron(-(transpose(zIminusA)\IA),zIminusAinv)*dzIminusA; - dTransferfct = dD + DerivABCD(C,dC,zIminusAinv,dzIminusAinv,B,dB); + dzIminusAinv = kron(-(transpose(zIminusA)\IA),zIminusAinv)*dzIminusA; %this takes long + dTransferfct = dD + DerivABCD(C,dC,zIminusAinv,dzIminusAinv,B,dB); %also long dTransferfct_conjt = K_obs_exo*conj(dTransferfct); - dOmega = (1/(2*pi))*DerivABCD(Transferfct,dTransferfct,Sigma_e0,dSig,Transferfct',dTransferfct_conjt); + dOmega = (1/(2*pi))*DerivABCD(Transferfct,dTransferfct,Varinov,dVarinov,Transferfct',dTransferfct_conjt); %also long if ig == 1 % add zero frequency once - G = dOmega'*dOmega; + dSPECTRUM = dOmega'*dOmega; else % due to symmetry to negative frequencies we can add positive frequencies twice - G = G + 2*(dOmega'*dOmega); + dSPECTRUM = dSPECTRUM + 2*(dOmega'*dOmega); end end %put back into tensor notation - dA = reshape(dA,endo_nbr,endo_nbr,totparam_nbr); - dB = reshape(dB,endo_nbr,exo_nbr,totparam_nbr); - dC = reshape(dC,obs_nbr,endo_nbr,totparam_nbr); - dD = reshape(dD,obs_nbr,exo_nbr,totparam_nbr); - dSig = reshape(dSig,exo_nbr,exo_nbr,totparam_nbr); + 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(endo_nbr); + invIzminusA = IzminusA\eye(size(A,1)); Transferfct = D + C*invIzminusA*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*dSig(:,:,j)*Transferfct'; + dOmega_tmp = Transferfct*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*M.Sigma_e*Transferfct' + Transferfct*dSig(:,:,j)*Transferfct' + Transferfct*M.Sigma_e*dTransferfct'; + dOmega_tmp = dTransferfct*Varinov*Transferfct' + Transferfct*dVarinov(:,:,j)*Transferfct' + Transferfct*Varinov*dTransferfct'; end dOmega(:,j) = (1/(2*pi))*dOmega_tmp(:); end if ig == 1 % add zero frequency once - G = dOmega'*dOmega; + dSPECTRUM = dOmega'*dOmega; else % due to symmetry to negative frequencies we can add positive frequencies twice - G = G + 2*(dOmega'*dOmega); + dSPECTRUM = dSPECTRUM + 2*(dOmega'*dOmega); end end end % Normalize Matrix and add steady state Jacobian, note that G is real and symmetric by construction - G = 2*pi*G./(2*length(freqs)-1) + [zeros(obs_nbr,stderrparam_nbr+corrparam_nbr) dYss(indvobs,:)]'*[zeros(obs_nbr,stderrparam_nbr+corrparam_nbr) dYss(indvobs,:)]; - G = real(G); + dSPECTRUM = 2*pi*dSPECTRUM./(2*length(freqs)-1) + dE_y'*dE_y; + dSPECTRUM = real(dSPECTRUM); else - G = []; + dSPECTRUM = []; end -%% Komunjer and Ng (2012) +%% Compute dMINIMAL if ~no_identification_minimal if obs_nbr < exo_nbr % Check whether criteria can be used @@ -342,36 +433,30 @@ if ~no_identification_minimal warning_KomunjerNg = [warning_KomunjerNg ' There are more shocks and measurement errors than observables, this is not implemented (yet).\n']; warning_KomunjerNg = [warning_KomunjerNg ' Skip identification analysis based on minimal state space system.\n']; fprintf(warning_KomunjerNg); - D = []; + dMINIMAL = []; else - % Derive and check minimal state - if isfield(oo.dr,'state_var') - state_var = oo.dr.state_var; %state variables in declaration order - else - % DR-order: static variables first, then purely backward variables, then mixed variables, finally purely forward variables. - % Inside each category, variables are arranged according to the declaration order. - % state variables are the purely backward variables and the mixed variables - state_var = transpose(oo.dr.order_var( (M.nstatic+1):(M.nstatic+M.npred+M.nboth) ) ); %state variables in declaration order - end - state_var_DR = oo.dr.inv_order_var(state_var); %state vector in DR order - minA = A(state_var_DR,state_var_DR); dminA = dA(state_var_DR,state_var_DR,:); - minB = B(state_var_DR,:); dminB = dB(state_var_DR,:,:); - minC = C(:,state_var_DR); dminC = dC(:,state_var_DR,:); - minD = D(:,:); dminD = dD(:,:,:); - [CheckCO,minnx,minA,minB,minC,minD,dminA,dminB,dminC,dminD] = get_minimal_state_representation(minA,minB,minC,minD,dminA,dminB,dminC,dminD); + % Derive and check minimal state vector of first-order + [CheckCO,minnx,minA,minB,minC,minD,dminA,dminB,dminC,dminD] = get_minimal_state_representation(oo.dr.ghx(oo.dr.pruned.indx,:),... %A + oo.dr.ghu(oo.dr.pruned.indx,:),... %B + oo.dr.ghx(oo.dr.pruned.indy,:),... %C + oo.dr.ghu(oo.dr.pruned.indy,:),... %D + oo.dr.derivs.dghx(oo.dr.pruned.indx,:,:),... %dA + oo.dr.derivs.dghu(oo.dr.pruned.indx,:,:),... %dB + oo.dr.derivs.dghx(oo.dr.pruned.indy,:,:),... %dC + oo.dr.derivs.dghu(oo.dr.pruned.indy,:,:)); %dD if CheckCO == 0 warning_KomunjerNg = 'WARNING: Komunjer and Ng (2011) failed:\n'; warning_KomunjerNg = [warning_KomunjerNg ' Conditions for minimality are not fullfilled:\n']; warning_KomunjerNg = [warning_KomunjerNg ' Skip identification analysis based on minimal state space system.\n']; fprintf(warning_KomunjerNg); %use sprintf to have line breaks - D = []; + dMINIMAL = []; else %reshape into Magnus-Neudecker Jacobians, i.e. dvec(X)/dp dminA = reshape(dminA,size(dminA,1)*size(dminA,2),size(dminA,3)); dminB = reshape(dminB,size(dminB,1)*size(dminB,2),size(dminB,3)); dminC = reshape(dminC,size(dminC,1)*size(dminC,2),size(dminC,3)); dminD = reshape(dminD,size(dminD,1)*size(dminD,2),size(dminD,3)); - dvechSig = reshape(dSig,size(dSig,1)*size(dSig,2),size(dSig,3)); + dvechSig = reshape(oo.dr.derivs.dSigma_e,exo_nbr*exo_nbr,totparam_nbr); indvechSig= find(tril(ones(exo_nbr,exo_nbr))); dvechSig = dvechSig(indvechSig,:); Inx = eye(minnx); @@ -387,64 +472,39 @@ if ~no_identification_minimal kron(Inu,minB); zeros(obs_nbr*minnx,exo_nbr^2); kron(Inu,minD); - -2*Enu*kron(M.Sigma_e,Inu)]; - D = full([KomunjerNg_DL KomunjerNg_DT KomunjerNg_DU]); - %add Jacobian of steady state - D = [[zeros(obs_nbr,stderrparam_nbr+corrparam_nbr) dYss(indvobs,:)], zeros(obs_nbr,minnx^2+exo_nbr^2); D]; + -2*Enu*kron(Sigma_e0,Inu)]; + dMINIMAL = full([KomunjerNg_DL KomunjerNg_DT KomunjerNg_DU]); + %add Jacobian of steady state (here we also allow for higher-order perturbation, i.e. only the mean provides additional restrictions + dMINIMAL = [dE_y zeros(obs_nbr,minnx^2+exo_nbr^2); dMINIMAL]; end end else - D = []; -end - -if nargout > 8 - options.ar = nlags; - nodecomposition = 1; - [Gamma_y,~] = th_autocovariances(oo.dr,oo.dr.order_var(indvobs),M,options,nodecomposition); %focus only on observed variables - sdy=sqrt(diag(Gamma_y{1})); %theoretical standard deviation - sy=sdy*sdy'; %cross products of standard deviations - if useautocorr - sy=sy-diag(diag(sy))+eye(obs_nbr); - Gamma_y{1}=Gamma_y{1}./sy; - else - for j=1:nlags - Gamma_y{j+1}=Gamma_y{j+1}.*sy; - end - end - %Ezz is vector of theoretical moments of VAROBS variables - MOMENTS = dyn_vech(Gamma_y{1}); - for j=1:nlags - MOMENTS = [MOMENTS; vec(Gamma_y{j+1})]; - end - MOMENTS = [oo.dr.ys(oo.dr.order_var(indvobs)); MOMENTS]; + dMINIMAL = []; end -function [dX] = DerivABCD(A,dA,B,dB,C,dC,D,dD) -% function [dX] = DerivABCD(A,dA,B,dB,C,dC,D,dD) +function [dX] = DerivABCD(X1,dX1,X2,dX2,X3,dX3,X4,dX4) +% function [dX] = DerivABCD(X1,dX1,X2,dX2,X3,dX3,X4,dX4) % ------------------------------------------------------------------------- -% Derivative of X(p)=A(p)*B(p)*C(p)*D(p) w.r.t to p +% Derivative of X(p)=X1(p)*X2(p)*X3(p)*X4(p) w.r.t to p % See Magnus and Neudecker (1999), p. 175 % ------------------------------------------------------------------------- -% Inputs: Matrices A,B,C,D, and the corresponding derivatives w.r.t p. -% Output: Derivative of product of A*B*C*D w.r.t. p +% Inputs: Matrices X1,X2,X3,X4, and the corresponding derivatives w.r.t p. +% Output: Derivative of product of X1*X2*X3*X4 w.r.t. p % ========================================================================= -nparam = size(dA,2); +nparam = size(dX1,2); % If one or more matrices are left out, they are set to zero if nargin == 4 - C=speye(size(B,2)); dC=spalloc(numel(C),nparam,0); - D=speye(size(C,2)); dD=spalloc(numel(D),nparam,0); + X3=speye(size(X2,2)); dX3=spalloc(numel(X3),nparam,0); + X4=speye(size(X3,2)); dX4=spalloc(numel(X4),nparam,0); elseif nargin == 6 - D=speye(size(C,2)); dD=spalloc(numel(D),nparam,0); + X4=speye(size(X3,2)); dX4=spalloc(numel(X4),nparam,0); end -dX1 = kron(transpose(D)*transpose(C)*transpose(B),speye(size(A,1)))*dA; -dX2 = kron(transpose(D)*transpose(C),A)*dB; -dX3 = kron(transpose(D),A*B)*dC; -dX4 = kron(speye(size(D,2)),A*B*C)*dD; -dX= dX1+dX2+dX3+dX4; +dX = kron(transpose(X4)*transpose(X3)*transpose(X2),speye(size(X1,1)))*dX1... + + kron(transpose(X4)*transpose(X3),X1)*dX2... + + kron(transpose(X4),X1*X2)*dX3... + + kron(speye(size(X4,2)),X1*X2*X3)*dX4; end %DerivABCD end -end%main function end - - +end%main function end \ No newline at end of file diff --git a/matlab/get_perturbation_params_derivs.m b/matlab/get_perturbation_params_derivs.m index 28b7fd119..d3efd45e4 100644 --- a/matlab/get_perturbation_params_derivs.m +++ b/matlab/get_perturbation_params_derivs.m @@ -260,8 +260,9 @@ if analytic_derivation_mode == -1 DERIVS.dghxss = reshape(dSig_gh(ind_ghxss,:), [M.endo_nbr M.nspred totparam_nbr]); %in tensor notation, wrt selected parameters DERIVS.dghuss = reshape(dSig_gh(ind_ghuss,:), [M.endo_nbr M.exo_nbr totparam_nbr]); %in tensor notation, wrt selected parameters end - % Parameter Jacobian of Om=ghu*Sigma_e*ghu' (wrt selected stderr, corr and model paramters) + % Parameter Jacobian of Om=ghu*Sigma_e*ghu' and Correlation_matrix (wrt selected stderr, corr and model paramters) DERIVS.dOm = zeros(M.endo_nbr,M.endo_nbr,totparam_nbr); %initialize in tensor notation + DERIVS.dCorrelation_matrix = zeros(M.exo_nbr,M.exo_nbr,totparam_nbr); %initialize if ~isempty(indpstderr) %derivatives of ghu wrt stderr parameters are zero by construction for jp=1:stderrparam_nbr DERIVS.dOm(:,:,jp) = oo.dr.ghu*DERIVS.dSigma_e(:,:,jp)*oo.dr.ghu'; @@ -270,6 +271,8 @@ if analytic_derivation_mode == -1 if ~isempty(indpcorr) %derivatives of ghu wrt corr parameters are zero by construction for jp=1:corrparam_nbr DERIVS.dOm(:,:,stderrparam_nbr+jp) = oo.dr.ghu*DERIVS.dSigma_e(:,:,stderrparam_nbr+jp)*oo.dr.ghu'; + DERIVS.dCorrelation_matrix(indpcorr(jp,1),indpcorr(jp,2),stderrparam_nbr+jp) = 1; + DERIVS.dCorrelation_matrix(indpcorr(jp,2),indpcorr(jp,1),stderrparam_nbr+jp) = 1;%symmetry end end if ~isempty(indpmodel) %derivatives of Sigma_e wrt model parameters are zero by construction @@ -1190,6 +1193,7 @@ end %% Store into structure DERIVS.dg1 = dg1; DERIVS.dSigma_e = dSigma_e; +DERIVS.dCorrelation_matrix = dCorrelation_matrix; DERIVS.dYss = dYss; DERIVS.dghx = cat(3,zeros(M.endo_nbr,M.nspred,stderrparam_nbr+corrparam_nbr), dghx); DERIVS.dghu = cat(3,zeros(M.endo_nbr,M.exo_nbr,stderrparam_nbr+corrparam_nbr), dghu); diff --git a/matlab/identification_analysis.m b/matlab/identification_analysis.m index 9f72aaafb..b7d326e2f 100644 --- a/matlab/identification_analysis.m +++ b/matlab/identification_analysis.m @@ -1,12 +1,15 @@ -function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide_lre, derivatives_info, info, options_ident] = identification_analysis(params, indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, init) -% [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide_lre, derivatives_info, info, options_ident] = identification_analysis(params, indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, init) +function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide_dynamic, derivatives_info, info, options_ident] = identification_analysis(params, indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, init) +% [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide_dynamic, derivatives_info, info, options_ident] = identification_analysis(params, indpmodel, indpstderr, indpcorr, options_ident, dataset_info, prior_exist, init) % ------------------------------------------------------------------------- % This function wraps all identification analysis, i.e. it % (1) wraps functions for the theoretical identification analysis based on moments (Iskrev, 2010), % spectrum (Qu and Tkachenko, 2012), minimal system (Komunjer and Ng, 2011), information matrix, -% reduced-form solution and linear rational expectation model (Ratto and Iskrev, 2011). +% reduced-form solution and dynamic model derivatives (Ratto and Iskrev, 2011). % (2) computes the identification strength based on moments (Ratto and Iskrev, 2011) % (3) checks which parameters are involved. +% If options_ident.order>1, then the identification analysis is based on +% Mutschler (2015), i.e. the pruned state space system and the corresponding +% moments, spectrum, reduced-form solution and dynamic model derivatives % ========================================================================= % INPUTS % * params [mc_sample_nbr by totparam_nbr] @@ -15,8 +18,8 @@ function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide % index of model parameters for which identification is checked % * indpstderr [stderrparam_nbr by 1] % index of stderr parameters for which identification is checked -% * indpcorr [corrparam_nbr by 1] -% index of corr parmeters for which identification is checked +% * indpcorr [corrparam_nbr by 2] +% matrix of corr parmeters for which identification is checked % * options_ident [structure] % identification options % * dataset_info [structure] @@ -25,23 +28,21 @@ function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide % 1: prior exists. Identification is checked for stderr, corr and model parameters as declared in estimated_params block % 0: prior does not exist. Identification is checked for all stderr and model parameters, but no corr parameters % * init [integer] -% flag for initialization of persistent vars. This is needed if one may want to re-initialize persistent -% variables even if they are not empty, e.g. by making more calls to identification in the same dynare mod file +% flag for initialization of persistent vars. This is needed if one may wants to make more calls to identification in the same dynare mod file % ------------------------------------------------------------------------- % OUTPUTS % * ide_moments [structure] -% identification results using theoretical moments (Iskrev, 2010) +% identification results using theoretical moments (Iskrev, 2010; Mutschler, 2015) % * ide_spectrum [structure] -% identification results using spectrum (Qu and Tkachenko, 2012) +% identification results using spectrum (Qu and Tkachenko, 2012; Mutschler, 2015) % * ide_minimal [structure] % identification results using minimal system (Komunjer and Ng, 2011) % * ide_hess [structure] % identification results using asymptotic Hessian (Ratto and Iskrev, 2011) % * ide_reducedform [structure] -% identification results using reduced form solution (Ratto and Iskrev, 2011) -% * ide_lre [structure] -% identification results using linear rational expectations model, -% i.e. steady state and Jacobian of dynamic model, (Ratto and Iskrev, 2011) +% identification results using steady state and reduced form solution (Ratto and Iskrev, 2011) +% * ide_dynamic [structure] +% identification results using steady state and dynamic model derivatives (Ratto and Iskrev, 2011) % * derivatives_info [structure] % info about derivatives, used in dsge_likelihood.m % * info [integer] @@ -57,7 +58,6 @@ function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide % * dseries % * dsge_likelihood.m % * dyn_vech -% * dynare_resolve % * ident_bruteforce % * identification_checks % * identification_checks_via_subsets @@ -65,6 +65,7 @@ function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide % * get_identification_jacobians (previously getJJ) % * matlab_ver_less_than % * prior_bounds +% * resol % * set_all_parameters % * simulated_moment_uncertainty % * stoch_simul @@ -89,23 +90,23 @@ function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide % ========================================================================= global oo_ M_ options_ bayestopt_ estim_params_ -persistent ind_J ind_G ind_D ind_dTAU ind_dLRE +persistent ind_dMOMENTS ind_dREDUCEDFORM ind_dDYNAMIC % persistence is necessary, because in a MC loop the numerical threshold used may provide vectors of different length, leading to crashes in MC loops %initialize output structures -ide_hess = struct(); %asymptotic/simulated information matrix -ide_reducedform = struct(); %reduced form solution -ide_lre = struct(); %linear rational expectation model -ide_moments = struct(); %Iskrev (2010)'s J -ide_spectrum = struct(); %Qu and Tkachenko (2012)'s G -ide_minimal = struct(); %Komunjer and Ng (2011)'s D +ide_hess = struct(); %Jacobian of asymptotic/simulated information matrix +ide_reducedform = struct(); %Jacobian of steady state and reduced form solution +ide_dynamic = struct(); %Jacobian of steady state and dynamic model derivatives +ide_moments = struct(); %Jacobian of first two moments (Iskrev, 2010; Mutschler, 2015) +ide_spectrum = struct(); %Gram matrix of Jacobian of spectral density plus Gram matrix of Jacobian of steady state (Qu and Tkachenko, 2012; Mutschler, 2015) +ide_minimal = struct(); %Jacobian of minimal system (Komunjer and Ng, 2011) derivatives_info = struct(); %storage for Jacobians used in dsge_likelihood.m for (analytical or simulated) asymptotic Gradient and Hession of likelihood totparam_nbr = length(params); %number of all parameters to be checked modparam_nbr = length(indpmodel); %number of model parameters to be checked stderrparam_nbr = length(indpstderr); %number of stderr parameters to be checked corrparam_nbr = size(indpcorr,1); %number of stderr parameters to be checked -indvobs = bayestopt_.mf2; % index of observable variables +indvobs = bayestopt_.mf2; %index of observable variables if ~isempty(estim_params_) %estimated_params block is available, so we are able to use set_all_parameters.m @@ -113,163 +114,138 @@ if ~isempty(estim_params_) end %get options (see dynare_identification.m for description of options) -nlags = options_ident.ar; -advanced = options_ident.advanced; -replic = options_ident.replic; -periods = options_ident.periods; -max_dim_cova_group = options_ident.max_dim_cova_group; -normalize_jacobians = options_ident.normalize_jacobians; -tol_deriv = options_ident.tol_deriv; -tol_rank = options_ident.tol_rank; -tol_sv = options_ident.tol_sv; -no_identification_strength = options_ident.no_identification_strength; -no_identification_reducedform = options_ident.no_identification_reducedform; -no_identification_moments = options_ident.no_identification_moments; -no_identification_minimal = options_ident.no_identification_minimal; -no_identification_spectrum = options_ident.no_identification_spectrum; -checks_via_subsets = options_ident.checks_via_subsets; +no_identification_strength = options_ident.no_identification_strength; +no_identification_reducedform = options_ident.no_identification_reducedform; +no_identification_moments = options_ident.no_identification_moments; +no_identification_minimal = options_ident.no_identification_minimal; +no_identification_spectrum = options_ident.no_identification_spectrum; +order = options_ident.order; +nlags = options_ident.ar; +advanced = options_ident.advanced; +replic = options_ident.replic; +periods = options_ident.periods; +max_dim_cova_group = options_ident.max_dim_cova_group; +normalize_jacobians = options_ident.normalize_jacobians; +checks_via_subsets = options_ident.checks_via_subsets; +tol_deriv = options_ident.tol_deriv; +tol_rank = options_ident.tol_rank; +tol_sv = options_ident.tol_sv; -[I,~] = find(M_.lead_lag_incidence'); %I is used to select nonzero columns of the Jacobian of endogenous variables in dynamic model files - -%Compute linear approximation and matrices A and B of the transition equation for all endogenous variables in DR order -[A, B, ~, info, M_, options_, oo_] = dynare_resolve(M_,options_,oo_); +%Compute linear approximation and fill dr structure +[oo_.dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_); if info(1) == 0 %no errors in solution - oo0 = oo_; %store oo_ structure - y0 = oo_.dr.ys(oo_.dr.order_var); %steady state in DR order - yy0 = oo_.dr.ys(I); %steady state of dynamic (endogenous and auxiliary variables) in DR order - ex0 = oo_.exo_steady_state'; %steady state of exogenous variables in declaration order - param0 = M_.params; %parameter values at which to evaluate dynamic, static and param_derivs files - Sigma_e0 = M_.Sigma_e; %store covariance matrix of exogenous shocks - - %compute vectorized reduced-form solution in DR order - tau = [y0; vec(A); dyn_vech(B*Sigma_e0*B')]; - - %compute vectorized linear rational expectations model in DR order - [~, g1 ] = feval([M_.fname,'.dynamic'], yy0, repmat(ex0, M_.maximum_exo_lag+M_.maximum_exo_lead+1), param0, oo_.dr.ys, 1); - %g1 is [endo_nbr by (dynamicvar_nbr+exo_nbr)] first derivative (wrt all endogenous, exogenous and auxiliary variables) of dynamic model equations, i.e. df/d[yy0;ex0], in DR order - lre = [y0; vec(g1)]; %add steady state in DR order and vectorize - - %Compute Jacobians for identification analysis either analytically or numerically dependent on analytic_derivation_mode in options_ident - [J, G, D, dTAU, dLRE, dA, dOm, dYss, MOMENTS] = get_identification_jacobians(A, B, estim_params_, M_, oo0, options_, options_ident, indpmodel, indpstderr, indpcorr, indvobs); - if isempty(D) + [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); + 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; options_ident.no_identification_minimal = 1; end - %store derivatives for use in dsge_likelihood.m - derivatives_info.DT = dA; - derivatives_info.DOm = dOm; - derivatives_info.DYss = dYss; - if init %check stationarity if ~no_identification_moments - ind_J = (find(max(abs(J'),[],1) > tol_deriv)); %index for non-zero derivatives - if isempty(ind_J) && any(any(isnan(J))) - error('There are NaN in the JJ matrix. Please check whether your model has units roots and you forgot to set diffuse_filter=1.' ) + ind_dMOMENTS = (find(max(abs(dMOMENTS'),[],1) > tol_deriv)); %index for non-zero rows + if isempty(ind_dMOMENTS) && any(any(isnan(dMOMENTS))) + error('There are NaN in the dMOMENTS matrix. Please check whether your model has units roots and you forgot to set diffuse_filter=1.' ) end if any(any(isnan(MOMENTS))) error('There are NaN''s in the theoretical moments: make sure that for non-stationary models stationary transformations of non-stationary observables are used for checking identification. [TIP: use first differences].') end end if ~no_identification_spectrum - ind_G = (find(max(abs(G'),[],1) > tol_deriv)); %index for non-zero derivatives - if isempty(ind_G) && any(any(isnan(G))) - warning_QuTkachenko = 'WARNING: There are NaN in Qu and Tkachenko (2012)''s G matrix. Please check whether your model has units roots (diffuse_filter=1 does not work yet for the G matrix).\n'; - warning_QuTkachenko = [warning_QuTkachenko ' Skip identification analysis based on spectrum.\n']; - fprintf(warning_QuTkachenko); - %reset options to neither display nor plot Qu and Tkachenko's G anymore + ind_dSPECTRUM = (find(max(abs(dSPECTRUM'),[],1) > tol_deriv)); %index for non-zero rows + if isempty(ind_dSPECTRUM) && any(any(isnan(dSPECTRUM))) + warning_SPECTRUM = 'WARNING: There are NaN in the dSPECTRUM matrix. Please check whether your model has units roots and your forgot to set diffuse_filter=1.\n'; + warning_SPECTRUM = [warning_SPECTRUM ' Skip identification analysis based on spectrum.\n']; + fprintf(warning_SPECTRUM); + %reset options to neither display nor plot dSPECTRUM anymore no_identification_spectrum = 1; options_ident.no_identification_spectrum = 1; end end if ~no_identification_minimal - ind_D = (find(max(abs(D'),[],1) > tol_deriv)); %index for non-zero derivatives - if isempty(ind_D) && any(any(isnan(D))) - warning_KomunjerNg = 'WARNING: There are NaN in Komunjer and Ng (2011)''s D matrix. Please check whether your model has units roots and you forgot to set diffuse_filter=1.\n'; - warning_KomunjerNg = [warning_KomunjerNg ' Skip identification analysis based on minimal system.\n']; - fprintf(warning_KomunjerNg); - %reset options to neither display nor plot Komunjer and Ng's D anymore + ind_dMINIMAL = (find(max(abs(dMINIMAL'),[],1) > tol_deriv)); %index for non-zero rows + if isempty(ind_dMINIMAL) && any(any(isnan(dMINIMAL))) + warning_MINIMAL = 'WARNING: There are NaN in the dMINIMAL matrix. Please check whether your model has units roots and you forgot to set diffuse_filter=1.\n'; + warning_MINIMAL = [warning_MINIMAL ' Skip identification analysis based on minimal system.\n']; + fprintf(warning_MINIMAL); + %reset options to neither display nor plot dMINIMAL anymore no_identification_minimal = 1; options_ident.no_identification_minimal = 1; end end if no_identification_moments && no_identification_minimal && no_identification_spectrum - %error if all three criteria fail + %display error if all three criteria fail error('identification_analyis: Stationarity condition(s) failed and/or diffuse_filter option missing'); end % Check order conditions if ~no_identification_moments %check order condition of Iskrev (2010) - while length(ind_J) < totparam_nbr && nlags < 10 + while length(ind_dMOMENTS) < totparam_nbr && nlags < 10 %Try to add lags to autocovariogram if order condition fails disp('The number of moments with non-zero derivative is smaller than the number of parameters') disp(['Try increasing ar = ', int2str(nlags+1)]) nlags = nlags + 1; - options_ident.no_identification_minimal = 1; %do not recompute D - options_ident.no_identification_spectrum = 1; %do not recompute G + options_ident.no_identification_minimal = 1; %do not recompute dMINIMAL + 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 - [J, ~, ~, dTAU, dLRE, dA, dOm, dYss, MOMENTS] = get_identification_jacobians(A, B, estim_params_, M_, oo0, options_, options_ident, indpmodel, indpstderr, indpcorr, indvobs); - ind_J = (find(max(abs(J'),[],1) > tol_deriv)); %new index with non-zero derivatives - %store derivatives for use in dsge_likelihood.m - derivatives_info.DT = dA; - derivatives_info.DOm = dOm; - derivatives_info.DYss = dYss; + [~, ~, ~, ~, ~, ~, 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 - options_ident.no_identification_minimal = no_identification_minimal; % reset option to original choice - options_ident.no_identification_spectrum = no_identification_spectrum; % reset option to original choice - if length(ind_J) < totparam_nbr - warning_Iskrev = 'WARNING: Order condition for Iskrev (2010) failed: There are not enough moments and too many parameters.\n'; - warning_Iskrev = [warning_Iskrev ' The number of moments with non-zero derivative is smaller than the number of parameters up to 10 lags.\n']; - warning_Iskrev = [warning_Iskrev ' Either reduce the list of parameters, or further increase ar, or increase number of observables.\n']; - warning_Iskrev = [warning_Iskrev ' Skip identification analysis based on moments.\n']; - fprintf(warning_Iskrev); - %reset options to neither display nor plot Iskrev's J anymore + options_ident.no_identification_minimal = no_identification_minimal; % reset option to original setting + options_ident.no_identification_spectrum = no_identification_spectrum; % reset option to original setting + if length(ind_dMOMENTS) < totparam_nbr + warning_MOMENTS = 'WARNING: Order condition for dMOMENTS failed: There are not enough moments and too many parameters.\n'; + warning_MOMENTS = [warning_MOMENTS ' The number of moments with non-zero derivative is smaller than the number of parameters up to 10 lags.\n']; + warning_MOMENTS = [warning_MOMENTS ' Either reduce the list of parameters, or further increase ar, or increase number of observables.\n']; + warning_MOMENTS = [warning_MOMENTS ' Skip identification analysis based on moments.\n']; + fprintf(warning_MOMENTS); + %reset options to neither display nor plot dMOMENTS anymore no_identification_moments = 1; options_ident.no_identification_moments = 1; end end if ~no_identification_minimal - if length(ind_D) < size(D,2) - warning_KomunjerNg = 'WARNING: Order condition for Komunjer and Ng (2011) failed: There are too many parameters or too few observable variables.\n'; - warning_KomunjerNg = [warning_KomunjerNg ' The number of minimal system elements with non-zero derivative is smaller than the number of parameters.\n']; - warning_KomunjerNg = [warning_KomunjerNg ' Either reduce the list of parameters, or increase number of observables.\n']; - warning_KomunjerNg = [warning_KomunjerNg ' Skip identification analysis based on minimal state space system.\n']; - fprintf(warning_KomunjerNg); - %resest options to neither display nor plot Komunjer and Ng's D anymore + if length(ind_dMINIMAL) < size(dMINIMAL,2) + warning_MINIMAL = 'WARNING: Order condition for dMINIMAL failed: There are too many parameters or too few observable variables.\n'; + warning_MINIMAL = [warning_MINIMAL ' The number of minimal system elements with non-zero derivative is smaller than the number of parameters.\n']; + warning_MINIMAL = [warning_MINIMAL ' Either reduce the list of parameters, or increase number of observables.\n']; + warning_MINIMAL = [warning_MINIMAL ' Skip identification analysis based on minimal state space system.\n']; + fprintf(warning_MINIMAL); + %resest options to neither display nor plot dMINIMAL anymore no_identification_minimal = 1; options_ident.no_identification_minimal = 1; end end - %There is no order condition for Qu and Tkachenko's G + %Note that there is no order condition for dSPECTRUM, as the matrix is always of dimension totparam_nbr by totparam_nbr if no_identification_moments && no_identification_minimal && no_identification_spectrum %error if all three criteria fail error('identification_analyis: Order condition(s) failed'); end if ~no_identification_reducedform - ind_dTAU = (find(max(abs(dTAU'),[],1) > tol_deriv)); %index with non-zero derivatives + ind_dREDUCEDFORM = (find(max(abs(dREDUCEDFORM'),[],1) > tol_deriv)); %index with non-zero rows end - ind_dLRE = (find(max(abs(dLRE'),[],1) > tol_deriv)); %index with non-zero derivatives + ind_dDYNAMIC = (find(max(abs(dDYNAMIC'),[],1) > tol_deriv)); %index with non-zero rows end - LRE = lre(ind_dLRE); - si_dLRE = (dLRE(ind_dLRE,:)); %focus only on non-zero derivatives + DYNAMIC = DYNAMIC(ind_dDYNAMIC); %focus only on non-zero entries + si_dDYNAMIC = (dDYNAMIC(ind_dDYNAMIC,:)); %focus only on non-zero rows if ~no_identification_reducedform - TAU = tau(ind_dTAU); - si_dTAU = (dTAU(ind_dTAU,:)); %focus only on non-zero derivatives + REDUCEDFORM = REDUCEDFORM(ind_dREDUCEDFORM); %focus only on non-zero entries + si_dREDUCEDFORM = (dREDUCEDFORM(ind_dREDUCEDFORM,:)); %focus only on non-zero rows end if ~no_identification_moments - MOMENTS = MOMENTS(ind_J); - si_J = (J(ind_J,:)); %focus only on non-zero derivatives - %% compute identification strength based on moments + MOMENTS = MOMENTS(ind_dMOMENTS); %focus only on non-zero entries + si_dMOMENTS = (dMOMENTS(ind_dMOMENTS,:)); %focus only on non-zero derivatives +%% MOMENTS IDENTIFICATION STRENGTH ANALYSIS if ~no_identification_strength && init %only for initialization of persistent vars - ide_strength_J = NaN(1,totparam_nbr); %initialize - ide_strength_J_prior = NaN(1,totparam_nbr); %initialize + ide_strength_dMOMENTS = NaN(1,totparam_nbr); %initialize + ide_strength_dMOMENTS_prior = NaN(1,totparam_nbr); %initialize ide_uncert_unnormaliz = NaN(1,totparam_nbr); %initialize if prior_exist offset_prior = 0; @@ -291,11 +267,10 @@ if info(1) == 0 %no errors in solution end try - %try to compute asymptotic Hessian for identification strength analysis based on moments + %try to compute asymptotic Hessian for identification strength analysis based on moments % reset some options for faster computations options_.irf = 0; options_.noprint = 1; - options_.order = 1; options_.SpectralDensity.trigger = 0; options_.periods = periods+100; if options_.kalman_algo > 2 @@ -306,17 +281,17 @@ if info(1) == 0 %no errors in solution [info, oo_, options_] = stoch_simul(M_, options_, oo_, options_.varobs); dataset_ = dseries(oo_.endo_simul(options_.varobs_id,100+1:end)',dates('1Q1'), options_.varobs); %get information on moments derivatives_info.no_DLIK = 1; - bounds = prior_bounds(bayestopt_, options_.prior_trunc); %reset bounds as lb and ub must only be operational during mode-finding + bounds = prior_bounds(bayestopt_, options_.prior_trunc); %reset bounds as lb and ub must only be operational during mode-finding + if options_.order > 1 + fprintf('IDENTIFICATION STRENGTH: Analytic computation of Hessian is not available for ''order>1''\n') + fprintf(' Identification strength is based on simulated moment uncertainty\n'); + end + %note that for order>1 we do not provide any information on DT,DYss,DOM in derivatives_info, such that dsge_likelihood creates an error. Therefore the computation will be based on simulated_moment_uncertainty for order>1. [fval, info, cost_flag, DLIK, AHess, ys, trend_coeff, M_, options_, bayestopt_, oo_] = dsge_likelihood(params', dataset_, dataset_info, options_, M_, estim_params_, bayestopt_, bounds, oo_, derivatives_info); %non-used output variables need to be set for octave for some reason %note that for the order of parameters in AHess we have: stderr parameters come first, corr parameters second, model parameters third. the order within these blocks corresponds to the order specified in the estimated_params block options_.analytic_derivation = analytic_derivation; %reset option AHess = -AHess; %take negative of hessian - if ischar(tol_rank) %if tolerance is specified as 'robust' - tol_rankAhess = max(size(AHess))*eps(norm(AHess)); - else - tol_rankAhess = tol_rank; - end - if min(eig(AHess))<-tol_rankAhess + if min(eig(AHess))<-tol_rank error('identification_analysis: Analytic Hessian is not positive semi-definite!') end ide_hess.AHess = AHess; %store asymptotic Hessian @@ -332,16 +307,16 @@ if info(1) == 0 %no errors in solution indok = find(max(ide_hess.indno,[],1)==0); ide_uncert_unnormaliz(indok) = sqrt(diag(inv(AHess(indok,indok))))'; ind1 = find(ide_hess.ind0); - cmm = si_J(:,ind1)*((AHess(ind1,ind1))\si_J(:,ind1)'); %covariance matrix of moments - temp1 = ((AHess(ind1,ind1))\si_dTAU(:,ind1)'); - diag_chh = sum(si_dTAU(:,ind1)'.*temp1)'; + cmm = si_dMOMENTS(:,ind1)*((AHess(ind1,ind1))\si_dMOMENTS(:,ind1)'); %covariance matrix of moments + temp1 = ((AHess(ind1,ind1))\si_dREDUCEDFORM(:,ind1)'); + diag_chh = sum(si_dREDUCEDFORM(:,ind1)'.*temp1)'; ind1 = ind1(ind1>stderrparam_nbr+corrparam_nbr); - clre = si_dLRE(:,ind1-stderrparam_nbr-corrparam_nbr)*((AHess(ind1,ind1))\si_dLRE(:,ind1-stderrparam_nbr-corrparam_nbr)'); - flag_score = 1; %this is only used for a title of a plot in plot_identification.m + cdynamic = si_dDYNAMIC(:,ind1-stderrparam_nbr-corrparam_nbr)*((AHess(ind1,ind1))\si_dDYNAMIC(:,ind1-stderrparam_nbr-corrparam_nbr)'); + flag_score = 1; %this is used for the title in plot_identification.m catch %Asymptotic Hessian via simulation - replic = max([replic, length(ind_J)*3]); - cmm = simulated_moment_uncertainty(ind_J, periods, replic,options_,M_,oo_); %covariance matrix of moments + replic = max([replic, length(ind_dMOMENTS)*3]); + cmm = simulated_moment_uncertainty(ind_dMOMENTS, periods, replic,options_,M_,oo_); %covariance matrix of moments sd = sqrt(diag(cmm)); cc = cmm./(sd*sd'); if isoctave || matlab_ver_less_than('8.3') @@ -353,7 +328,7 @@ if info(1) == 0 %no errors in solution [VV,DD,WW] = eig(cc); end id = find(diag(DD)>tol_deriv); - siTMP = si_J./repmat(sd,[1 totparam_nbr]); + siTMP = si_dMOMENTS./repmat(sd,[1 totparam_nbr]); MIM = (siTMP'*VV(:,id))*(DD(id,id)\(WW(:,id)'*siTMP)); clear siTMP; ide_hess.AHess = MIM; %store asymptotic hessian @@ -368,122 +343,123 @@ if info(1) == 0 %no errors in solution end indok = find(max(ide_hess.indno,[],1)==0); ind1 = find(ide_hess.ind0); - temp1 = ((MIM(ind1,ind1))\si_dTAU(:,ind1)'); - diag_chh = sum(si_dTAU(:,ind1)'.*temp1)'; + temp1 = ((MIM(ind1,ind1))\si_dREDUCEDFORM(:,ind1)'); + diag_chh = sum(si_dREDUCEDFORM(:,ind1)'.*temp1)'; ind1 = ind1(ind1>stderrparam_nbr+corrparam_nbr); - clre = si_dLRE(:,ind1-stderrparam_nbr-corrparam_nbr)*((MIM(ind1,ind1))\si_dLRE(:,ind1-stderrparam_nbr-corrparam_nbr)'); + cdynamic = si_dDYNAMIC(:,ind1-stderrparam_nbr-corrparam_nbr)*((MIM(ind1,ind1))\si_dDYNAMIC(:,ind1-stderrparam_nbr-corrparam_nbr)'); if ~isempty(indok) ide_uncert_unnormaliz(indok) = (sqrt(diag(inv(tildaM(indok,indok))))./deltaM(indok))'; %sqrt(diag(inv(MIM(indok,indok))))'; end - flag_score = 0; %this is only used for a title of a plot + flag_score = 0; %this is used for the title in plot_identification.m end % end of computing sample information matrix for identification strength measure - ide_strength_J(indok) = (1./(ide_uncert_unnormaliz(indok)'./abs(params(indok)'))); %this is s_i in Ratto and Iskrev (2011, p.13) - ide_strength_J_prior(indok) = (1./(ide_uncert_unnormaliz(indok)'./normaliz_prior_std(indok)')); %this is s_i^{prior} in Ratto and Iskrev (2011, p.14) + ide_strength_dMOMENTS(indok) = (1./(ide_uncert_unnormaliz(indok)'./abs(params(indok)'))); %this is s_i in Ratto and Iskrev (2011, p.13) + ide_strength_dMOMENTS_prior(indok) = (1./(ide_uncert_unnormaliz(indok)'./normaliz_prior_std(indok)')); %this is s_i^{prior} in Ratto and Iskrev (2011, p.14) sensitivity_zero_pos = find(isinf(deltaM)); deltaM_prior = deltaM.*abs(normaliz_prior_std'); %this is \Delta_i^{prior} in Ratto and Iskrev (2011, p.14) deltaM = deltaM.*abs(params'); %this is \Delta_i in Ratto and Iskrev (2011, p.14) - quant = si_J./repmat(sqrt(diag(cmm)),1,totparam_nbr); + quant = si_dMOMENTS./repmat(sqrt(diag(cmm)),1,totparam_nbr); if size(quant,1)==1 - si_Jnorm = abs(quant).*normaliz_prior_std; + si_dMOMENTSnorm = abs(quant).*normaliz_prior_std; else - si_Jnorm = vnorm(quant).*normaliz_prior_std; + si_dMOMENTSnorm = vnorm(quant).*normaliz_prior_std; end iy = find(diag_chh); - ind_dTAU = ind_dTAU(iy); - si_dTAU = si_dTAU(iy,:); + ind_dREDUCEDFORM = ind_dREDUCEDFORM(iy); + si_dREDUCEDFORM = si_dREDUCEDFORM(iy,:); if ~isempty(iy) - quant = si_dTAU./repmat(sqrt(diag_chh(iy)),1,totparam_nbr); + quant = si_dREDUCEDFORM./repmat(sqrt(diag_chh(iy)),1,totparam_nbr); if size(quant,1)==1 - si_dTAUnorm = abs(quant).*normaliz_prior_std; + si_dREDUCEDFORMnorm = abs(quant).*normaliz_prior_std; else - si_dTAUnorm = vnorm(quant).*normaliz_prior_std; + si_dREDUCEDFORMnorm = vnorm(quant).*normaliz_prior_std; end else - si_dTAUnorm = []; + si_dREDUCEDFORMnorm = []; end - diag_clre = diag(clre); - iy = find(diag_clre); - ind_dLRE = ind_dLRE(iy); - si_dLRE = si_dLRE(iy,:); + diag_cdynamic = diag(cdynamic); + iy = find(diag_cdynamic); + ind_dDYNAMIC = ind_dDYNAMIC(iy); + si_dDYNAMIC = si_dDYNAMIC(iy,:); if ~isempty(iy) - quant = si_dLRE./repmat(sqrt(diag_clre(iy)),1,modparam_nbr); + quant = si_dDYNAMIC./repmat(sqrt(diag_cdynamic(iy)),1,modparam_nbr); if size(quant,1)==1 - si_dLREnorm = abs(quant).*normaliz_prior_std(stderrparam_nbr+corrparam_nbr+1:end); + si_dDYNAMICnorm = abs(quant).*normaliz_prior_std(stderrparam_nbr+corrparam_nbr+1:end); else - si_dLREnorm = vnorm(quant).*normaliz_prior_std(stderrparam_nbr+corrparam_nbr+1:end); + si_dDYNAMICnorm = vnorm(quant).*normaliz_prior_std(stderrparam_nbr+corrparam_nbr+1:end); end else - si_dLREnorm=[]; + si_dDYNAMICnorm=[]; end %store results of identification strength - ide_hess.ide_strength_J = ide_strength_J; - ide_hess.ide_strength_J_prior = ide_strength_J_prior; + ide_hess.ide_strength_dMOMENTS = ide_strength_dMOMENTS; + ide_hess.ide_strength_dMOMENTS_prior = ide_strength_dMOMENTS_prior; ide_hess.deltaM = deltaM; ide_hess.deltaM_prior = deltaM_prior; ide_hess.sensitivity_zero_pos = sensitivity_zero_pos; ide_hess.identified_parameter_indices = indok; ide_hess.flag_score = flag_score; - ide_lre.si_dLREnorm = si_dLREnorm; - ide_moments.si_Jnorm = si_Jnorm; - ide_reducedform.si_dTAUnorm = si_dTAUnorm; + ide_dynamic.si_dDYNAMICnorm = si_dDYNAMICnorm; + ide_moments.si_dMOMENTSnorm = si_dMOMENTSnorm; + ide_reducedform.si_dREDUCEDFORMnorm = si_dREDUCEDFORMnorm; end %end of identification strength analysis end - %% Theoretical identification analysis based on linear rational expectations model, i.e. steady state and Jacobian of dynamic model equations +%% Normalization of Jacobians +% For Dynamic, ReducedForm, Moment and Minimal Jacobian: rescale each row by its largest element in absolute value +% For Spectrum: transform into correlation-type matrices (as above with AHess) if normalize_jacobians - norm_dLRE = max(abs(si_dLRE),[],2); - norm_dLRE = norm_dLRE(:,ones(size(dLRE,2),1)); + norm_dDYNAMIC = max(abs(si_dDYNAMIC),[],2); + norm_dDYNAMIC = norm_dDYNAMIC(:,ones(size(dDYNAMIC,2),1)); else - norm_dLRE = 1; + norm_dDYNAMIC = 1; end - % store into structure (not everything is really needed) - ide_lre.ind_dLRE = ind_dLRE; - ide_lre.norm_dLRE = norm_dLRE; - ide_lre.si_dLRE = si_dLRE; - ide_lre.dLRE = dLRE; - ide_lre.LRE = LRE; + % store into structure (not everything is used later on) + ide_dynamic.ind_dDYNAMIC = ind_dDYNAMIC; + ide_dynamic.norm_dDYNAMIC = norm_dDYNAMIC; + ide_dynamic.si_dDYNAMIC = si_dDYNAMIC; + ide_dynamic.dDYNAMIC = dDYNAMIC; + ide_dynamic.DYNAMIC = DYNAMIC; - %% Theoretical identification analysis based on steady state and reduced-form-solution if ~no_identification_reducedform if normalize_jacobians - norm_dTAU = max(abs(si_dTAU),[],2); - norm_dTAU = norm_dTAU(:,ones(totparam_nbr,1)); + norm_dREDUCEDFORM = max(abs(si_dREDUCEDFORM),[],2); + norm_dREDUCEDFORM = norm_dREDUCEDFORM(:,ones(totparam_nbr,1)); else - norm_dTAU = 1; + norm_dREDUCEDFORM = 1; end - % store into structure (not everything is really needed) - ide_reducedform.ind_dTAU = ind_dTAU; - ide_reducedform.norm_dTAU = norm_dTAU; - ide_reducedform.si_dTAU = si_dTAU; - ide_reducedform.dTAU = dTAU; - ide_reducedform.TAU = TAU; + % store into structure (not everything is used later on) + ide_reducedform.ind_dREDUCEDFORM = ind_dREDUCEDFORM; + ide_reducedform.norm_dREDUCEDFORM = norm_dREDUCEDFORM; + ide_reducedform.si_dREDUCEDFORM = si_dREDUCEDFORM; + ide_reducedform.dREDUCEDFORM = dREDUCEDFORM; + ide_reducedform.REDUCEDFORM = REDUCEDFORM; end - %% Theoretical identification analysis based on Iskrev (2010)'s method, i.e. first two moments if ~no_identification_moments if normalize_jacobians - norm_J = max(abs(si_J),[],2); - norm_J = norm_J(:,ones(totparam_nbr,1)); + norm_dMOMENTS = max(abs(si_dMOMENTS),[],2); + norm_dMOMENTS = norm_dMOMENTS(:,ones(totparam_nbr,1)); else - norm_J = 1; + norm_dMOMENTS = 1; end - % store into structure (not everything is really needed) - ide_moments.ind_J = ind_J; - ide_moments.norm_J = norm_J; - ide_moments.si_J = si_J; - ide_moments.J = J; - ide_moments.MOMENTS = MOMENTS; + % store into structure (not everything is used later on) + ide_moments.ind_dMOMENTS = ind_dMOMENTS; + ide_moments.norm_dMOMENTS = norm_dMOMENTS; + ide_moments.si_dMOMENTS = si_dMOMENTS; + ide_moments.dMOMENTS = dMOMENTS; + ide_moments.MOMENTS = MOMENTS; if advanced - % here we do not normalize (i.e. set normJ=1) as the OLS in ident_bruteforce is very sensitive to normJ - [ide_moments.pars, ide_moments.cosnJ] = ident_bruteforce(J(ind_J,:), max_dim_cova_group, options_.TeX, options_ident.name_tex, options_ident.tittxt, options_ident.tol_deriv); + % here we do not normalize (i.e. we set norm_dMOMENTS=1) as the OLS in ident_bruteforce is very sensitive to norm_dMOMENTS + [ide_moments.pars, ide_moments.cosndMOMENTS] = ident_bruteforce(dMOMENTS(ind_dMOMENTS,:), max_dim_cova_group, options_.TeX, options_ident.name_tex, options_ident.tittxt, options_ident.tol_deriv); end - %here idea is to have the unnormalized S and V, which is then used in plot_identification.m and for prior_mc > 1 - [~, S, V] = svd(J(ind_J,:),0); + + %here we focus on the unnormalized S and V, which is then used in plot_identification.m and for prior_mc > 1 + [~, S, V] = svd(dMOMENTS(ind_dMOMENTS,:),0); S = diag(S); - S = [S;zeros(size(J,2)-length(ind_J),1)]; + S = [S;zeros(size(dMOMENTS,2)-length(ind_dMOMENTS),1)]; if totparam_nbr > 8 ide_moments.S = S([1:4, end-3:end]); ide_moments.V = V(:,[1:4, end-3:end]); @@ -493,62 +469,62 @@ if info(1) == 0 %no errors in solution end end - %% Theoretical identification analysis based on Komunjer and Ng (2011)'s method, i.e. steady state and observational equivalent spectral densities within a minimal system if ~no_identification_minimal if normalize_jacobians - norm_D = max(abs(D(ind_D,:)),[],2); - norm_D = norm_D(:,ones(size(D,2),1)); + ind_dMINIMAL = (find(max(abs(dMINIMAL'),[],1) > tol_deriv)); %index for non-zero rows + norm_dMINIMAL = max(abs(dMINIMAL(ind_dMINIMAL,:)),[],2); + norm_dMINIMAL = norm_dMINIMAL(:,ones(size(dMINIMAL,2),1)); else - norm_D = 1; + norm_dMINIMAL = 1; end - % store into structure - ide_minimal.ind_D = ind_D; - ide_minimal.norm_D = norm_D; - ide_minimal.D = D; + % store into structure (not everything is used later on) + ide_minimal.ind_dMINIMAL = ind_dMINIMAL; + ide_minimal.norm_dMINIMAL = norm_dMINIMAL; + ide_minimal.dMINIMAL = dMINIMAL; end - %% Theoretical identification analysis based on Qu and Tkachenko (2012)'s method, i.e. steady state and spectral density if ~no_identification_spectrum if normalize_jacobians - tilda_G = zeros(size(G)); - delta_G = sqrt(diag(G(ind_G,ind_G))); - tilda_G(ind_G,ind_G) = G(ind_G,ind_G)./((delta_G)*(delta_G')); - norm_G = max(abs(G(ind_G,:)),[],2); - norm_G = norm_G(:,ones(size(G,2),1)); + ind_dSPECTRUM = (find(max(abs(dSPECTRUM'),[],1) > tol_deriv)); %index for non-zero rows + tilda_dSPECTRUM = zeros(size(dSPECTRUM)); + delta_dSPECTRUM = sqrt(diag(dSPECTRUM(ind_dSPECTRUM,ind_dSPECTRUM))); + tilda_dSPECTRUM(ind_dSPECTRUM,ind_dSPECTRUM) = dSPECTRUM(ind_dSPECTRUM,ind_dSPECTRUM)./((delta_dSPECTRUM)*(delta_dSPECTRUM')); + norm_dSPECTRUM = max(abs(dSPECTRUM(ind_dSPECTRUM,:)),[],2); + norm_dSPECTRUM = norm_dSPECTRUM(:,ones(size(dSPECTRUM,2),1)); else - tilda_G = G; - norm_G = 1; + tilda_dSPECTRUM = dSPECTRUM; + norm_dSPECTRUM = 1; end - % store into structure - ide_spectrum.ind_G = ind_G; - ide_spectrum.tilda_G = tilda_G; - ide_spectrum.G = G; - ide_spectrum.norm_G = norm_G; + % store into structure (not everything is used later on) + ide_spectrum.ind_dSPECTRUM = ind_dSPECTRUM; + ide_spectrum.norm_dSPECTRUM = norm_dSPECTRUM; + ide_spectrum.tilda_dSPECTRUM = tilda_dSPECTRUM; + ide_spectrum.dSPECTRUM = dSPECTRUM; end - %% Perform identification checks, i.e. find out which parameters are involved +%% Perform identification checks, i.e. find out which parameters are involved if checks_via_subsets % identification_checks_via_subsets is only for debugging - [ide_lre, ide_reducedform, ide_moments, ide_spectrum, ide_minimal] = ... - identification_checks_via_subsets(ide_lre, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, totparam_nbr, modparam_nbr, options_ident); + [ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal] = ... + identification_checks_via_subsets(ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, totparam_nbr, modparam_nbr, options_ident); else - [ide_lre.cond, ide_lre.rank, ide_lre.ind0, ide_lre.indno, ide_lre.ino, ide_lre.Mco, ide_lre.Pco, ide_lre.jweak, ide_lre.jweak_pair] = ... - identification_checks(dLRE(ind_dLRE,:)./norm_dLRE, 1, tol_rank, tol_sv, modparam_nbr); + [ide_dynamic.cond, ide_dynamic.rank, ide_dynamic.ind0, ide_dynamic.indno, ide_dynamic.ino, ide_dynamic.Mco, ide_dynamic.Pco, ide_dynamic.jweak, ide_dynamic.jweak_pair] = ... + identification_checks(dDYNAMIC(ind_dDYNAMIC,:)./norm_dDYNAMIC, 1, tol_rank, tol_sv, modparam_nbr); if ~no_identification_reducedform [ide_reducedform.cond, ide_reducedform.rank, ide_reducedform.ind0, ide_reducedform.indno, ide_reducedform.ino, ide_reducedform.Mco, ide_reducedform.Pco, ide_reducedform.jweak, ide_reducedform.jweak_pair] = ... - identification_checks(dTAU(ind_dTAU,:)./norm_dTAU, 1, tol_rank, tol_sv, totparam_nbr); + identification_checks(dREDUCEDFORM(ind_dREDUCEDFORM,:)./norm_dREDUCEDFORM, 1, tol_rank, tol_sv, totparam_nbr); end if ~no_identification_moments [ide_moments.cond, ide_moments.rank, ide_moments.ind0, ide_moments.indno, ide_moments.ino, ide_moments.Mco, ide_moments.Pco, ide_moments.jweak, ide_moments.jweak_pair] = ... - identification_checks(J(ind_J,:)./norm_J, 1, tol_rank, tol_sv, totparam_nbr); + identification_checks(dMOMENTS(ind_dMOMENTS,:)./norm_dMOMENTS, 1, tol_rank, tol_sv, totparam_nbr); end if ~no_identification_minimal [ide_minimal.cond, ide_minimal.rank, ide_minimal.ind0, ide_minimal.indno, ide_minimal.ino, ide_minimal.Mco, ide_minimal.Pco, ide_minimal.jweak, ide_minimal.jweak_pair] = ... - identification_checks(D(ind_D,:)./norm_D, 2, tol_rank, tol_sv, totparam_nbr); + identification_checks(dMINIMAL(ind_dMINIMAL,:)./norm_dMINIMAL, 2, tol_rank, tol_sv, totparam_nbr); end if ~no_identification_spectrum [ide_spectrum.cond, ide_spectrum.rank, ide_spectrum.ind0, ide_spectrum.indno, ide_spectrum.ino, ide_spectrum.Mco, ide_spectrum.Pco, ide_spectrum.jweak, ide_spectrum.jweak_pair] = ... - identification_checks(tilda_G, 3, tol_rank, tol_sv, totparam_nbr); + identification_checks(tilda_dSPECTRUM, 3, tol_rank, tol_sv, totparam_nbr); end end end \ No newline at end of file diff --git a/matlab/identification_checks.m b/matlab/identification_checks.m index 4b4ba057c..4c7c23102 100644 --- a/matlab/identification_checks.m +++ b/matlab/identification_checks.m @@ -48,7 +48,9 @@ function [condX, rankX, ind0, indno, ixno, Mco, Pco, jweak, jweak_pair] = identi % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . % ========================================================================= - +if issparse(X) + X = full(X); +end if nargin < 3 || isempty(tol_rank) tol_rank = 1.e-10; %tolerance level used for rank computations end diff --git a/matlab/identification_checks_via_subsets.m b/matlab/identification_checks_via_subsets.m index 3cb9d3ece..cb5589f6b 100644 --- a/matlab/identification_checks_via_subsets.m +++ b/matlab/identification_checks_via_subsets.m @@ -1,5 +1,5 @@ -function [ide_lre, ide_reducedform, ide_moments, ide_spectrum, ide_minimal] = identification_checks_via_subsets(ide_lre, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, totparam_nbr, modparam_nbr, options_ident) -%[ide_lre, ide_reducedform, ide_moments, ide_spectrum, ide_minimal] = identification_checks_via_subsets(ide_lre, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, totparam_nbr, modparam_nbr, options_ident) +function [ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal] = identification_checks_via_subsets(ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, totparam_nbr, modparam_nbr, options_ident) +%[ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal] = identification_checks_via_subsets(ide_dynamic, ide_reducedform, ide_moments, ide_spectrum, ide_minimal, totparam_nbr, modparam_nbr, options_ident) % ------------------------------------------------------------------------- % Finds problematic sets of paramters via checking the necessary rank condition % of the Jacobians for all possible combinations of parameters. The rank is @@ -70,14 +70,14 @@ function [ide_lre, ide_reducedform, ide_moments, ide_spectrum, ide_minimal] = id % ========================================================================= %% initialize output objects and get options -no_identification_lre = 0; %always compute lre +no_identification_dynamic = 0; %always compute dynamic no_identification_reducedform = options_ident.no_identification_reducedform; no_identification_moments = options_ident.no_identification_moments; no_identification_spectrum = options_ident.no_identification_spectrum; no_identification_minimal = options_ident.no_identification_minimal; tol_rank = options_ident.tol_rank; max_dim_subsets_groups = options_ident.max_dim_subsets_groups; -lre_problpars = cell(1,max_dim_subsets_groups); +dynamic_problpars = cell(1,max_dim_subsets_groups); reducedform_problpars = cell(1,max_dim_subsets_groups); moments_problpars = cell(1,max_dim_subsets_groups); spectrum_problpars = cell(1,max_dim_subsets_groups); @@ -87,298 +87,298 @@ indtotparam = 1:totparam_nbr; %initialize index of parameters %% Prepare Jacobians and check rank % initialize linear rational expectations model -if ~no_identification_lre - dLRE = ide_lre.dLRE; - dLRE(ide_lre.ind_dLRE,:) = dLRE(ide_lre.ind_dLRE,:)./ide_lre.norm_dLRE; %normalize +if ~no_identification_dynamic + dDYNAMIC = ide_dynamic.dDYNAMIC; + dDYNAMIC(ide_dynamic.ind_dDYNAMIC,:) = dDYNAMIC(ide_dynamic.ind_dDYNAMIC,:)./ide_dynamic.norm_dDYNAMIC; %normalize if totparam_nbr > modparam_nbr - dLRE = [zeros(size(ide_lre.dLRE,1),totparam_nbr-modparam_nbr) dLRE]; %add derivatives wrt stderr and corr parameters + dDYNAMIC = [zeros(size(ide_dynamic.dDYNAMIC,1),totparam_nbr-modparam_nbr) dDYNAMIC]; %add derivatives wrt stderr and corr parameters end - rank_dLRE = rank(dLRE,tol_rank); %compute rank with imposed tolerance level - ide_lre.rank = rank_dLRE; + rank_dDYNAMIC = rank(full(dDYNAMIC),tol_rank); %compute rank with imposed tolerance level + ide_dynamic.rank = rank_dDYNAMIC; % check rank criteria for full Jacobian - if rank_dLRE == totparam_nbr + if rank_dDYNAMIC == totparam_nbr % all parameters are identifiable - no_identification_lre = 1; %skip in the following - indparam_dLRE = []; + no_identification_dynamic = 1; %skip in the following + indparam_dDYNAMIC = []; else % there is lack of identification - indparam_dLRE = indtotparam; %initialize for nchoosek + indparam_dDYNAMIC = indtotparam; %initialize for nchoosek end else - indparam_dLRE = []; %empty for nchoosek + indparam_dDYNAMIC = []; %empty for nchoosek end % initialize for reduced form solution criteria if ~no_identification_reducedform - dTAU = ide_reducedform.dTAU; - dTAU(ide_reducedform.ind_dTAU,:) = dTAU(ide_reducedform.ind_dTAU,:)./ide_reducedform.norm_dTAU; %normalize - rank_dTAU = rank(dTAU,tol_rank); %compute rank with imposed tolerance level - ide_reducedform.rank = rank_dTAU; + dREDUCEDFORM = ide_reducedform.dREDUCEDFORM; + dREDUCEDFORM(ide_reducedform.ind_dREDUCEDFORM,:) = dREDUCEDFORM(ide_reducedform.ind_dREDUCEDFORM,:)./ide_reducedform.norm_dREDUCEDFORM; %normalize + rank_dREDUCEDFORM = rank(full(dREDUCEDFORM),tol_rank); %compute rank with imposed tolerance level + ide_reducedform.rank = rank_dREDUCEDFORM; % check rank criteria for full Jacobian - if rank_dTAU == totparam_nbr + if rank_dREDUCEDFORM == totparam_nbr % all parameters are identifiable no_identification_reducedform = 1; %skip in the following - indparam_dTAU = []; + indparam_dREDUCEDFORM = []; else % there is lack of identification - indparam_dTAU = indtotparam; %initialize for nchoosek + indparam_dREDUCEDFORM = indtotparam; %initialize for nchoosek end else - indparam_dTAU = []; %empty for nchoosek + indparam_dREDUCEDFORM = []; %empty for nchoosek end % initialize for moments criteria if ~no_identification_moments - J = ide_moments.J; - J(ide_moments.ind_J,:) = J(ide_moments.ind_J,:)./ide_moments.norm_J; %normalize - rank_J = rank(J,tol_rank); %compute rank with imposed tolerance level - ide_moments.rank = rank_J; + dMOMENTS = ide_moments.dMOMENTS; + dMOMENTS(ide_moments.ind_dMOMENTS,:) = dMOMENTS(ide_moments.ind_dMOMENTS,:)./ide_moments.norm_dMOMENTS; %normalize + rank_dMOMENTS = rank(full(dMOMENTS),tol_rank); %compute rank with imposed tolerance level + ide_moments.rank = rank_dMOMENTS; % check rank criteria for full Jacobian - if rank_J == totparam_nbr + if rank_dMOMENTS == totparam_nbr % all parameters are identifiable no_identification_moments = 1; %skip in the following - indparam_J = []; + indparam_dMOMENTS = []; else % there is lack of identification - indparam_J = indtotparam; %initialize for nchoosek + indparam_dMOMENTS = indtotparam; %initialize for nchoosek end else - indparam_J = []; %empty for nchoosek + indparam_dMOMENTS = []; %empty for nchoosek end % initialize for spectrum criteria if ~no_identification_spectrum - G = ide_spectrum.tilda_G; %tilda G is normalized G matrix in identification_analysis.m + dSPECTRUM = ide_spectrum.tilda_dSPECTRUM; %tilda dSPECTRUM is normalized dSPECTRUM matrix in identification_analysis.m %alternative normalization - %G = ide_spectrum.G; - %G(ide_spectrum.ind_G,:) = G(ide_spectrum.ind_G,:)./ide_spectrum.norm_G; %normalize - rank_G = rank(G,tol_rank); %compute rank with imposed tolerance level - ide_spectrum.rank = rank_G; + %dSPECTRUM = ide_spectrum.dSPECTRUM; + %dSPECTRUM(ide_spectrum.ind_dSPECTRUM,:) = dSPECTRUM(ide_spectrum.ind_dSPECTRUM,:)./ide_spectrum.norm_dSPECTRUM; %normalize + rank_dSPECTRUM = rank(full(dSPECTRUM),tol_rank); %compute rank with imposed tolerance level + ide_spectrum.rank = rank_dSPECTRUM; % check rank criteria for full Jacobian - if rank_G == totparam_nbr + if rank_dSPECTRUM == totparam_nbr % all parameters are identifiable no_identification_spectrum = 1; %skip in the following - indparam_G = []; + indparam_dSPECTRUM = []; else % lack of identification - indparam_G = indtotparam; %initialize for nchoosek + indparam_dSPECTRUM = indtotparam; %initialize for nchoosek end else - indparam_G = []; %empty for nchoosek + indparam_dSPECTRUM = []; %empty for nchoosek end % initialize for minimal system criteria if ~no_identification_minimal - D = ide_minimal.D; - D(ide_minimal.ind_D,:) = D(ide_minimal.ind_D,:)./ide_minimal.norm_D; %normalize - D_par = D(:,1:totparam_nbr); %part of D that is dependent on parameters - D_rest = D(:,(totparam_nbr+1):end); %part of D that is independent of parameters - rank_D = rank(D,tol_rank); %compute rank via SVD see function below - ide_minimal.rank = rank_D; - D_fixed_rank_nbr = size(D_rest,2); + dMINIMAL = ide_minimal.dMINIMAL; + dMINIMAL(ide_minimal.ind_dMINIMAL,:) = dMINIMAL(ide_minimal.ind_dMINIMAL,:)./ide_minimal.norm_dMINIMAL; %normalize + dMINIMAL_par = dMINIMAL(:,1:totparam_nbr); %part of dMINIMAL that is dependent on parameters + dMINIMAL_rest = dMINIMAL(:,(totparam_nbr+1):end); %part of dMINIMAL that is independent of parameters + rank_dMINIMAL = rank(full(dMINIMAL),tol_rank); %compute rank via SVD see function below + ide_minimal.rank = rank_dMINIMAL; + dMINIMAL_fixed_rank_nbr = size(dMINIMAL_rest,2); % check rank criteria for full Jacobian - if rank_D == totparam_nbr + D_fixed_rank_nbr + if rank_dMINIMAL == totparam_nbr + dMINIMAL_fixed_rank_nbr % all parameters are identifiable no_identification_minimal = 1; %skip in the following - indparam_D = []; + indparam_dMINIMAL = []; else % lack of identification - indparam_D = indtotparam; %initialize for nchoosek + indparam_dMINIMAL = indtotparam; %initialize for nchoosek end else - indparam_D = []; %empty for nchoosek + indparam_dMINIMAL = []; %empty for nchoosek end %% Check single parameters for j=1:totparam_nbr - if ~no_identification_lre + if ~no_identification_dynamic % Columns correspond to single parameters, i.e. full rank would be equal to 1 - if rank(dLRE(:,j),tol_rank) == 0 - lre_problpars{1} = [lre_problpars{1};j]; + if rank(full(dDYNAMIC(:,j)),tol_rank) == 0 + dynamic_problpars{1} = [dynamic_problpars{1};j]; end end if ~no_identification_reducedform % Columns correspond to single parameters, i.e. full rank would be equal to 1 - if rank(dTAU(:,j),tol_rank) == 0 + if rank(full(dREDUCEDFORM(:,j)),tol_rank) == 0 reducedform_problpars{1} = [reducedform_problpars{1};j]; end end if ~no_identification_moments % Columns correspond to single parameters, i.e. full rank would be equal to 1 - if rank(J(:,j),tol_rank) == 0 + if rank(full(dMOMENTS(:,j)),tol_rank) == 0 moments_problpars{1} = [moments_problpars{1};j]; end end if ~no_identification_spectrum % Diagonal values correspond to single parameters, absolute value needs to be greater than tolerance level - if abs(G(j,j)) < tol_rank + if abs(dSPECTRUM(j,j)) < tol_rank spectrum_problpars{1} = [spectrum_problpars{1};j]; end end if ~no_identification_minimal - % Columns of D_par correspond to single parameters, needs to be augmented with D_rest (part that is independent of parameters), - % full rank would be equal to 1+D_fixed_rank_nbr - if rank([D_par(:,j) D_rest],tol_rank) == D_fixed_rank_nbr + % Columns of dMINIMAL_par correspond to single parameters, needs to be augmented with dMINIMAL_rest (part that is independent of parameters), + % full rank would be equal to 1+dMINIMAL_fixed_rank_nbr + if rank(full([dMINIMAL_par(:,j) dMINIMAL_rest]),tol_rank) == dMINIMAL_fixed_rank_nbr minimal_problpars{1} = [minimal_problpars{1};j]; end end end % Check whether lack of identification is only due to single parameters -if ~no_identification_lre - if size(lre_problpars{1},1) == (totparam_nbr - rank_dLRE) +if ~no_identification_dynamic + if size(dynamic_problpars{1},1) == (totparam_nbr - rank_dDYNAMIC) %found all nonidentified parameter sets - no_identification_lre = 1; %skip in the following + no_identification_dynamic = 1; %skip in the following else %still parameter sets that are nonidentified - indparam_dLRE(lre_problpars{1}) = []; %remove single unidentified parameters from higher-order sets of indparam + indparam_dDYNAMIC(dynamic_problpars{1}) = []; %remove single unidentified parameters from higher-order sets of indparam end end if ~no_identification_reducedform - if size(reducedform_problpars{1},1) == (totparam_nbr - rank_dTAU) + if size(reducedform_problpars{1},1) == (totparam_nbr - rank_dREDUCEDFORM) %found all nonidentified parameter sets no_identification_reducedform = 1; %skip in the following else %still parameter sets that are nonidentified - indparam_dTAU(reducedform_problpars{1}) = []; %remove single unidentified parameters from higher-order sets of indparam + indparam_dREDUCEDFORM(reducedform_problpars{1}) = []; %remove single unidentified parameters from higher-order sets of indparam end end if ~no_identification_moments - if size(moments_problpars{1},1) == (totparam_nbr - rank_J) + if size(moments_problpars{1},1) == (totparam_nbr - rank_dMOMENTS) %found all nonidentified parameter sets no_identification_moments = 1; %skip in the following else %still parameter sets that are nonidentified - indparam_J(moments_problpars{1}) = []; %remove single unidentified parameters from higher-order sets of indparam + indparam_dMOMENTS(moments_problpars{1}) = []; %remove single unidentified parameters from higher-order sets of indparam end end if ~no_identification_spectrum - if size(spectrum_problpars{1},1) == (totparam_nbr - rank_G) + if size(spectrum_problpars{1},1) == (totparam_nbr - rank_dSPECTRUM) %found all nonidentified parameter sets no_identification_spectrum = 1; %skip in the following else %still parameter sets that are nonidentified - indparam_G(spectrum_problpars{1}) = []; %remove single unidentified parameters from higher-order sets of indparam + indparam_dSPECTRUM(spectrum_problpars{1}) = []; %remove single unidentified parameters from higher-order sets of indparam end end if ~no_identification_minimal - if size(minimal_problpars{1},1) == (totparam_nbr + D_fixed_rank_nbr - rank_D) + if size(minimal_problpars{1},1) == (totparam_nbr + dMINIMAL_fixed_rank_nbr - rank_dMINIMAL) %found all nonidentified parameter sets no_identification_minimal = 1; %skip in the following else %still parameter sets that are nonidentified - indparam_D(minimal_problpars{1}) = []; %remove single unidentified parameters from higher-order sets of indparam + indparam_dMINIMAL(minimal_problpars{1}) = []; %remove single unidentified parameters from higher-order sets of indparam end end %% check higher order (j>1) parameter sets -%get common parameter indices from which to sample higher-order sets using nchoosek (we do not want to run nchoosek three times), most of the times indparamJ, indparamG, and indparamD are equal anyways -indtotparam = unique([indparam_dLRE indparam_dTAU indparam_J indparam_G indparam_D]); +%get common parameter indices from which to sample higher-order sets using nchoosek (we do not want to run nchoosek three times), most of the times indparamdMOMENTS, indparamdSPECTRUM, and indparamdMINIMAL are equal anyways +indtotparam = unique([indparam_dDYNAMIC indparam_dREDUCEDFORM indparam_dMOMENTS indparam_dSPECTRUM indparam_dMINIMAL]); for j=2:min(length(indtotparam),max_dim_subsets_groups) % Check j-element subsets h = dyn_waitbar(0,['Brute force collinearity for ' int2str(j) ' parameters.']); %Step1: get all possible unique subsets of j elements - if ~no_identification_lre || ~no_identification_reducedform || ~no_identification_moments || ~no_identification_spectrum || ~no_identification_minimal + if ~no_identification_dynamic || ~no_identification_reducedform || ~no_identification_moments || ~no_identification_spectrum || ~no_identification_minimal indexj=nchoosek(int16(indtotparam),j); % int16 speeds up nchoosek % One could also use a mex version of nchoosek to speed this up, e.g.VChooseK from https://de.mathworks.com/matlabcentral/fileexchange/26190-vchoosek end %Step 2: remove already problematic sets and initialize rank vector - if ~no_identification_lre - indexj_dLRE = RemoveProblematicParameterSets(indexj,lre_problpars); - rankj_dLRE = zeros(size(indexj_dLRE,1),1); + if ~no_identification_dynamic + indexj_dDYNAMIC = RemoveProblematicParameterSets(indexj,dynamic_problpars); + rankj_dDYNAMIC = zeros(size(indexj_dDYNAMIC,1),1); else - indexj_dLRE = []; + indexj_dDYNAMIC = []; end if ~no_identification_reducedform - indexj_dTAU = RemoveProblematicParameterSets(indexj,reducedform_problpars); - rankj_dTAU = zeros(size(indexj_dTAU,1),1); + indexj_dREDUCEDFORM = RemoveProblematicParameterSets(indexj,reducedform_problpars); + rankj_dREDUCEDFORM = zeros(size(indexj_dREDUCEDFORM,1),1); else - indexj_dTAU = []; + indexj_dREDUCEDFORM = []; end if ~no_identification_moments - indexj_J = RemoveProblematicParameterSets(indexj,moments_problpars); - rankj_J = zeros(size(indexj_J,1),1); + indexj_dMOMENTS = RemoveProblematicParameterSets(indexj,moments_problpars); + rankj_dMOMENTS = zeros(size(indexj_dMOMENTS,1),1); else - indexj_J = []; + indexj_dMOMENTS = []; end if ~no_identification_spectrum - indexj_G = RemoveProblematicParameterSets(indexj,spectrum_problpars); - rankj_G = zeros(size(indexj_G,1),1); + indexj_dSPECTRUM = RemoveProblematicParameterSets(indexj,spectrum_problpars); + rankj_dSPECTRUM = zeros(size(indexj_dSPECTRUM,1),1); else - indexj_G = []; + indexj_dSPECTRUM = []; end if ~no_identification_minimal - indexj_D = RemoveProblematicParameterSets(indexj,minimal_problpars); - rankj_D = zeros(size(indexj_D,1),1); + indexj_dMINIMAL = RemoveProblematicParameterSets(indexj,minimal_problpars); + rankj_dMINIMAL = zeros(size(indexj_dMINIMAL,1),1); else - indexj_D = []; + indexj_dMINIMAL = []; end %Step3: Check rank criteria on submatrices - k_dLRE=0; k_dTAU=0; k_J=0; k_G=0; k_D=0; %initialize counters - maxk = max([size(indexj_dLRE,1), size(indexj_dTAU,1), size(indexj_J,1), size(indexj_D,1), size(indexj_G,1)]); + k_dDYNAMIC=0; k_dREDUCEDFORM=0; k_dMOMENTS=0; k_dSPECTRUM=0; k_dMINIMAL=0; %initialize counters + maxk = max([size(indexj_dDYNAMIC,1), size(indexj_dREDUCEDFORM,1), size(indexj_dMOMENTS,1), size(indexj_dMINIMAL,1), size(indexj_dSPECTRUM,1)]); for k=1:maxk - if ~no_identification_lre - k_dLRE = k_dLRE+1; - if k_dLRE <= size(indexj_dLRE,1) - dLRE_j = dLRE(:,indexj_dLRE(k_dLRE,:)); % pick columns that correspond to parameter subset - rankj_dLRE(k_dLRE,1) = rank(dLRE_j,tol_rank); %compute rank with imposed tolerance + if ~no_identification_dynamic + k_dDYNAMIC = k_dDYNAMIC+1; + if k_dDYNAMIC <= size(indexj_dDYNAMIC,1) + dDYNAMIC_j = dDYNAMIC(:,indexj_dDYNAMIC(k_dDYNAMIC,:)); % pick columns that correspond to parameter subset + rankj_dDYNAMIC(k_dDYNAMIC,1) = rank(full(dDYNAMIC_j),tol_rank); %compute rank with imposed tolerance end end if ~no_identification_reducedform - k_dTAU = k_dTAU+1; - if k_dTAU <= size(indexj_dTAU,1) - dTAU_j = dTAU(:,indexj_dTAU(k_dTAU,:)); % pick columns that correspond to parameter subset - rankj_dTAU(k_dTAU,1) = rank(dTAU_j,tol_rank); %compute rank with imposed tolerance + k_dREDUCEDFORM = k_dREDUCEDFORM+1; + if k_dREDUCEDFORM <= size(indexj_dREDUCEDFORM,1) + dREDUCEDFORM_j = dREDUCEDFORM(:,indexj_dREDUCEDFORM(k_dREDUCEDFORM,:)); % pick columns that correspond to parameter subset + rankj_dREDUCEDFORM(k_dREDUCEDFORM,1) = rank(full(dREDUCEDFORM_j),tol_rank); %compute rank with imposed tolerance end end if ~no_identification_moments - k_J = k_J+1; - if k_J <= size(indexj_J,1) - J_j = J(:,indexj_J(k_J,:)); % pick columns that correspond to parameter subset - rankj_J(k_J,1) = rank(J_j,tol_rank); %compute rank with imposed tolerance + k_dMOMENTS = k_dMOMENTS+1; + if k_dMOMENTS <= size(indexj_dMOMENTS,1) + dMOMENTS_j = dMOMENTS(:,indexj_dMOMENTS(k_dMOMENTS,:)); % pick columns that correspond to parameter subset + rankj_dMOMENTS(k_dMOMENTS,1) = rank(full(dMOMENTS_j),tol_rank); %compute rank with imposed tolerance end end if ~no_identification_minimal - k_D = k_D+1; - if k_D <= size(indexj_D,1) - D_j = [D_par(:,indexj_D(k_D,:)) D_rest]; % pick columns in D_par that correspond to parameter subset and augment with parameter-indepdendet part D_rest - rankj_D(k_D,1) = rank(D_j,tol_rank); %compute rank with imposed tolerance + k_dMINIMAL = k_dMINIMAL+1; + if k_dMINIMAL <= size(indexj_dMINIMAL,1) + dMINIMAL_j = [dMINIMAL_par(:,indexj_dMINIMAL(k_dMINIMAL,:)) dMINIMAL_rest]; % pick columns in dMINIMAL_par that correspond to parameter subset and augment with parameter-indepdendet part dMINIMAL_rest + rankj_dMINIMAL(k_dMINIMAL,1) = rank(full(dMINIMAL_j),tol_rank); %compute rank with imposed tolerance end end if ~no_identification_spectrum - k_G = k_G+1; - if k_G <= size(indexj_G,1) - G_j = G(indexj_G(k_G,:),indexj_G(k_G,:)); % pick rows and columns that correspond to parameter subset - rankj_G(k_G,1) = rank(G_j,tol_rank); % Compute rank with imposed tol + k_dSPECTRUM = k_dSPECTRUM+1; + if k_dSPECTRUM <= size(indexj_dSPECTRUM,1) + dSPECTRUM_j = dSPECTRUM(indexj_dSPECTRUM(k_dSPECTRUM,:),indexj_dSPECTRUM(k_dSPECTRUM,:)); % pick rows and columns that correspond to parameter subset + rankj_dSPECTRUM(k_dSPECTRUM,1) = rank(full(dSPECTRUM_j),tol_rank); % Compute rank with imposed tol end end dyn_waitbar(k/maxk,h) end %Step 4: Compare rank conditions for all possible subsets. If rank condition is violated, then the corresponding numbers of the parameters are stored - if ~no_identification_lre - lre_problpars{j} = indexj_dLRE(rankj_dLRE < j,:); + if ~no_identification_dynamic + dynamic_problpars{j} = indexj_dDYNAMIC(rankj_dDYNAMIC < j,:); end if ~no_identification_reducedform - reducedform_problpars{j} = indexj_dTAU(rankj_dTAU < j,:); + reducedform_problpars{j} = indexj_dREDUCEDFORM(rankj_dREDUCEDFORM < j,:); end if ~no_identification_moments - moments_problpars{j} = indexj_J(rankj_J < j,:); + moments_problpars{j} = indexj_dMOMENTS(rankj_dMOMENTS < j,:); end if ~no_identification_minimal - minimal_problpars{j} = indexj_D(rankj_D < (j+D_fixed_rank_nbr),:); + minimal_problpars{j} = indexj_dMINIMAL(rankj_dMINIMAL < (j+dMINIMAL_fixed_rank_nbr),:); end if ~no_identification_spectrum - spectrum_problpars{j} = indexj_G(rankj_G < j,:); + spectrum_problpars{j} = indexj_dSPECTRUM(rankj_dSPECTRUM < j,:); end % % Optional Step 5: % remove redundant 2-sets, eg. if the problematic sets are [(p1,p2);(p1,p3);(p2,p3)], then the unique problematic parameter sets are actually only [(p1,p2),(p1,p3)] % if j == 2 -% for jj=1:max([size(lre_problpars{2},1), size(reducedform_problpars{2},1), size(moments_problpars{2},1), size(spectrum_problpars{2},1), size(minimal_problpars{2},1)]) -% if jj <= size(lre_problpars{2},1) -% lre_problpars{2}(lre_problpars{2}(jj,2)==lre_problpars{2}(:,1)) = lre_problpars{2}(jj,1); +% for jj=1:max([size(dynamic_problpars{2},1), size(reducedform_problpars{2},1), size(moments_problpars{2},1), size(spectrum_problpars{2},1), size(minimal_problpars{2},1)]) +% if jj <= size(dynamic_problpars{2},1) +% dynamic_problpars{2}(dynamic_problpars{2}(jj,2)==dynamic_problpars{2}(:,1)) = dynamic_problpars{2}(jj,1); % end % if jj <= size(reducedform_problpars{2},1) % reducedform_problpars{2}(reducedform_problpars{2}(jj,2)==reducedform_problpars{2}(:,1)) = reducedform_problpars{2}(jj,1); @@ -393,13 +393,13 @@ for j=2:min(length(indtotparam),max_dim_subsets_groups) % Check j-element subset % minimal_problpars{2}(minimal_problpars{2}(jj,2)==minimal_problpars{2}(:,1)) = minimal_problpars{2}(jj,1); % end % end -% lre_problpars{2} = unique(lre_problpars{2},'rows'); +% dynamic_problpars{2} = unique(dynamic_problpars{2},'rows'); % reducedform_problpars{2} = unique(reducedform_problpars{2},'rows'); % moments_problpars{2} = unique(moments_problpars{2},'rows'); % spectrum_problpars{2} = unique(spectrum_problpars{2},'rows'); % minimal_problpars{2} = unique(minimal_problpars{2},'rows'); % % in indparam we replace the second parameter of problematic 2-sets by the observational equivalent first parameter to speed up nchoosek -% idx2 = unique([lre_problpars{2}; reducedform_problpars{2}; moments_problpars{2}; spectrum_problpars{2}; minimal_problpars{2}],'rows'); +% idx2 = unique([dynamic_problpars{2}; reducedform_problpars{2}; moments_problpars{2}; spectrum_problpars{2}; minimal_problpars{2}],'rows'); % if ~isempty(idx2) % indtotparam(ismember(indtotparam,idx2(:,2))) = []; % end @@ -408,10 +408,10 @@ for j=2:min(length(indtotparam),max_dim_subsets_groups) % Check j-element subset end %% Save output variables -if ~isempty(lre_problpars{1}) - lre_problpars{1}(ismember(lre_problpars{1},1:(totparam_nbr-modparam_nbr))) = []; % get rid of stderr and corr variables for lre +if ~isempty(dynamic_problpars{1}) + dynamic_problpars{1}(ismember(dynamic_problpars{1},1:(totparam_nbr-modparam_nbr))) = []; % get rid of stderr and corr variables for dynamic end -ide_lre.problpars = lre_problpars; +ide_dynamic.problpars = dynamic_problpars; ide_reducedform.problpars = reducedform_problpars; ide_moments.problpars = moments_problpars; ide_spectrum.problpars = spectrum_problpars; diff --git a/matlab/plot_identification.m b/matlab/plot_identification.m index 48fff502e..e534b4317 100644 --- a/matlab/plot_identification.m +++ b/matlab/plot_identification.m @@ -47,43 +47,43 @@ if nargin <11 end [SampleSize, nparam]=size(params); -si_Jnorm = idemoments.si_Jnorm; -si_dTAUnorm = idemodel.si_dTAUnorm; -si_dLREnorm = idelre.si_dLREnorm; +si_dMOMENTSnorm = idemoments.si_dMOMENTSnorm; +si_dTAUnorm = idemodel.si_dREDUCEDFORMnorm; +si_dLREnorm = idelre.si_dDYNAMICnorm; tittxt1=regexprep(tittxt, ' ', '_'); tittxt1=strrep(tittxt1, '.', ''); if SampleSize == 1 - si_J = idemoments.si_J; + si_dMOMENTS = idemoments.si_dMOMENTS; hh = dyn_figure(options_.nodisplay,'Name',[tittxt, ' - Identification using info from observables']); subplot(211) - mmm = (idehess.ide_strength_J); + mmm = (idehess.ide_strength_dMOMENTS); [ss, is] = sort(mmm); - if ~all(isnan(idehess.ide_strength_J_prior)) - bar(log([idehess.ide_strength_J(:,is)' idehess.ide_strength_J_prior(:,is)'])) + if ~all(isnan(idehess.ide_strength_dMOMENTS_prior)) + bar(log([idehess.ide_strength_dMOMENTS(:,is)' idehess.ide_strength_dMOMENTS_prior(:,is)'])) else - bar(log([idehess.ide_strength_J(:,is)' ])) + bar(log([idehess.ide_strength_dMOMENTS(:,is)' ])) end hold on - plot((1:length(idehess.ide_strength_J(:,is)))-0.15,log([idehess.ide_strength_J(:,is)']),'o','MarkerSize',7,'MarkerFaceColor',[0 0 0],'MarkerEdgeColor','none') - plot((1:length(idehess.ide_strength_J_prior(:,is)))+0.15,log([idehess.ide_strength_J_prior(:,is)']),'o','MarkerSize',7,'MarkerFaceColor',[0 0 0],'MarkerEdgeColor','none') - if any(isinf(log(idehess.ide_strength_J(idehess.identified_parameter_indices)))) + plot((1:length(idehess.ide_strength_dMOMENTS(:,is)))-0.15,log([idehess.ide_strength_dMOMENTS(:,is)']),'o','MarkerSize',7,'MarkerFaceColor',[0 0 0],'MarkerEdgeColor','none') + plot((1:length(idehess.ide_strength_dMOMENTS_prior(:,is)))+0.15,log([idehess.ide_strength_dMOMENTS_prior(:,is)']),'o','MarkerSize',7,'MarkerFaceColor',[0 0 0],'MarkerEdgeColor','none') + if any(isinf(log(idehess.ide_strength_dMOMENTS(idehess.identified_parameter_indices)))) %-Inf, i.e. 0 strength - inf_indices=find(isinf(log(idehess.ide_strength_J(idehess.identified_parameter_indices))) & log(idehess.ide_strength_J(idehess.identified_parameter_indices))<0); + inf_indices=find(isinf(log(idehess.ide_strength_dMOMENTS(idehess.identified_parameter_indices))) & log(idehess.ide_strength_dMOMENTS(idehess.identified_parameter_indices))<0); inf_pos=ismember(is,idehess.identified_parameter_indices(inf_indices)); plot(find(inf_pos)-0.15,zeros(sum(inf_pos),1),'o','MarkerSize',7,'MarkerFaceColor',[1 0 0],'MarkerEdgeColor',[0 0 0]) %+Inf, i.e. Inf strength - inf_indices=find(isinf(log(idehess.ide_strength_J(idehess.identified_parameter_indices))) & log(idehess.ide_strength_J(idehess.identified_parameter_indices))>0); + inf_indices=find(isinf(log(idehess.ide_strength_dMOMENTS(idehess.identified_parameter_indices))) & log(idehess.ide_strength_dMOMENTS(idehess.identified_parameter_indices))>0); inf_pos=ismember(is,idehess.identified_parameter_indices(inf_indices)); plot(find(inf_pos)-0.15,zeros(sum(inf_pos),1),'o','MarkerSize',7,'MarkerFaceColor',[1 1 1],'MarkerEdgeColor',[0 0 0]) end - if any(isinf(log(idehess.ide_strength_J_prior(idehess.identified_parameter_indices)))) + if any(isinf(log(idehess.ide_strength_dMOMENTS_prior(idehess.identified_parameter_indices)))) %-Inf, i.e. 0 strength - inf_indices=find(isinf(log(idehess.ide_strength_J_prior(idehess.identified_parameter_indices))) & log(idehess.ide_strength_J_prior(idehess.identified_parameter_indices))<0); + inf_indices=find(isinf(log(idehess.ide_strength_dMOMENTS_prior(idehess.identified_parameter_indices))) & log(idehess.ide_strength_dMOMENTS_prior(idehess.identified_parameter_indices))<0); inf_pos=ismember(is,idehess.identified_parameter_indices(inf_indices)); plot(find(inf_pos)+0.15,zeros(sum(inf_pos),1),'o','MarkerSize',7,'MarkerFaceColor',[1 0 0],'MarkerEdgeColor',[0 0 0]) %+Inf, i.e. 0 strength - inf_indices=find(isinf(log(idehess.ide_strength_J_prior(idehess.identified_parameter_indices))) & log(idehess.ide_strength_J_prior(idehess.identified_parameter_indices))>0); + inf_indices=find(isinf(log(idehess.ide_strength_dMOMENTS_prior(idehess.identified_parameter_indices))) & log(idehess.ide_strength_dMOMENTS_prior(idehess.identified_parameter_indices))>0); inf_pos=ismember(is,idehess.identified_parameter_indices(inf_indices)); plot(find(inf_pos)+0.15,zeros(sum(inf_pos),1),'o','MarkerSize',7,'MarkerFaceColor',[1 1 1],'MarkerEdgeColor',[0 0 0]) end @@ -93,7 +93,7 @@ if SampleSize == 1 for ip=1:nparam text(ip,dy(1),name{is(ip)},'rotation',90,'HorizontalAlignment','right','interpreter','none') end - if ~all(isnan(idehess.ide_strength_J_prior)) + if ~all(isnan(idehess.ide_strength_dMOMENTS_prior)) legend('relative to param value','relative to prior std','Location','Best') else legend('relative to param value','Location','Best') @@ -161,12 +161,12 @@ if SampleSize == 1 skipline() disp('Plotting advanced diagnostics') end - if all(isnan([si_Jnorm';si_dTAUnorm';si_dLREnorm'])) + if all(isnan([si_dMOMENTSnorm';si_dTAUnorm';si_dLREnorm'])) fprintf('\nIDENTIFICATION: Skipping sensitivity plot, because standard deviation of parameters is NaN, possibly due to the use of ML.\n') else hh = dyn_figure(options_.nodisplay,'Name',[tittxt, ' - Sensitivity plot']); subplot(211) - mmm = (si_Jnorm)'./max(si_Jnorm); + mmm = (si_dMOMENTSnorm)'./max(si_dMOMENTSnorm); mmm1 = (si_dTAUnorm)'./max(si_dTAUnorm); mmm=[mmm mmm1]; mmm1 = (si_dLREnorm)'./max(si_dLREnorm); @@ -199,7 +199,7 @@ if SampleSize == 1 end end % identificaton patterns - for j=1:size(idemoments.cosnJ,2) + for j=1:size(idemoments.cosndMOMENTS,2) pax=NaN(nparam,nparam); % fprintf('\n') % disp(['Collinearity patterns with ', int2str(j) ,' parameter(s)']) @@ -212,10 +212,10 @@ if SampleSize == 1 namx=[namx ' ' sprintf('%-15s','--')]; else namx=[namx ' ' sprintf('%-15s',name{dumpindx})]; - pax(i,dumpindx)=idemoments.cosnJ(i,j); + pax(i,dumpindx)=idemoments.cosndMOMENTS(i,j); end end - % fprintf('%-15s [%s] %10.3f\n',name{i},namx,idemoments.cosnJ(i,j)) + % fprintf('%-15s [%s] %10.3f\n',name{i},namx,idemoments.cosndMOMENTS(i,j)) end hh = dyn_figure(options_.nodisplay,'Name',[tittxt,' - Collinearity patterns with ', int2str(j) ,' parameter(s)']); imagesc(pax,[0 1]); @@ -337,9 +337,9 @@ if SampleSize == 1 else hh = dyn_figure(options_.nodisplay,'Name',['MC sensitivities']); subplot(211) - mmm = (idehess.ide_strength_J); + mmm = (idehess.ide_strength_dMOMENTS); [ss, is] = sort(mmm); - mmm = mean(si_Jnorm)'; + mmm = mean(si_dMOMENTSnorm)'; mmm = mmm./max(mmm); if advanced mmm1 = mean(si_dTAUnorm)'; diff --git a/matlab/prodmom.m b/matlab/prodmom.m new file mode 100644 index 000000000..230085083 --- /dev/null +++ b/matlab/prodmom.m @@ -0,0 +1,80 @@ +% +% prodmom.m Date: 4/29/2006 +% This Matlab program 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). +% V only needs to be positive semidefinite. +% V: variance-covariance matrix of X +% ii: vector of i_j +% nu: power of X_{i_j} +% Reference: 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 +% Note that there is a typo in Eq.(46.25), there should be an extra rho in front +% of the equation. +% Usage: prodmom(V,[i1 i2 ... ir],[nu1 nu2 ... nur]) +% Example: To get E[X_2X_4^3X_7^2], use prodmom(V,[2 4 7],[1 3 2]) +% +function y = prodmom(V,ii,nu); +if nargin<3 + nu = ones(size(ii)); +end +s = sum(nu); +if s==0 + y = 1; + return +end +if rem(s,2)==1 + y = 0; + return +end +nuz = nu==0; +nu(nuz) = []; +ii(nuz) = []; +m = length(ii); +V = V(ii,ii); +s2 = s/2; +% +% Use univariate normal results +% +if m==1 + y = V^s2*prod([1:2:s-1]); + 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)); + y = V(1,1)^(nu(1)/2)*V(2,2)^(nu(2)/2)*bivmom(nu,rho); + return +end +% +% Regular case +% +[nu,inu] = sort(nu,2,'descend'); +V = V(inu,inu); % Extract only the relevant part of V +x = zeros(1,m); +V = V./2; +nu2 = nu./2; +p = 2; +q = nu2*V*nu2'; +y = 0; +for i=1:fix(prod(nu+1)/2) + y = y+p*q^s2; + for j=1:m + if x(j) 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 +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); +end + +if order > 1 + % 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); + end + for j = 1:size(COMBOS4,1) + E_u_u_u_u(j) = prodmom(E_uu, 1:u_nbr, COMBOS4(j,:)); + 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))); + 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[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,:); + + 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(:); + + 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)); + + 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)); + + 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); + + 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); + + 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)); + + dd(1:y_nbr , jp) = 0.5*dghs2(indy,jp) + 0.5*(dghuu(indy,:,jp)*E_uu(:) + ghuu(indy,:)*vec(dE_uu(:,:,jp))); + + 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))); + 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); + + 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))]; + 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 + + 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,:)); + + 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)); + end + end + + 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_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_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); + + 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); + + 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(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(:); + + 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; + + + 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)); + + 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); + + 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 + + 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))); + + 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)); + + 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 + + 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 + + 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 + + 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 + + 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 + + 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)'; + 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); + end + end + end +end + + +E_y = Yss(indy,:) + C*E_z + d; +Om_z = B*Varinov*transpose(B); +Om_y = D*Varinov*transpose(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 + end + end +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; +end diff --git a/matlab/quadruplication.m b/matlab/quadruplication.m new file mode 100644 index 000000000..462857ff0 --- /dev/null +++ b/matlab/quadruplication.m @@ -0,0 +1,83 @@ +% 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 +% +% Inputs: +% p: size of vector +% Outputs: +% QP: quadruplication matrix +% QPinv: Moore-Penrose inverse of QP +% +function [QP,QPinv] = quadruplication(p,progress,sparseflag) + +if nargin <2 + progress =0; +end +if nargin < 3 + sparseflag = 1; +end +reverseStr = ''; counti=1; +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 +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 +end + +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); + n = ii + (jj-1)*p + (kk-1)*p^2 + (ll-1)*p^3; + m = mue(p,i,j,k,l); + QP(n,m)=1; + if nargout > 1 + if i==j && j==k && k==l + QPinv(m,n)=1; + elseif i==j && j==k && k>l + QPinv(m,n)=1/4; + elseif i>j && j==k && k==l + QPinv(m,n)=1/4; + elseif i==j && j>k && k==l + QPinv(m,n) = 1/6; + elseif i>j && j>k && k==l + QPinv(m,n) = 1/12; + elseif i>j && j==k && k>l + QPinv(m,n) = 1/12; + elseif i==j && j>k && k>l + QPinv(m,n) = 1/12; + elseif i>j && j>k && k>l + QPinv(m,n) = 1/24; + end + end + end + counti = counti+1; + end + end + end +end +%QPinv = (transpose(QP)*QP)\transpose(QP); + +function m = mue(p,i,j,k,l) + 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 + + +end \ No newline at end of file diff --git a/matlab/uperm.m b/matlab/uperm.m new file mode 100644 index 000000000..af0d1f6f1 --- /dev/null +++ b/matlab/uperm.m @@ -0,0 +1,27 @@ +% By Willi Mutschler, September 26, 2016. Email: willi@mutschler.eu +function p = uperm(a) +[u, ~, J] = unique(a); +p = u(up(J, length(a))); + + +function p = up(J, n) +ktab = histcounts(J,1:max(J)); +l = n; +p = zeros(1, n); +s = 1; +for i=1:length(ktab) + k = ktab(i); + c = nchoosek(1:l, k); + m = size(c,1); + [t, ~] = find(~p.'); + t = reshape(t, [], s); + c = t(c,:)'; + s = s*m; + r = repmat((1:s)',[1 k]); + q = accumarray([r(:) c(:)], i, [s n]); + p = repmat(p, [m 1]) + q; + l = l - k; +end +end + +end % uperm \ No newline at end of file diff --git a/tests/Makefile.am b/tests/Makefile.am index 92fabc336..b278925dc 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -190,6 +190,7 @@ MODFILES = \ identification/as2007/as2007_kronflags.mod \ identification/as2007/as2007_QT.mod \ identification/as2007/as2007_QT_equal_autocorr.mod \ + identification/as2007/as2007_order_1_2_3.mod \ identification/BrockMirman/BrockMirman.mod \ identification/cgg/cgg_criteria_differ.mod \ identification/ident_unit_root/ident_unit_root.mod \ diff --git a/tests/identification/as2007/as2007_QT.mod b/tests/identification/as2007/as2007_QT.mod index 51318e1cc..1837bd236 100644 --- a/tests/identification/as2007/as2007_QT.mod +++ b/tests/identification/as2007/as2007_QT.mod @@ -76,14 +76,14 @@ identification(parameter_set=calibration, load('G_QT'); %note that this is computed using replication files of Qu and Tkachenko (2012) temp = load([M_.dname filesep 'identification' filesep M_.fname '_identif']); -G_dynare = temp.ide_spectrum_point.G; +G_dynare = temp.ide_spectrum_point.dSPECTRUM; % Compare signs if ~isequal(sign(G_dynare),sign(G_QT)) error('signs of normalized G are note equal'); end % Compare normalized versions -tilda_G_dynare = temp.ide_spectrum_point.tilda_G; +tilda_G_dynare = temp.ide_spectrum_point.tilda_dSPECTRUM; ind_G_QT = (find(max(abs(G_QT'),[],1) > temp.store_options_ident.tol_deriv)); tilda_G_QT = zeros(size(G_QT)); delta_G_QT = sqrt(diag(G_QT(ind_G_QT,ind_G_QT))); diff --git a/tests/identification/as2007/as2007_order_1_2_3.mod b/tests/identification/as2007/as2007_order_1_2_3.mod new file mode 100644 index 000000000..a8fcbfd2d --- /dev/null +++ b/tests/identification/as2007/as2007_order_1_2_3.mod @@ -0,0 +1,109 @@ +%this is the mod file used in replication files of An and Schorfheide (2007) +% modified to include some obvious and artificial identification failures +% and to check whether all kronflags are working +% created by Willi Mutschler (willi@mutschler.eu) + +var y R g z c dy p YGR INFL INT; +varobs y R g z c dy p YGR INFL INT; +varexo e_r e_g e_z; +parameters sigr sigg sigz tau phi psi1 psi2 rhor rhog rhoz rrst pist gamst nu cyst dumpy dumpyrhog; + +rrst = 1.0000; +pist = 3.2000; +gamst= 0.5500; +tau = 2.0000; +nu = 0.1000; +kap = 0.3300; +phi = tau*(1-nu)/nu/kap/exp(pist/400)^2; +cyst = 0.8500; +psi1 = 1.5000; +psi2 = 0.1250; +rhor = 0.7500; +rhog = 0.9500; +rhoz = 0.9000; +sigr = 0.2; +sigg = 0.6; +sigz = 0.3; +dumpy = 0; +dumpyrhog = 1; + +model; +#pist2 = exp(pist/400); +#rrst2 = exp(rrst/400); +#bet = 1/rrst2; +#gst = 1/cyst; +#cst = (1-nu)^(1/tau); +#yst = cst*gst; +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+dy(+1)+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) + sigr*e_r; +g = dumpyrhog*rhog*g(-1) + sigg*e_g; +z = rhoz*z(-1) + sigz*e_z; +YGR = gamst+100*(dy+z); +INFL = pist+400*p; +INT = pist+rrst+4*gamst+400*R; +dy = y - y(-1); +end; + +shocks; +var e_r = 0.6^2; +var e_g = 0.5^2; +var e_z = 0.4^2; +corr e_r, e_g = 0.3; +corr e_r, e_z = 0.2; +corr e_z, e_g = 0.1; +end; + +steady_state_model; +z=0; g=0; c=0; y=0; p=0; R=0; dy=0; +YGR=gamst; INFL=pist; INT=pist+rrst+4*gamst; +end; + +estimated_params; +tau, 2, 1e-5, 10, gamma_pdf, 2, 0.5; + +%these parameters do not enter the linearized solution +cyst, 0.85, 1e-5, 0.99999, beta_pdf, 0.85, 0.1; +sigg, 0.6, 1e-8, 5, inv_gamma_pdf, 0.4, 4; +rhoz, 0.9, 1e-5, 0.99999, beta_pdf, 0.66, 0.15; +corr e_r,e_g, 0.3, 1e-8, 5, inv_gamma_pdf, 0.4, 4; +corr e_z,e_g, 0.3, 1e-8, 5, inv_gamma_pdf, 0.4, 4; +corr e_z,e_r, 0.3, 1e-8, 5, inv_gamma_pdf, 0.4, 4; + +%these parameters could only be identified from the steady state of YGR INFL and INT, however, we observer y pi R instead +rrst, 1, 1e-5, 10, gamma_pdf, 0.8, 0.5; +gamst, 0.55, -5, 5, normal_pdf, 0.4, 0.2; +dumpy, 0, -10, 10, normal_pdf, 0, 1; + +%these parameters jointly determine the slope kappa of the linearized new keynesian phillips curve +pist, 3.2, 1e-5, 20, gamma_pdf, 4, 2; +nu, 0.1, 1e-5, 0.99999, beta_pdf, 0.1, .05; +phi, 50, 1e-5, 100, gamma_pdf, 50, 20; + +%these parameters are pairwise collinear as one should not use both formulations for the standard error of a shock +sigz, 0.3, 1e-8, 5, inv_gamma_pdf, 0.4, 4; +stderr e_z, 0.3, 1e-8, 5, inv_gamma_pdf, 0.4, 4; + +%these parameters are pairwise collinear as they are multiplicative +rhog, 0.95, 1e-5, 0.99999, beta_pdf, 0.8, 0.1; +dumpyrhog, 1, -10, 10, normal_pdf, 1, 1; + +%these parameters are jointly not identified due to the specification of the Taylor rule +psi1, 1.5, 1e-5, 10, gamma_pdf, 1.5, 0.25; +psi2, 0.125, 1e-5, 10, gamma_pdf, 0.5, 0.25; +rhor, 0.75, 1e-5, 0.99999, beta_pdf, 0.5, 0.2; +stderr e_r, 0.2, 1e-8, 5, inv_gamma_pdf, 0.3, 4; + +end; + +steady; +check; +stoch_simul(order=3,irf=0,periods=0); %needed for identification(order=3) + +@#for ORDER in [1, 2, 3] +@#for KRONFLAG in [-1, -2, 0] +fprintf('*** ORDER = @{ORDER} WITH ANALYTIC_DERIVATION_MODE=@{KRONFLAG} ***\n') +identification(order=@{ORDER}, parameter_set=calibration, grid_nbr=10,analytic_derivation_mode=@{KRONFLAG}); +@#endfor +@#endfor \ No newline at end of file diff --git a/tests/identification/rbc_ident/rbc_ident_varexo_only.mod b/tests/identification/rbc_ident/rbc_ident_varexo_only.mod index f66beb54c..53d210194 100644 --- a/tests/identification/rbc_ident/rbc_ident_varexo_only.mod +++ b/tests/identification/rbc_ident/rbc_ident_varexo_only.mod @@ -96,7 +96,7 @@ identification(advanced=1); temp=load([M_.dname filesep 'identification' filesep M_.fname '_ML_Starting_value_identif']) temp_comparison=load(['rbc_ident_std_as_structural_par' filesep 'identification' filesep 'rbc_ident_std_as_structural_par' '_ML_Starting_value_identif']) -if max(abs(temp.ide_hess_point.ide_strength_J - temp_comparison.ide_hess_point.ide_strength_J))>1e-8 ||... +if max(abs(temp.ide_hess_point.ide_strength_dMOMENTS - temp_comparison.ide_hess_point.ide_strength_dMOMENTS))>1e-8 ||... max(abs(temp.ide_hess_point.deltaM - temp_comparison.ide_hess_point.deltaM))>1e-8 ||... max(abs(temp.ide_moments_point.MOMENTS- temp_comparison.ide_moments_point.MOMENTS))>1e-8 ||... max(abs(temp.ide_moments_point.MOMENTS- temp_comparison.ide_moments_point.MOMENTS))>1e-8