diff --git a/matlab/estimation/likelihood_quadratic_approximation.m b/matlab/estimation/likelihood_quadratic_approximation.m
new file mode 100644
index 000000000..ac281596c
--- /dev/null
+++ b/matlab/estimation/likelihood_quadratic_approximation.m
@@ -0,0 +1,141 @@
+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 function 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);
+        size(xx)
+        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
+
+end