Finished identification order=1|2|3

Note that I still need to do a code clean up (provide some licenses for functions from other people) and to double check order=3. There is also much room for speed and memory improvement, but the code works fine for now. I will also provide more information to the merge request soon about the detailed changes for future reference.

matlab/Q6_plication.m

+% 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

matlab/allVL1.m

+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

matlab/bivmom.m

+%
+% 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
+

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

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

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

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);

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 <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

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)';

matlab/prodmom.m

+%
+% 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;