1function [r,flag] = smm_objective(xparams,sample_moments,weighting_matrix,options,parallel)
2% Evaluates the objective of the Simulated Moments Method.
3%
4% INPUTS:
5%  xparams          [double]  p*1 vector of estimated parameters.
6%  sample_moments   [double]  n*1 vector of sample moments (n>=p).
7%  weighting_matrix [double]  n*n symetric, positive definite matrix.
8%  options          [      ]  Structure defining options for SMM.
9%  parallel         [      ]  Structure defining the parallel mode settings (optional).
10%
11% OUTPUTS:
12%  r                [double]  scalar, the value of the objective function.
13%  junk             [      ]  empty matrix.
14%
15% SPECIAL REQUIREMENTS
16%  The user has to provide a file where the moment conditions are defined.
17
18% Copyright (C) 2010-2019 Dynare Team
19%
20% This file is part of Dynare.
21%
22% Dynare is free software: you can redistribute it and/or modify
23% it under the terms of the GNU General Public License as published by
24% the Free Software Foundation, either version 3 of the License, or
25% (at your option) any later version.
26%
27% Dynare is distributed in the hope that it will be useful,
28% but WITHOUT ANY WARRANTY; without even the implied warranty of
29% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
30% GNU General Public License for more details.
31%
32% You should have received a copy of the GNU General Public License
33% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
34
35global M_ options_ oo_
36persistent mainStream mainState
37persistent priorObjectiveValue
38
39flag = 1;
40
41if nargin<5
42    if isempty(mainStream)
43        if matlab_ver_less_than('7.12')
44            mainStream = RandStream.getDefaultStream;
45        else
46            mainStream = RandStream.getGlobalStream;
47        end
48        mainState  = mainStream.State;
49    else
50        mainStream.State = mainState;
51    end
52end
53
54if isempty(priorObjectiveValue)
55    priorObjectiveValue = Inf;
56end
57
58penalty = 0;
59for i=1:options.estimated_parameters.nb
60    if ~isnan(options.estimated_parameters.upper_bound(i)) && xparams(i)>options.estimated_parameters.upper_bound(i)
61        penalty = penalty + (xparams(i)-options.estimated_parameters.upper_bound(i))^2;
62    end
63    if ~isnan(options.estimated_parameters.lower_bound(i)) && xparams(i)<options.estimated_parameters.lower_bound(i)
64        penalty = penalty + (xparams(i)-options.estimated_parameters.lower_bound(i))^2;
65    end
66end
67
68if penalty>0
69    flag = 0;
70    r = priorObjectiveValue + penalty;
71    return
72end
73
74save('estimated_parameters.mat','xparams');
75
76% Check for local determinacy of the deterministic steady state.
77noprint = options_.noprint; options_.noprint = 1;
78[~,local_determinacy_and_stability,info] = check(M_,options_,oo_); options_.noprint = noprint;
79if ~local_determinacy_and_stability
80    r = priorObjectiveValue * (1+info(2));
81    flag = 0;
82    return
83end
84
85simulated_moments = zeros(size(sample_moments));
86
87% Just to be sure that things don't mess up with persistent variables...
88clear perfect_foresight_simulation;
89
90if nargin<5
91    for s = 1:options.number_of_simulated_sample
92        time_series = extended_path([],options.simulated_sample_size,1);
93        data = time_series(options.observed_variables_idx,options.burn_in_periods+1:options.simulated_sample_size);
94        eval(['tmp = ' options.moments_file_name '(data);'])
95        simulated_moments = simulated_moments + tmp;
96        simulated_moments = simulated_moments / options.number_of_simulated_sample;
97    end
98else% parallel mode.
99    if ~isunix
100        error('The parallel version of SMM estimation is not implemented for non unix platforms!')
101    end
102    job_number = 1;% Remark. First job is for the master.
103    [~,hostname] = unix('hostname --fqdn');
104    hostname = deblank(hostname);
105    for i=1:length(parallel)
106        machine = deblank(parallel(i).machine);
107        if ~strcmpi(hostname,machine)
108            % For the slaves on a remote computer.
109            unix(['scp estimated_parameters.mat ' , parallel(i).login , '@' , machine , ':' parallel(i).folder ' > /dev/null']);
110        else
111            if ~strcmpi(pwd,parallel(i).folder)
112                % For the slaves on this computer but not in the same directory as the master.
113                unix(['cp estimated_parameters.mat ' , parallel(i).folder]);
114            end
115        end
116        for j=1:parallel(i).number_of_jobs
117            if (strcmpi(hostname,machine) && j>1) || ~strcmpi(hostname,machine)
118                job_number = job_number + 1;
119                unix(['ssh -A ' parallel(i).login '@' machine ' ./call_matlab_session.sh job' int2str(job_number) '.m &']);
120            end
121        end
122    end
123    % Finally the Master do its job
124    tStartMasterJob = clock;
125    eval('job1;')
126    tElapsedMasterJob = etime(clock, tStartMasterJob);
127    TimeLimit = tElapsedMasterJob*1.2;
128    % Master waits for the  slaves' output...
129    tStart = clock;
130    tElapsed = 0;
131    while tElapsed<TimeLimit
132        if ( length(dir('./intermediary_results_from_master_and_slaves/simulated_moments_slave_*.dat'))==job_number )
133            break
134        end
135        tElapsed = etime(clock, tStart);
136    end
137    try
138        tmp = zeros(length(sample_moments),1);
139        for i=1:job_number
140            simulated_moments = load(['./intermediary_results_from_master_and_slaves/simulated_moments_slave_' int2str(i) '.dat'],'-ascii');
141            tmp = tmp + simulated_moments;
142        end
143        simulated_moments = tmp / job_number;
144    catch
145        r = priorObjectiveValue*1.1;
146        flag = 0;
147        return
148    end
149end
150
151r = transpose(simulated_moments-sample_moments)*weighting_matrix*(simulated_moments-sample_moments);
152priorObjectiveValue = r;
153
154if (options.optimization_routine>0) && exist('optimization_path.mat')
155    load('optimization_path.mat');
156    new_state = [ r; xparams];
157    estimated_parameters_optimization_path = [ estimated_parameters_optimization_path , new_state ];
158    save('optimization_path.mat','estimated_parameters_optimization_path');
159end