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