[WIP] Approximation of the likelihood with a second order multivariate polynomial.

[skip ci]

doc/manual/source/conf.py

@@ -1,6 +1,6 @@
 # -*- coding: utf-8 -*-
 
-# Copyright © 2018-2023 Dynare Team
+# Copyright © 2018-2024 Dynare Team
 #
 # This file is part of Dynare.
 #
@@ -34,7 +34,7 @@ html_static_path = ['_static']
 master_doc = 'index'
 
 project = u'Dynare'
-copyright = u'1996–2023 Dynare Team'
+copyright = u'1996–2024 Dynare Team'
 author = u'Dynare Team'
 
 add_function_parentheses = False

doc/manual/source/index.rst

@@ -26,7 +26,7 @@ The following people used to be members of the team:
 * Ferhat Mihoubi
 * George Perendia
 
-Copyright © 1996-2023, Dynare Team.
+Copyright © 1996-2024, Dynare Team.
 
 Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.3 or any later version published by the Free Software Foundation; with no Invariant Sections, no Front-Cover Texts, and no Back-Cover Texts.
 

license.txt

@@ -1,6 +1,6 @@
 Format: https://www.debian.org/doc/packaging-manuals/copyright-format/1.0/
 Upstream-Name: Dynare
-Upstream-Contact: Dynare Team, whose members in 2023 are:
+Upstream-Contact: Dynare Team, whose members in 2024 are:
                   - Stéphane Adjemian <stephane.adjemian@univ-lemans.fr>
                   - Michel Juillard <michel.juillard@mjui.fr>
                   - Frédéric Karamé <frederic.karame@univ-lemans.fr>
@@ -23,7 +23,7 @@ Upstream-Contact: Dynare Team, whose members in 2023 are:
 Source: https://www.dynare.org
 
 Files: *
-Copyright: 1996-2023 Dynare Team
+Copyright: 1996-2024 Dynare Team
 License: GPL-3+
 
 Files: matlab/+occbin/IVF_core.m

matlab/+bvar/forecast.m

@@ -119,7 +119,7 @@ OutputDirectoryName = CheckPath('graphs',M_.dname);
 dyn_graph=bvar.graph_init(sprintf('BVAR forecasts (nlags = %d)', nlags), ny, {'b-' 'g-' 'g-' 'r-' 'r-'});
 
 for i = 1:ny
