Commit 24a2aa2a authored by Sébastien Villemot's avatar Sébastien Villemot

Merge branch 'identification_higher_order' into 'master'

Finished identification order=1|2|3

See merge request !1689
parents d40b7752 45e9771e
% By Willi Mutschler, September 26, 2016. Email: willi@mutschler.eu
% Quadruplication Matrix as defined by
% Meijer (2005) - Matrix algebra for higher order moments. Linear Algebra and its Applications, 410,pp. 112134
%
% 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
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
%
% 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
......@@ -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
......
......@@ -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
This diff is collapsed.
This diff is collapsed.
......@@ -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);
......
This diff is collapsed.
......@@ -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 <http://www.gnu.org/licenses/>.
% =========================================================================
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
......
This diff is collapsed.
......@@ -88,33 +88,53 @@ else
end
%% compute Kalman transition matrices and steady state with updated parameters
[A, B, ~, ~, M, options, oo] = dynare_resolve(M, options, oo);
ys = oo.dr.ys; %steady state of model variables in declaration order
y0 = ys(oo.dr.order_var); %steady state of model variables in DR order
[~,info,M,options,oo] = resol(0,M,options,oo);
options = rmfield(options,'options_ident');
oo.dr = pruned_state_space_system(M, options, oo.dr);
A = oo.dr.pruned.A;
B = oo.dr.pruned.B;
C = oo.dr.pruned.C(indvar,:);
D = oo.dr.pruned.D(indvar,:);
Om_z = oo.dr.pruned.Om_z;
Om_y = oo.dr.pruned.Om_y(indvar,indvar);
Varinov = oo.dr.pruned.Varinov;
obs_nbr = size(C,1);
%% out = [vech(cov(Y_t,Y_t)); vec(cov(Y_t,Y_{t-1}); ...; vec(cov(Y_t,Y_{t-nlags})] of indvar variables, in DR order. This is Iskrev (2010)'s J matrix.
if outputflag == 1
% Denote Ezz0 = E_t(z_t * z_t'), then the following Lyapunov equation defines the autocovariagrom: Ezz0 -A*Ezz*A' = B*Sig_e*B'
Ezz0 = lyapunov_symm(A,B*M.Sigma_e*B',options.lyapunov_fixed_point_tol,options.qz_criterium,options.lyapunov_complex_threshold,[],options.debug);
indzeros = find(abs(Ezz0) < 1e-12); %set small values to zero
Ezz0(indzeros) = 0;
[Ezz0,u] = lyapunov_symm(A, Om_z, options.lyapunov_fixed_point_tol, options.qz_criterium, options.lyapunov_complex_threshold, 1, options.debug);
stationary_vars = (1:size(C,1))';
if ~isempty(u)
x = abs(C*u);
stationary_vars = find(all(x < options.Schur_vec_tol,2));
end
Eyy0 = NaN*ones(obs_nbr,obs_nbr);
Eyy0(stationary_vars,stationary_vars) = C(stationary_vars,:)*Ezz0*C(stationary_vars,:)' + Om_y(stationary_vars,stationary_vars);
indzeros = find(abs(Eyy0) < 1e-12); %find values that are numerical zero
Eyy0(indzeros) = 0;
if useautocorr
sy = sqrt(diag(Ezz0));
sy = sy*sy';
sy = sqrt(diag(Ezz0)); %theoretical standard deviation
sy = sy(stationary_vars);
sy = sy*sy'; %cross products of standard deviations
sy0 = sy-diag(diag(sy))+eye(length(sy));
Ezz0corr = Ezz0./sy0;
out = dyn_vech(Ezz0corr(indvar,indvar)); %focus only on unique terms
Eyy0corr = NaN*ones(size(C,1),size(C,1));
Eyy0corr(stationary_vars,stationary_vars) = Eyy0./sy0;
out = dyn_vech(Eyy0corr); %focus only on unique terms
else
out = dyn_vech(Ezz0(indvar,indvar)); %focus only on unique terms
out = dyn_vech(Eyy0); %focus only on unique terms
end
% compute autocovariances/autocorrelations of lagged observed variables
for ii = 1:nlags
Ezzii = A^(ii)*Ezz0;
tmpEyyi = A*Ezz0*C(stationary_vars,:)' + B*Varinov*D(stationary_vars,:)';
Ai = eye(size(A,1)); %this is A^0
for i = 1:nlags
Eyyi = NaN*ones(obs_nbr,obs_nbr);
Eyyi(stationary_vars,stationary_vars) = C(stationary_vars,:)*Ai*tmpEyyi;
if useautocorr
Ezzii = Ezzii./sy;
Eyyi = Eyyi./sy;
end
out = [out;vec(Ezzii(indvar,indvar))];
end
out = [out;vec(Eyyi)];
Ai = Ai*A; %note that this is A^(i-1)
end
end
%% out = vec(g_omega). This is needed for Qu and Tkachenko (2012)'s G matrix.
......@@ -122,15 +142,13 @@ if outputflag == 2
% This computes the spectral density g_omega where the interval [-pi;\pi] is discretized by grid_nbr points
freqs = (0 : pi/(grid_nbr/2):pi);% we focus only on positive values including the 0 frequency
tpos = exp( sqrt(-1)*freqs); %Fourier frequencies
C = A(indvar,:);
D = B(indvar,:);
IA = eye(size(A,1));
var_nbr = length(indvar);
var_nbr = size(C,1);
out = zeros(var_nbr^2*length(freqs),1);
kk = 0;
for ig = 1:length(freqs)
Transferfct = D + C*((tpos(ig)*IA-A)\B);
g_omega = (1/(2*pi))*(Transferfct*M.Sigma_e*Transferfct'); % note that ' is the conjugate transpose
g_omega = (1/(2*pi))*(Transferfct*Varinov*Transferfct'); % note that ' is the conjugate transpose
kk = kk+1;
out(1 + (kk-1)*var_nbr^2 : kk*var_nbr^2) = g_omega(:);
end
......
......@@ -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);
%