Fix header doc.

matlab/estimation/likelihood_quadratic_approximation.m

@@ -0,0 +1,191 @@
+function [f, df, d2f, R2] = likelihood_quadratic_approximation(particles, likelihoodvalues)
+
+% Approximate the shape of the likelihood function with a multivariate second order polynomial.
+%
+%
+% INPUTS
+% - particles            [double]   n×p matrix of (p) particles around the estimated posterior mode.
+% - likelihoodvalues     [double]   p×1 vector of corresponding values for the likelihood (or posterior kernel).
+%
+% OUTPUTS
+% - f                    [handle]   function handle for the approximated likelihood.
+% - df                   [handle]   function handle for the gradient of the approximated likelihood.
+% - d2f                  [handle]   Hessian matrix of the approximated likelihood (constant since we consider a second order multivariate polynomial)
+% - R2                   [double]   scalar, goodness of fit measure.
+%
+% REMARKS
+% [1] Function f takes a n×m matrix as input argument (the approximated likelihood is evaluated in m points) and returns a m×1 vector.
+% [2] Funtion df takes a n×1 vector as input argument (the point where the gradient is computed) and returns a n×1 vector.
+
+% Copyright © 2024 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
+
+    n = rows(particles);      % Number of parmaeters
+    p = columns(particles);   % Number of particles
+    q = 1 + n + n*(n+1)/2;    % Number of regressors (with a constant)
+
+    if p<=q
+        error('Quadratic approximation requires more than %u particles.', q)
+    end
+
+    %
+    % Build the set of regressors.
+    %
+
+    X = NaN(p, q);
+    X(:,1) = 1;                                   % zero order term
+    X(:,2:n+1) = transpose(particles);            % first order terms
+    X(:,n+2:end) = crossproducts(particles);      % second order terms
+
+    %
+    % Perform the regression
+    %
+
+    parameters = X\likelihoodvalues(:);
+
+    %
+    % Return a function to evaluate the approximation at x (a n×1 vector).
+    %
+
+    f = @(X) parameters(1) + transpose(X)*parameters(2:n+1) + crossproducts(X)*parameters(n+2:end);
+
+    if nargout>1
+        %
+        % Return a function to evaluate the gradient of the approximation at x (a n×1 vector)
+        %
+
+        df = @(X) parameters(2:n+1) + dcrossproducts(X)*parameters(n+2:end);
+
+        if nargout>2
+            %
+            % Return the hessian matrix of the approximation.
+            %
+
+            d2f = NaN(n,n);
+
+            h = 1;
+            for i=1:n
+                for j=i:n
+                    d2f(i,j) = parameters(n+1+h);
+                    if ~isequal(j, i)
+                        d2f(j,i) = d2f(i,j);
+                    end
+                    h = h+1;
+                end
+            end
+
+            if nargout>3
+                %
+                % Return a measure of fit goodness
+                %
+
+                R2 = 1-sum((likelihoodvalues(:)-X*parameters).^2)/sum(demean(likelihoodvalues(:)).^2);
+
+            end
+
+        end
+    end
+
+
+    function XX = crossproducts(X)
+    %                         n          n
+    % XX*ones(1,(n+1)*n/2) =  ∑  xᵢ² + 2 ∑ xᵢxⱼ
+    %                        i=1        i=1
+    %                                   j>i
+        XX = NaN(columns(X), n*(n+1)/2);
+        column = 1;
+        for i=1:n
+            XX(:,column) = transpose(X(i,:).*X(i,:));
+            column = column+1;
+            for j=i+1:n
+                XX(:,column) = 2*transpose(X(i,:).*X(j,:));
+                column = column+1;
+            end
+        end
+    end
+
+
+    function xx = dcrossproducts(x)
+        xx = zeros(n, n*(n+1)/2);
+        for i = 1:n
+            base = (i-1)*n-sum(0:i-2);
+            incol = 1;
+            xx(i,base+incol) = 2*x(i);
+            for j = i+1:n
+                incol = incol+1;
+                xx(i,incol) = 2*x(j);
+            end
+            for j=1:i-1
+                base = (j-1)*n-sum(0:j-2)+1;
+                colid = base+i-j;
+                xx(i,colid) = 2*x(j);
+            end
+        end
+    end
+
+    return % --*-- Unit tests --*--
+
+    %@test:1
+    % Create data
+    X = randn(10,1000);
+    Y = 1 + rand(1,10)*X.^2+0.01*randn(1,1000);
+
+    % Perform approximation
+    try
+        [f,df, d2f,R2] = likelihood_quadratic_approximation(X,Y);
+        t(1) = true;
+    catch
+        t(1) = false;
+    end
+
+    % Test returned arguments
+    if t(1)
+        try
+            y = f(randn(10,100));
+            t(2) = true;
+            if ~(rows(y)==100 && columns(y)==1)
+                t(2) = false;
+            end
+        catch
+            t(2) = false;
+        end
+        try
+            dy = df(zeros(10,1));
+            t(3) = true;
+            if ~(rows(dy)==10 && columns(dy)==1)
+                t(3) = false;
+            end
+        catch
+            t(3) = false;
+        end
+        t(4) = true;
+        if ~(rows(d2f)==10 && columns(d2f)==10)
+            t(4) = false
+        end
+        if ~(rows(d2f)==10 && columns(d2f)==10)
+            t(4) = false
+        end
+        t(4) = issymmetric(d2f);
+        t(4) = ispd(d2f);
+        t(5) = isscalar(R2);
+        t(5) = (R2>0) & (R2<1); % Note that in a nonlinear model nothing ensures that these inequalities are satisfied.
+    end
+
+    T = all(t);
+    %@eof:1
+
+end

matlab/utilities/general/ispd.m

@@ -1,30 +1,15 @@
 function [test, penalty] = ispd(A)
 
-%@info:
-%! @deftypefn {Function File} {[@var{test}, @var{penalty}  =} ispd (@var{A})
-%! @anchor{ispd}
-%! @sp 1
-%! Tests if the square matrix @var{A} is positive definite.
-%! @sp 2
-%! @strong{Inputs}
-%! @sp 1
-%! @table @ @var
-%! @item A
-%! A square matrix.
-%! @end table
-%! @sp 2
-%! @strong{Outputs}
-%! @sp 1
-%! @table @ @var
-%! @item test
-%! Integer scalar equal to 1 if @var{A} is a positive definite sqquare matrix, 0 otherwise.
-%! @item penalty
-%! Absolute value of the uum of the negative eigenvalues of A. This output argument is optional.
-%! @end table
-%! @end deftypefn
-%@eod:
+% Test if the square matrix A is positive definite.
+%
+% INPUTS
+% - A       [double]   n×n matrix.
+%
+% OUTPUTS
+% - test    [logical]  scalar, true if and only if matrix A is positive definite (and symmetric...)
+% - penalty [double]   scalar, absolute value of the uum of the negative eigenvalues of A. This output argument is optional.
 
-% Copyright © 2007-2022 Dynare Team
+% Copyright © 2007-2024 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -54,6 +39,7 @@ if nargout>1
         if isoctave && any(any(~isfinite(A))) % workaround for https://savannah.gnu.org/bugs/index.php?63082
             penalty = 1;
         else
+            % TODO: the penalty is only concerned with negative eigenvalues, we should also consider the case of non symmetric matrix A.
             a = diag(eig(A));
             k = find(a<0);
             if k > 0
@@ -61,4 +47,4 @@ if nargout>1
             end
         end
     end
-end
+end