-    dyn_graph=dynare_graph(dyn_graph,[ sims_no_shock_median(:, i) ...
+    dyn_graph=plot_graph(dyn_graph,[ sims_no_shock_median(:, i) ...
                         sims_no_shock_up_conf(:, i) sims_no_shock_down_conf(:, i) ...
                         sims_with_shocks_up_conf(:, i) sims_with_shocks_down_conf(:, i) ], ...
                            options_.varobs{i});
@@ -183,3 +183,31 @@ for i = 1:length(options_.varobs)
         oo_.bvar.forecast.rmse.(name) = rmse(i);
     end
 end
+
+
+function dyn_graph=plot_graph(dyn_graph,y,tit,x)
+% function plot_graph(dyn_graph, y,tit,x)
+
+if nargin < 4
+    x = (1:size(y,1))';
+end
+nplot = dyn_graph.plot_nbr + 1;
+if nplot > dyn_graph.max_nplot
+    figure('Name',dyn_graph.figure_name);
+    nplot = 1;
+end
+dyn_graph.plot_nbr = nplot;
+subplot(dyn_graph.nr,dyn_graph.nc,nplot);
+
+line_types = dyn_graph.line_types;
+line_type = line_types{1};
+for i=1:size(y,2)
+    if length(line_types) > 1
+        line_type = line_types{i};
+    end
+
+    plot(x,y(:,i),line_type);
+    hold on
+end
+title(tit);
+hold off

matlab/dynare_graph.m

@@ -1,54 +0,0 @@
-function dyn_graph=dynare_graph(dyn_graph,y,tit,x)
-% function dynare_graph(y,tit,x)
-% graphs
-%
-% INPUT
-%   figure_name: name of the figures
-%   colors: line colors
-%
-% OUTPUT
-%   dyn_graph: structure with figure information
-%
-% SPECIAL REQUIREMENT
-%   none
-
-% Copyright © 2006-2017 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/>.
-
-if nargin < 4
-    x = (1:size(y,1))';
-end
-nplot = dyn_graph.plot_nbr + 1;
-if nplot > dyn_graph.max_nplot
-    figure('Name',dyn_graph.figure_name);
-    nplot = 1;
-end
-dyn_graph.plot_nbr = nplot;
-subplot(dyn_graph.nr,dyn_graph.nc,nplot);
-
-line_types = dyn_graph.line_types;
-line_type = line_types{1};
-for i=1:size(y,2)
-    if length(line_types) > 1
-        line_type = line_types{i};
-    end
-
-    plot(x,y(:,i),line_type);
-    hold on
-end
-title(tit);
-hold off

matlab/dynare_graph_close.m

@@ -1,29 +0,0 @@
-function dynare_graph_close()
-% function dynare_graph_close()
-% close a figure
-%
-% INPUT
-%   none
-%
-% OUTPUT
-%   none
-%
-% SPECIAL REQUIREMENT
-%   none
-
-% Copyright © 2006-2017 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/>.

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

matlab/forecast_graphs.m

@@ -76,7 +76,7 @@ for j= 1:nvar
             fprintf(fidTeX,' \n');
         end
         n_fig =n_fig+1;
-        eval(['hh_fig=dyn_figure(options_.nodisplay,''Name'',''Forecasts (' int2str(n_fig) ')'');']);
+        hh_fig=dyn_figure(options_.nodisplay,'Name',['Forecasts (' int2str(n_fig) ')']);
         m = 1;
     end
     subplot(nr, nc, m);
@@ -138,7 +138,7 @@ if isfield(oo_.forecast,'HPDinf_ME')
                 fprintf(fidTeX,' \n');
             end
             n_fig =n_fig+1;
-            eval(['hh_fig=dyn_figure(options_.nodisplay,''Name'',''Forecasts (' int2str(n_fig) ')'');']);
+            hh_fig=dyn_figure(options_.nodisplay,'Name',['Forecasts including ME (' int2str(n_fig) ')']);
             m = 1;
         end
         subplot(nr,nc,m);

matlab/utilities/dataset/makedataset.m

@@ -242,8 +242,14 @@ if newdatainterface
     end
 else
     % ...  or check that nobs is smaller than the number of observations in dataset_.
-    if nobs>dataset_.nobs
-        error('makedataset: nobs (%s) cannot be greater than the last date in the dataset (%s)!', num2str(nobs), num2str(dataset_.nobs))
+    if FIRSTOBS>dataset_.dates(1)
+        if FIRSTOBS+nobs-1>dataset_.dates(end)
+            error('makedataset: Given first_obs=%u and %u total observations in the dataset, the current nobs of %s must not be greater than %s!', options_.first_obs, dataset_.nobs, num2str(nobs), num2str(dataset_.nobs-find(dataset_.dates==FIRSTOBS)+1))
+        end
+    else
+        if nobs>dataset_.nobs
+            error('makedataset: nobs (%s) cannot be greater than the last date in the dataset (%s)!', num2str(nobs), num2str(dataset_.nobs))
+        end
     end
 end
 

windows/deps/Makefile

@@ -157,6 +157,7 @@ mingw64: tarballs/mingw-w64-x86_64-gcc-$(MINGW64_GCC_VERSION)-any.pkg.tar.zst ta
 	touch $@
 
 tarballs/mingw-w64-x86_64-%-any.pkg.tar.zst:
+	mkdir -p tarballs
 	wget $(WGET_OPTIONS) -O $@ http://repo.msys2.org/mingw/x86_64/$(notdir $@)
 
 clean-msys2: