1function mode_check(fun,x,hessian_mat,DynareDataset,DatasetInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,BoundsInfo,DynareResults)
2% Checks the estimated ML mode or Posterior mode.
3
4%@info:
5%! @deftypefn {Function File} mode_check (@var{fun}, @var{x}, @var{hessian_mat}, @var{DynareDataset}, @var{DynareOptions}, @var{Model}, @var{EstimatedParameters}, @var{BayesInfo}, @var{DynareResults})
6%! @anchor{mode_check}
7%! @sp 1
8%! Checks the estimated ML mode or Posterior mode by plotting sections of the likelihood/posterior kernel.
9%! Each plot shows the variation of the likelihood implied by the variations of a single parameter, ceteris paribus)
10%! @sp 2
11%! @strong{Inputs}
12%! @sp 1
13%! @table @ @var
14%! @item fun
15%! Objective function.
16%! @item x
17%! Estimated mode.
18%! @item start
19%! Hessian of the objective function at the estimated mode @var{x}.
20%! @item DynareDataset
21%! Structure specifying the dataset used for estimation (dataset_).
22%! @item DynareOptions
23%! Structure defining dynare's options (options_).
24%! @item Model
25%! Structure specifying the (estimated) model (M_).
26%! @item EstimatedParameters
27%! Structure specifying the estimated parameters (estimated_params_).
28%! @item BayesInfo
29%! Structure containing information about the priors used for estimation (bayestopt_).
30%! @item DynareResults
31%! Structure gathering the results (oo_).
32%! @end table
33%! @sp 2
34%! @strong{Outputs}
35%! @sp 2
36%! @strong{This function is called by:}
37%! @sp 2
38%! @strong{This function calls:}
39%! The objective function (@var{func}).
40%! @end deftypefn
41%@eod:
42
43% Copyright (C) 2003-2017 Dynare Team
44%
45% This file is part of Dynare.
46%
47% Dynare is free software: you can redistribute it and/or modify
48% it under the terms of the GNU General Public License as published by
49% the Free Software Foundation, either version 3 of the License, or
50% (at your option) any later version.
51%
52% Dynare is distributed in the hope that it will be useful,
53% but WITHOUT ANY WARRANTY; without even the implied warranty of
54% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
55% GNU General Public License for more details.
56%
57% You should have received a copy of the GNU General Public License
58% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
59
60TeX = DynareOptions.TeX;
61if ~isempty(hessian_mat)
62    [ s_min, k ] = min(diag(hessian_mat));
63end
64
65fval = feval(fun,x,DynareDataset,DatasetInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,BoundsInfo,DynareResults);
66
67if ~isempty(hessian_mat)
68    skipline()
69    disp('MODE CHECK')
70    skipline()
71    fprintf('Fval obtained by the minimization routine (minus the posterior/likelihood)): %f', fval);
72    skipline()
73    if s_min<eps
74        disp(sprintf('Most negative variance %f for parameter %d (%s = %f)', s_min, k , BayesInfo.name{k}, x(k)))
75    end
76end
77
78[nbplt,nr,nc,lr,lc,nstar] = pltorg(length(x));
79
80if ~exist([Model.fname filesep 'graphs'],'dir')
81    mkdir(Model.fname,'graphs');
82end
83if TeX && any(strcmp('eps',cellstr(DynareOptions.graph_format)))
84    fidTeX = fopen([Model.fname, '/graphs/', Model.fname '_CheckPlots.tex'],'w');
85    fprintf(fidTeX,'%% TeX eps-loader file generated by mode_check.m (Dynare).\n');
86    fprintf(fidTeX,['%% ' datestr(now,0) '\n']);
87    fprintf(fidTeX,' \n');
88end
89
90ll = DynareOptions.mode_check.neighbourhood_size;
91if isinf(ll)
92    DynareOptions.mode_check.symmetric_plots = false;
93end
94
95mcheck = struct('cross',struct(),'emode',struct());
96
97for plt = 1:nbplt
98    if TeX
99        NAMES = [];
100        TeXNAMES = [];
101    end
102    hh = dyn_figure(DynareOptions.nodisplay,'Name','Mode check plots');
103    for k=1:min(nstar,length(x)-(plt-1)*nstar)
104        subplot(nr,nc,k)
105        kk = (plt-1)*nstar+k;
106        [name,texname] = get_the_name(kk,TeX,Model,EstimatedParameters,DynareOptions);
107        if TeX
108            if isempty(NAMES)
109                NAMES = name;
110                TeXNAMES = texname;
111            else
112                NAMES = char(NAMES,name);
113                TeXNAMES = char(TeXNAMES,texname);
114            end
115        end
116        xx = x;
117        if x(kk)~=0 || ~isinf(BoundsInfo.lb(kk)) || ~isinf(BoundsInfo.lb(kk))
118            l1 = max(BoundsInfo.lb(kk),(1-sign(x(kk))*ll)*x(kk)); m1 = 0; %lower bound
119            l2 = min(BoundsInfo.ub(kk),(1+sign(x(kk))*ll)*x(kk)); %upper bound
120        else
121            %size info for 0 parameter is missing, use prior standard
122            %deviation
123            upper_bound=BoundsInfo.lb(kk);
124            if isinf(upper_bound)
125                upper_bound=-1e-6*DynareOptions.huge_number;
126            end
127            lower_bound=BoundsInfo.ub(kk);
128            if isinf(lower_bound)
129                lower_bound=-1e-6*DynareOptions.huge_number;
130            end
131            l1 = max(lower_bound,-BayesInfo.p2(kk)); m1 = 0; %lower bound
132            l2 = min(upper_bound,BayesInfo.p2(kk)); %upper bound
133        end
134        binding_lower_bound=0;
135        binding_upper_bound=0;
136        if isequal(x(kk),BoundsInfo.lb(kk))
137            binding_lower_bound=1;
138            bound_value=BoundsInfo.lb(kk);
139        elseif isequal(x(kk),BoundsInfo.ub(kk))
140            binding_upper_bound=1;
141            bound_value=BoundsInfo.ub(kk);
142        end
143        if DynareOptions.mode_check.symmetric_plots && ~binding_lower_bound && ~binding_upper_bound
144            if l2<(1+ll)*x(kk) %test whether upper bound is too small due to prior binding
145                l1 = x(kk) - (l2-x(kk)); %adjust lower bound to become closer
146                m1 = 1;
147            end
148            if ~m1 && (l1>(1-ll)*x(kk)) && (x(kk)+(x(kk)-l1)<BoundsInfo.ub(kk)) % if lower bound was truncated and using difference from lower bound does not violate upper bound
149                l2 = x(kk) + (x(kk)-l1); %set upper bound to same distance as lower bound
150            end
151        end
152        z1 = l1:((x(kk)-l1)/(DynareOptions.mode_check.number_of_points/2)):x(kk);
153        z2 = x(kk):((l2-x(kk))/(DynareOptions.mode_check.number_of_points/2)):l2;
154        z  = union(z1,z2);
155        if ~DynareOptions.mode_check.nolik
156            y = zeros(length(z),2);
157            dy = priordens(xx,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4);
158        end
159        for i=1:length(z)
160            xx(kk) = z(i);
161            [fval, info, exit_flag] = feval(fun,xx,DynareDataset,DatasetInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,BoundsInfo,DynareResults);
162            if exit_flag
163                y(i,1) = fval;
164            else
165                y(i,1) = NaN;
166                if DynareOptions.debug
167                    fprintf('mode_check:: could not solve model for parameter %s at value %4.3f, error code: %u\n',name,z(i),info(1))
168                end
169            end
170            if ~DynareOptions.mode_check.nolik
171                lnprior = priordens(xx,BayesInfo.pshape,BayesInfo.p6,BayesInfo.p7,BayesInfo.p3,BayesInfo.p4);
172                y(i,2)  = (y(i,1)+lnprior-dy);
173            end
174        end
175        mcheck.cross = setfield(mcheck.cross, name, [transpose(z), -y]);
176        mcheck.emode = setfield(mcheck.emode, name, x(kk));
177        fighandle=plot(z,-y);
178        hold on
179        yl=get(gca,'ylim');
180        plot( [x(kk) x(kk)], yl, 'c', 'LineWidth', 1)
181        NaN_index = find(isnan(y(:,1)));
182        zNaN = z(NaN_index);
183        yNaN = yl(1)*ones(size(NaN_index));
184        plot(zNaN,yNaN,'o','MarkerEdgeColor','r','MarkerFaceColor','r','MarkerSize',6);
185        title(name,'interpreter','none')
186        axis tight
187        if binding_lower_bound || binding_upper_bound
188            xl=get(gca,'xlim');
189            plot( [bound_value bound_value], yl, 'r--', 'LineWidth', 1)
190            xlim([xl(1)-0.5*binding_lower_bound*(xl(2)-xl(1)) xl(2)+0.5*binding_upper_bound*(xl(2)-xl(1))])
191        end
192        hold off
193        drawnow
194    end
195    if ~DynareOptions.mode_check.nolik
196        if isoctave
197            axes('outerposition',[0.3 0.93 0.42 0.07],'box','on'),
198        else
199            axes('position',[0.3 0.01 0.42 0.05],'box','on'),
200        end
201        line_color=get(fighandle,'color');
202        plot([0.48 0.68],[0.5 0.5],'color',line_color{2})
203        hold on, plot([0.04 0.24],[0.5 0.5],'color',line_color{1})
204        set(gca,'xlim',[0 1],'ylim',[0 1],'xtick',[],'ytick',[])
205        text(0.25,0.5,'log-post')
206        text(0.69,0.5,'log-lik kernel')
207    end
208    dyn_saveas(hh,[Model.fname, '/graphs/', Model.fname '_CheckPlots' int2str(plt) ],DynareOptions.nodisplay,DynareOptions.graph_format);
209    if TeX && any(strcmp('eps',cellstr(DynareOptions.graph_format)))
210        % TeX eps loader file
211        fprintf(fidTeX,'\\begin{figure}[H]\n');
212        for jj = 1:min(nstar,length(x)-(plt-1)*nstar)
213            fprintf(fidTeX,'\\psfrag{%s}[1][][0.5][0]{%s}\n',deblank(NAMES(jj,:)),deblank(TeXNAMES(jj,:)));
214        end
215        fprintf(fidTeX,'\\centering \n');
216        fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s_CheckPlots%s}\n',DynareOptions.figures.textwidth*min(k/nc,1),[Model.fname, '/graphs/',Model.fname],int2str(plt));
217        fprintf(fidTeX,'\\caption{Check plots.}');
218        fprintf(fidTeX,'\\label{Fig:CheckPlots:%s}\n',int2str(plt));
219        fprintf(fidTeX,'\\end{figure}\n');
220        fprintf(fidTeX,' \n');
221    end
222end
223if TeX && any(strcmp('eps',cellstr(DynareOptions.graph_format)))
224    fclose(fidTeX);
225end
226
227OutputDirectoryName = CheckPath('modecheck',Model.dname);
228save([OutputDirectoryName '/check_plot_data.mat'],'mcheck');
229